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.

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.

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.

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.

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.

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.

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.

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.

Nonadiabatic Molecular Dynamics on Graphics Processing Units: Performance and Application to Rotary Molecular Motors.

Nonadiabatic molecular dynamics (NAMD) simulations of molecular systems require the efficient evaluation of excited-state properties, such as energies, gradients, and nonadiabatic coupling vectors. Here, we investigate the use of graphics processing units (GPUs) in addition to central processing units (CPUs) to efficiently calculate these properties at the time-dependent density functional theory (TDDFT) level of theory. Our implementation in the FermiONs++ program package uses the J-engine and a preselective screening procedure for the calculation of Coulomb and exchange kernels, respectively. We observe good speed-ups for small and large molecular systems (comparable to those observed in ground-state calculations) and reduced (down to sublinear) scaling behavior with respect to the system size (depending on the spatial locality of the investigated excitation). As a first illustrative application, we present efficient NAMD simulations of a series of newly designed light-driven rotary molecular motors and compare their S1 lifetimes. Although all four rotors show different S1 excitation energies, their ability to rotate upon excitation is conserved, making the series an interesting starting point for rotary molecular motors with tunable excitation energies.

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.

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.

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.

Efficient and Linear-Scaling Seminumerical Method for Local Hybrid Density Functionals.

Local hybrid functionals, that is, functionals with local dependence on the exact exchange energy density, generalize the popular class of global hybrid functionals and extend the applicability of density functional theory to electronic structures that require an accurate description of static correlation. However, the higher computational cost compared to conventional Kohn-Sham density functional theory restrained their widespread application. Here, we present a low-prefactor, linear-scaling method to evaluate the local hybrid exchange-correlation potential as well as the corresponding nuclear forces by employing a seminumerical integration scheme. In the seminumerical scheme, one integration in the electron repulsion integrals is carried out analytically and the other one is carried out numerically, employing an integration grid. A high computational efficiency is achieved by combining the preLinK method [ J. Kussmann and C. Ochsenfeld, J. Chem. Phys. 2013 138, 134114 ] with explicit screening of integrals for batches of grid points to minimize the screening overhead. This new method, denoted as preLinX, provides an 8-fold performance increase for a DNA fragment containing four base pairs as compared to existing S- and P-junction-based methods. In this way, our method allows for the evaluation of local hybrid functionals at a cost similar to that of global hybrid functionals. The linear-scaling behavior, efficiency, accuracy, and multinode parallelization of the presented method is demonstrated for large systems containing more than 1000 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.

Efficient and Accurate Born-Oppenheimer Molecular Dynamics for Large Molecular Systems.

An efficient scheme for the calculation of Born-Oppenheimer molecular dynamics (BOMD) simulations is introduced. It combines the corrected small basis set Hartree-Fock (HF-3c) method by Sure and Grimme [J. Comput. Chem. 2013, 43, 1672], extended Lagrangian BOMD (XL-BOMD) by Niklasson et al. [J. Chem. Phys. 2009, 130, 214109], and the calculation of the two electron integrals on graphics processing units (GPUs) [J. Chem. Phys. 2013, 138, 134114; J. Chem. Theory Comput. 2015, 11, 918]. To explore the parallel performance of our strong scaling implementation of the method, we present timings and extract, as its validation and first illustrative application, high-quality vibrational spectra from simulated trajectories of ß-carotene, paclitaxel, and liquid water (up to 500 atoms). We conclude that the presented BOMD scheme may be used as a cost-efficient and reliable tool for computing vibrational spectra and thermodynamics of large molecular systems including explicit solvent molecules containing 500 atoms and more. Simulating 50 ps of maitotoxin (nearly 500 atoms) employing time steps of 0.5 fs requires ∼3 weeks on 12 CPUs (Intel Xeon E5 2620 v3) with 24 GPUs (AMD FirePro 3D W8100).

Employing OpenCL to Accelerate Ab Initio Calculations on Graphics Processing Units.

We present an extension of our graphics processing units (GPU)-accelerated quantum chemistry package to employ OpenCL compute kernels, which can be executed on a wide range of computing devices like CPUs, Intel Xeon Phi, and AMD GPUs. Here, we focus on the use of AMD GPUs and discuss differences as compared to CUDA-based calculations on NVIDIA GPUs. First illustrative timings are presented for hybrid density functional theory calculations using serial as well as parallel compute environments. The results show that AMD GPUs are as fast or faster than comparable NVIDIA GPUs and provide a viable alternative for quantum chemical applications.

Hybrid CPU/GPU Integral Engine for Strong-Scaling Ab Initio Methods.

We present a parallel integral algorithm for two-electron contributions occurring in Hartree-Fock and hybrid density functional theory that allows for a strong scaling parallelization on inhomogeneous compute clusters. With a particular focus on graphic processing units, we show that our approach allows an efficient use of CPUs and graphics processing units (GPUs) simultaneously, although the different architectures demand conflictive strategies in order to ensure efficient program execution. Furthermore, we present a general strategy to use large basis sets like quadruple-ζ split valence on GPUs and investigate the balance between CPUs and GPUs depending on l-quantum numbers of the corresponding basis functions. Finally, we present first illustrative calculations using a hybrid CPU/GPU environment and demonstrate the strong-scaling performance of our parallelization strategy also for pure CPU-based calculations.

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.

Nuclear Magnetic Shieldings of Stacked Aromatic and Antiaromatic Molecules.

Nuclear magnetic shieldings have been calculated at the density functional theory (DFT) level for stacks of benzene, hexadehydro[12]annulene, dodecadehydro[18]annulene, and hexabenzocoronene. The magnetic shieldings due to the ring currents in the adjacent molecules have been estimated by calculating nucleus independent molecular shieldings for the monomer in the atomic positions of neighbor molecules. The calculations show that the independent shielding model works reasonably well for the 1H NMR shieldings of benzene and hexadehydro[12]annulene, whereas for the larger molecules and for the 13C NMR shieldings the interaction between the molecules leads to shielding effects that are at least of the same size as the ring current contributions from the adjacent molecules. A better agreement is obtained when the nearest neighbors are also considered at full quantum mechanical (QM) level. The calculations suggest that the nearest solvent molecules must be included in the quantum mechanical system, at least when estimating solvent shifts at the molecular mechanics (MM) level. Current density calculations show that the stacking does not significantly affect the ring current strengths of the individual molecules, whereas the shape of the ring current for a single molecule differs from that of the stacked molecules.