Accelerating seminumerical Fock-exchange calculations using mixed single- and double-precision arithmethic.

We investigate the applicability of single-precision (fp32) floating point operations within our linear-scaling, seminumerical exchange method sn-LinK [Laqua et al., J. Chem. Theory Comput. 16, 1456 (2020)] and find that the vast majority of the three-center-one-electron (3c1e) integrals can be computed with reduced numerical precision with virtually no loss in overall accuracy. This leads to a near doubling in performance on central processing units (CPUs) compared to pure fp64 evaluation. Since the cost of evaluating the 3c1e integrals is less significant on graphic processing units (GPUs) compared to CPU, the performance gains from accelerating 3c1e integrals alone is less impressive on GPUs. Therefore, we also investigate the possibility of employing only fp32 operations to evaluate the exchange matrix within the self-consistent-field (SCF) followed by an accurate one-shot evaluation of the exchange energy using mixed fp32/fp64 precision. This still provides very accurate (1.8 µEh maximal error) results while providing a sevenfold speedup on a typical "gaming" GPU (GTX 1080Ti). We also propose the use of incremental exchange-builds to further reduce these errors. The proposed SCF scheme (i-sn-LinK) requires only one mixed-precision exchange matrix calculation, while all other exchange-matrix builds are performed with only fp32 operations. Compared to pure fp64 evaluation, this leads to 4-7× speedups for the whole SCF procedure without any significant deterioration of the results or the convergence behavior.

Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package.

This article summarizes technical advances contained in the fifth major release of the Q-Chem quantum chemistry program package, covering developments since 2015. A comprehensive library of exchange-correlation functionals, along with a suite of correlated many-body methods, continues to be a hallmark of the Q-Chem software. The many-body methods include novel variants of both coupled-cluster and configuration-interaction approaches along with methods based on the algebraic diagrammatic construction and variational reduced density-matrix methods. Methods highlighted in Q-Chem 5 include a suite of tools for modeling core-level spectroscopy, methods for describing metastable resonances, methods for computing vibronic spectra, the nuclear-electronic orbital method, and several different energy decomposition analysis techniques. High-performance capabilities including multithreaded parallelism and support for calculations on graphics processing units are described. Q-Chem boasts a community of well over 100 active academic developers, and the continuing evolution of the software is supported by an "open teamware" model and an increasingly modular design.

A Fermi smearing variant of the Tamm-Dancoff approximation for nonadiabatic dynamics involving S1-S0 transitions: Validation and application to azobenzene.

The main shortcoming of time-dependent density functional theory (TDDFT) regarding its use for nonadiabatic molecular dynamics (NAMD) is its incapability to describe conical intersections involving the ground state. To overcome this problem, we combine Fermi smearing (FS) DFT with a fractional-occupation variant of the Tamm-Dancoff approximation (TDA) of TDDFT in the generalized gradient approximation. The resulting method (which we denote as FS-TDA) gives access to ground- and excited-state energies, gradients, and nonadiabatic coupling vectors, which are physically correct even in the vicinity of S1-S0 conical intersections. This is shown for azobenzene, a widely used photoswitch, via single point calculations and NAMD simulations of its cis-trans photoisomerization. We conclude that FS-TDA may be used as an efficient alternative to investigate these processes.

Identifying Free Energy Hot-Spots in Molecular Transformations.

The free energy is one of the central quantities in material and natural sciences. While being well-established, e.g., in drug design or catalyst optimization, computational methods lack a straightforward way to gain deeper insights into the calculated free energy, and thus the underlying chemical or physical processes. Here, we present a generally applicable, spectrum-based ansatz that tackles this shortcoming by identifying contributions from specific atoms or groups to the vibrational free energy. We illustrate this in studies of the bromodomain-inhibitor binding and the anomeric effect in glucose providing quantitative evidence in line with chemical intuition in both cases. For the latter example we also report an experimental infrared spectrum and find excellent agreement with our simulated spectra.

Calculating free energies from the vibrational density of states function: Validation and critical assessment.

We explore and show the usefulness of the density of states function for computing vibrational free energies and free energy differences between small systems. Therefore, we compare this density of states integration method (DSI) to more established schemes such as Bennett's Acceptance Ratio method (BAR), the Normal Mode Analysis (NMA), and the Quasiharmonic Analysis (QHA). The strengths and shortcomings of all methods are highlighted with three numerical examples. Furthermore, the free energy of the ionization of ammonia and the mutation from serine to cysteine are computed using extensive ab initio molecular dynamics simulations. We conclude that DSI improves upon the other frequency-based methods (NMA and QHA) regarding the treatment of anharmonicity and yielding results comparable to BAR in all cases without the need for alchemical transformations. Low-frequency modes lead to larger errors indicating that long simulation times might be required for larger systems. In addition, we introduce the use of DSI for the localization of the vibrational free energy to specific atoms or residues, leading to insights into the underlying process, a unique feature that is only offered by this method.

An improved molecular partitioning scheme for numerical quadratures in density functional theory.

We present a modification to Becke's molecular partitioning scheme [A. D. Becke, J. Chem. Phys. 88, 2547 (1988)] that provides substantially better accuracy for weakly bound complexes and allows for a faster and linear scaling grid generation without introducing a cutoff error. We present the accuracy of our new partitioning scheme for atomization energies of small molecules and for interaction energies of van der Waals complexes. Furthermore, the efficiency and scaling behavior of the grid generation are demonstrated for large molecular systems with up to 1707 atoms.

Communication: Density functional theory model for multi-reference systems based on the exact-exchange hole normalization.

The correct description of multi-reference electronic ground states within Kohn-Sham density functional theory (DFT) requires an ensemble-state representation, employing fractionally occupied orbitals. However, the use of fractional orbital occupation leads to non-normalized exact-exchange holes, resulting in large fractional-spin errors for conventional approximative density functionals. In this communication, we present a simple approach to directly include the exact-exchange-hole normalization into DFT. Compared to conventional functionals, our model strongly improves the description for multi-reference systems, while preserving the accuracy in the single-reference case. We analyze the performance of our proposed method at the example of spin-averaged atoms and spin-restricted bond dissociation energy surfaces.

Screening methods for linear-scaling short-range hybrid calculations on CPU and GPU architectures.

We present screening schemes that allow for efficient, linear-scaling short-range exchange calculations employing Gaussian basis sets for both CPU and GPU architectures. They are based on the LinK [C. Ochsenfeld et al., J. Chem. Phys. 109, 1663 (1998)] and PreLinK [J. Kussmann and C. Ochsenfeld, J. Chem. Phys. 138, 134114 (2013)] methods, but account for the decay introduced by the attenuated Coulomb operator in short-range hybrid density functionals. Furthermore, we discuss the implementation of short-range electron repulsion integrals on GPUs. The introduction of our screening methods allows for speedups of up to a factor 7.8 as compared to the underlying linear-scaling algorithm, while retaining full numerical control over the accuracy. With the increasing number of short-range hybrid functionals, our new schemes will allow for significant computational savings on CPU and GPU architectures.

Computation of indirect nuclear spin-spin couplings with reduced complexity in pure and hybrid density functional approximations.

We present a (sub)linear-scaling algorithm to determine indirect nuclear spin-spin coupling constants at the Hartree-Fock and Kohn-Sham density functional levels of theory. Employing efficient integral algorithms and sparse algebra routines, an overall (sub)linear scaling behavior can be obtained for systems with a non-vanishing HOMO-LUMO gap. Calculations on systems with over 1000 atoms and 20 000 basis functions illustrate the performance and accuracy of our reference implementation. Specifically, we demonstrate that linear algebra dominates the runtime of conventional algorithms for 10 000 basis functions and above. Attainable speedups of our method exceed 6 × in total runtime and 10 × in the linear algebra steps for the tested systems. Furthermore, a convergence study of spin-spin couplings of an aminopyrazole peptide upon inclusion of the water environment is presented: using the new method it is shown that large solvent spheres are necessary to converge spin-spin coupling values.

A reduced-scaling density matrix-based method for the computation of the vibrational Hessian matrix at the self-consistent field level.

An analytical method to calculate the molecular vibrational Hessian matrix at the self-consistent field level is presented. By analysis of the multipole expansions of the relevant derivatives of Coulomb-type two-electron integral contractions, we show that the effect of the perturbation on the electronic structure due to the displacement of nuclei decays at least as r(-2) instead of r(-1). The perturbation is asymptotically local, and the computation of the Hessian matrix can, in principle, be performed with O(N) complexity. Our implementation exhibits linear scaling in all time-determining steps, with some rapid but quadratic-complexity steps remaining. Sample calculations illustrate linear or near-linear scaling in the construction of the complete nuclear Hessian matrix for sparse systems. For more demanding systems, scaling is still considerably sub-quadratic to quadratic, depending on the density of the underlying electronic structure.

QM/MM Free Energy Calculations of Long-Range Biological Protonation Dynamics by Adaptive and Focused Sampling.

Water-mediated proton transfer reactions are central for catalytic processes in a wide range of biochemical systems, ranging from biological energy conversion to chemical transformations in the metabolism. Yet, the accurate computational treatment of such complex biochemical reactions is highly challenging and requires the application of multiscale methods, in particular hybrid quantum/classical (QM/MM) approaches combined with free energy simulations. Here, we combine the unique exploration power of new advanced sampling methods with density functional theory (DFT)-based QM/MM free energy methods for multiscale simulations of long-range protonation dynamics in biological systems. In this regard, we show that combining multiple walkers/well-tempered metadynamics with an extended system adaptive biasing force method (MWE) provides a powerful approach for exploration of water-mediated proton transfer reactions in complex biochemical systems. We compare and combine the MWE method also with QM/MM umbrella sampling and explore the sampling of the free energy landscape with both geometric (linear combination of proton transfer distances) and physical (center of excess charge) reaction coordinates and show how these affect the convergence of the potential of mean force (PMF) and the activation free energy. We find that the QM/MM-MWE method can efficiently explore both direct and water-mediated proton transfer pathways together with forward and reverse hole transfer mechanisms in the highly complex proton channel of respiratory Complex I, while the QM/MM-US approach shows a systematic convergence of selected long-range proton transfer pathways. In this regard, we show that the PMF along multiple proton transfer pathways is recovered by combining the strengths of both approaches in a QM/MM-MWE/focused US (FUS) scheme and reveals new mechanistic insight into the proton transfer principles of Complex I. Our findings provide a promising basis for the quantitative multiscale simulations of long-range proton transfer reactions in biological systems.

Pre-selective screening for matrix elements in linear-scaling exact exchange calculations.

We present a simple but accurate preselection method based on Schwarz integral estimates to determine the significant elements of the exact exchange matrix before its evaluation, thus providing an asymptotical linear-scaling behavior for non-metallic systems. Our screening procedure proves to be highly suitable for exchange matrix calculations on massively parallel computing architectures, such as graphical processing units, for which we present a first linear-scaling exchange matrix evaluation algorithm.

Efficient distance-including integral screening in linear-scaling Møller-Plesset perturbation theory.

Efficient estimates for the preselection of two-electron integrals in atomic-orbital based Møller-Plesset perturbation theory (AO-MP2) theory are presented, which allow for evaluating the AO-MP2 energy with computational effort that scales linear with molecular size for systems with a significant HOMO-LUMO gap. The estimates are based on our recently introduced QQR approach [S. A. Maurer, D. S. Lambrecht, D. Flaig, and C. Ochsenfeld, J. Chem. Phys. 136, 144107 (2012)], which exploits the asympotic decay of the integral values with increasing bra-ket separation as deduced from the multipole expansion and combines this decay behavior with the common Schwarz bound to a tight and simple estimate. We demonstrate on a diverse selection of benchmark systems that our AO-MP2 method in combination with the QQR-type estimates produces reliable results for systems with both localized and delocalized electronic structure, while in the latter case the screening essentially reverts to the common Schwarz screening. For systems with localized electronic structure, our AO-MP2 method shows an early onset of linear scaling as demonstrated on DNA systems. The favorable scaling behavior allows to compute systems with more than 1000 atoms and 10,000 basis functions on a single core that are clearly not accessible with conventional MP2 methods. Furthermore, our AO-MP2 method is particularly suited for parallelization and we present benchmark calculations on a protein-DNA repair complex comprising 2025 atoms and 20,371 basis functions.

Factors That Determine the Variation of Equilibrium and Kinetic Properties of QM/MM Enzyme Simulations: QM Region, Conformation, and Boundary Condition.

To analyze the impact of various technical details on the results of quantum mechanical (QM)/molecular mechanical (MM) enzyme simulations, including the QM region size, catechol-O-methyltransferase (COMT) is studied as a model system using an approximate QM/MM method (DFTB3/CHARMM). The results show that key equilibrium and kinetic properties for methyl transfer in COMT exhibit limited variations with respect to the size of the QM region, which ranges from ∼100 to ∼500 atoms in this study. With extensive sampling, local and global structural characteristics of the enzyme are largely conserved across the studied QM regions, while the nature of the transition state (e.g., secondary kinetic isotope effect) and reaction exergonicity are largely maintained. Deviations in the free energy profile with different QM region sizes are similar in magnitude to those observed with changes in other simulation protocols, such as different initial enzyme conformations and boundary conditions. Electronic structural properties, such as the covariance matrix of residual charge fluctuations, appear to exhibit rather long-range correlations, especially when the peptide backbone is included in the QM region; this observation holds when a range-separated DFT approach is used as the QM region, suggesting that delocalization error is unlikely the origin. Overall, the analyses suggest that multiple simulation details determine the results of QM/MM enzyme simulations with comparable contributions.

Accelerating Hybrid Density Functional Theory Molecular Dynamics Simulations by Seminumerical Integration, Resolution-of-the-Identity Approximation, and Graphics Processing Units.

The computationally very demanding evaluation of the 4-center-2-electron (4c2e) integrals and their respective integral derivatives typically represents the major bottleneck within hybrid Kohn-Sham density functional theory molecular dynamics simulations. Building upon our previous works on seminumerical exact-exchange (sn-LinK) [Laqua, H., Thompsons, T. H., Kussmann, J., Ochsenfeld, C., J. Chem. Theory Comput. 2020, 16, 1465] and resolution-of-the-identity Coulomb (RI-J) [Kussmann, J., Laqua, H., Ochsenfeld, C., J. Chem. Theory Comput. 2021, 17, 1512], the expensive 4c2e integral evaluation can be avoided entirely, resulting in a highly efficient electronic structure theory method, allowing for fast ab initio molecular dynamics (AIMD) simulations even with large basis sets. Moreover, we propose to combine the final self-consistent field (SCF) step with the subsequent nuclear forces evaluation, providing the forces at virtually no additional cost after a converged SCF calculation, reducing the total runtime of an AIMD simulation by about another 25%. In addition, multiple independent MD trajectories can be computed concurrently on a single node, leading to a greatly increased utilization of the available hardware─especially when combined with graphics processing unit acceleration─improving the overall throughput by up to another 5 times in this way. With all of those optimizations combined, our proposed method provides nearly 3 orders of magnitude faster execution times than traditional 4c2e integral-based methods. To demonstrate the practical utility of the approach, quantum-mechanical/molecular-mechanical dynamics simulations on double-stranded DNA were performed, investigating the relative hydrogen bond strength between adenine-thymine and guanine-cytosine base pairs. In addition, this illustrative application also contains a general accuracy assessment of the introduced approximations (integration grids, resolution-of-the-identity) within AIMD simulations, serving as a protocol on how to apply these new methods to practical problems.

Nuclei-selected NMR shielding calculations: a sublinear-scaling quantum-chemical method.

An ab initio method for the direct calculation of NMR shieldings for selected nuclei at the Hartree-Fock and density-functional theory level is presented. Our method shows a computational effort scaling only sublinearly with molecular size, as it is motivated by the physical consideration that the chemical shielding is dominated by its local environment. The key feature of our method is to avoid the conventionally performed calculation of all NMR shieldings but instead to solve directly for specific nuclear shieldings. This has important implications not only for the study of large molecules, but also for the simulation of solvent effects and molecular dynamics, since often just a few shieldings are of interest. Our theory relies on two major aspects both necessary to provide a sublinear scaling behavior: First, an alternative expression for the shielding tensor is derived, which involves the response density matrix with respect to the nuclear magnetic moment instead of the response to the external magnetic field. Second, as unphysical long-range contributions occur within the description of distributed gauge origin methods that do not influence the final expectation value, we present a screening procedure to truncate the B-field dependent basis set, which is crucial in order to ensure an early onset of the sublinear scaling. The screening is in line with the r(-2) distance decay of Biot-Savarts law for induced magnetic fields. Our present truncation relies on the introduced concept of "individual gauge shielding contributions" applied to a reformulated shielding tensor, the latter consisting of gauge-invariant terms. The presented method is generally applicable and shows typical speed-ups of about one order of magnitude; moreover, due to the reduced scaling behavior of O(1) as compared to O(N), the wins become larger with increasing system size. We illustrate the validity of our method for several test systems, including ring-current dominated systems and biomolecules with more than 1000 atoms.

Highly Efficient Resolution-of-Identity Density Functional Theory Calculations on Central and Graphics Processing Units.

We present an efficient method to evaluate Coulomb potential matrices using the resolution of identity approximation and semilocal exchange-correlation potentials on central (CPU) and graphics processing units (GPU). The new GPU-based RI-algorithm shows a high performance and ensures the favorable scaling with increasing basis set size as the conventional CPU-based method. Furthermore, our method is based on the J-engine algorithm [White; , Head-Gordon, J. Chem. Phys. 1996, 7, 2620], which allows for further optimizations that also provide a significant improvement of the corresponding CPU-based algorithm. Due to the increased performance for the Coulomb evaluation, the calculation of the exchange-correlation potential of density functional theory on CPUs quickly becomes a bottleneck to the overall computational time. Hence, we also present a GPU-based algorithm to evaluate the exchange-correlation terms, which results in an overall high-performance method for density functional calculations. The algorithms to evaluate the potential and nuclear derivative terms are discussed, and their performance on CPUs and GPUs is demonstrated for illustrative calculations.

Monte Carlo simulations of single- and multistep enzyme-catalyzed reaction sequences: effects of diffusion, cell size, enzyme fluctuations, colocalization, and segregation.

Following the discovery of slow fluctuations in the catalytic activity of an enzyme in single-molecule experiments, it has been shown that the classical Michaelis-Menten (MM) equation relating the average enzymatic velocity and the substrate concentration may hold even for slowly fluctuating enzymes. In many cases, the average velocity is that given by the MM equation with time-averaged values of the fluctuating rate constants and the effect of enzyme fluctuations is simply averaged out. The situation is quite different for a sequence of reactions. For colocalization of a pair of enzymes in a sequence to be effective in promoting reaction, the second must be active when the first is active or soon after. If the enzymes are slowly varying and only rarely active, the product of the first reaction may diffuse away before the second enzyme is active, and colocalization may have little value. Even for single-step reactions the interplay of reaction and diffusion with enzyme fluctuations leads to added complexities, but for multistep reactions the interplay of reaction and diffusion, cell size, compartmentalization, enzyme fluctuations, colocalization, and segregation is far more complex than for single-step reactions. In this paper, we report the use of stochastic simulations at the level of whole cells to explore, understand, and predict the behavior of single- and multistep enzyme-catalyzed reaction systems exhibiting some of these complexities. Results for single-step reactions confirm several earlier observations by others. The MM relationship, with altered constants, is found to hold for single-step reactions slowed by diffusion. For single-step reactions, the distribution of enzymes in a regular grid is slightly more effective than a random distribution. Fluctuations of enzyme activity, with average activity fixed, have no observed effects for simple single-step reactions slowed by diffusion. Two-step sequential reactions are seen to be slowed by segregation of the enzymes for each step, and results of the calculations suggest limits for cell size. Colocalization of enzymes for a two-step sequence is seen to promote reaction, and rates fall rapidly with increasing distance between enzymes. Low frequency fluctuations of the activities of colocalized enzymes, with average activities fixed, can greatly reduce reaction rates for sequential reactions.

Combining Graphics Processing Units, Simplified Time-Dependent Density Functional Theory, and Finite-Difference Couplings to Accelerate Nonadiabatic Molecular Dynamics.

Starting from our recently published implementation of nonadiabatic molecular dynamics (NAMD) on graphics processing units (GPUs), we explore further approaches to accelerate ab initio NAMD calculations at the time-dependent density functional theory (TDDFT) level of theory. We employ (1) the simplified TDDFT schemes of Grimme et al. and (2) the Hammes-Schiffer-Tully approach to obtain nonadiabatic couplings from finite-difference calculations. The resulting scheme delivers an accurate physical picture while virtually eliminating the two computationally most demanding steps of the algorithm. Combined with our GPU-based integral routines for SCF, TDDFT, and TDDFT derivative calculations, NAMD simulations of systems of a few hundreds of atoms at a reasonable time scale become accessible on a single compute node. To demonstrate this and to present a first, illustrative example, we perform TDDFT/MM-NAMD simulations of the rhodopsin protein.

Highly Efficient, Linear-Scaling Seminumerical Exact-Exchange Method for Graphic Processing Units.

We present a highly efficient and asymptotically linear-scaling graphic processing unit accelerated seminumerical exact-exchange method (sn-LinK). We go beyond our previous central processing unit-based method (Laqua, H.; Kussmann, J.; Ochsenfeld, C. J. Chem. Theory Comput. 2018, 14, 3451-3458) by employing our recently developed integral bounds (Thompson, T. H.; Ochsenfeld, C. J. Chem. Phys. 2019, 150, 044101) and high-accuracy numerical integration grid (Laqua, H.; Kussmann, J.; Ochsenfeld, C. J. Chem. Phys. 2018, 149, 204111). The accuracy is assessed for several established test sets, providing errors significantly below 1mEh for the smallest grid. Moreover, a comprehensive performance analysis for large molecules between 62 and 1347 atoms is provided, revealing the outstanding performance of our method, in particular, for large basis sets such as the polarized quadruple-zeta level with diffuse functions.