Achieving an Order of Magnitude Speedup in Hybrid-Functional- and Plane-Wave-Based Ab Initio Molecular Dynamics: Applications to Proton-Transfer Reactions in Enzymes and in Solution

Complex-Valued K-Means Clustering of Interpolative Separable Density Fitting Algorithm for Large-Scale Hybrid Functional Enabled Ab Initio Molecular Dynamics Simulations within Plane Waves

K-means clustering, as a classic unsupervised machine learning algorithm, is the key step to select the interpolation sampling points in interpolative separable density fitting (ISDF) decomposition for hybrid functional electronic structure calculations. Real-valued K-means clustering for accelerating the ISDF decomposition has been demonstrated for large-scale hybrid functional enabled ab initio molecular dynamics (hybrid AIMD) simulations within plane-wave basis sets where the Kohn-Sham orbitals are real-valued. However, it is unclear whether such K-means clustering works for complex-valued Kohn-Sham orbitals. Here, we propose an improved weight function defined as the sum of the square modulus of complex-valued Kohn-Sham orbitals in K-means clustering for hybrid AIMD simulations. Numerical results demonstrate that the K-means algorithm with a new weight function yields smoother and more delocalized interpolation sampling points, resulting in smoother energy potential, smaller energy drift, and longer time steps for hybrid AIMD simulations compared to the previous weight function used in the real-valued K-means algorithm. In particular, we find that this improved algorithm can obtain more accurate oxygen-oxygen radial distribution functions in liquid water molecules and a more accurate power spectrum in crystal silicon dioxide compared to the previous K-means algorithm. Finally, we describe a massively parallel implementation of this ISDF decomposition to accelerate large-scale complex-valued hybrid AIMD simulations containing thousands of atoms (2,744 atoms), which can scale up to 5,504 CPU cores on modern supercomputers.

Can simple 'molecular' corrections outperform projector augmented-wave density functional theory in the prediction of 35 Cl electric field gradient tensor parameters for chlorine-containing crystalline systems?

Many-body expansion (MBE) fragment approaches have been applied to accurately compute nuclear magnetic resonance (NMR) parameters in crystalline systems. Recent examples demonstrate that electric field gradient (EFG) tensor parameters can be accurately calculated for 14 N and 17 O. A key additional development is the simple molecular correction (SMC) approach, which uses two one-body fragment (i.e., isolated molecule) calculations to adjust NMR parameter values established using 'benchmark' projector augmented-wave (PAW) density functional theory (DFT) values. Here, we apply a SMC using the hybrid PBE0 exchange-correlation (XC) functional to see if this can improve the accuracy of calculated 35 Cl EFG tensor parameters. We selected eight organic and two inorganic crystal structures and considered 15 chlorine sites. We find that this SMC improves the accuracy of computed values for both the 35 Cl quadrupolar coupling constant (CQ ) and the asymmetry parameter ( η Q ) by approximately 30% compared with benchmark PAW DFT values. We also assessed a SMC that offers local improvements not only in terms of the quality of the XC functional but simultaneously in the quality of the description of relativistic effects via the inclusion of spin-orbit effects. As the inorganic systems considered contain heavy atoms bonded to the chlorine atoms, we find further improvements in the accuracy of calculated 35 Cl EFG tensor parameters when both a hybrid functional and spin-orbit effects are included in the SMC. On the contrary, for chlorine-containing organics, the inclusion of spin-orbit relativistic effects using a SMC does not improve the accuracy of computed 35 Cl EFG tensor parameters.

Speeding-up Hybrid Functional-Based Ab Initio Molecular Dynamics Using Multiple Time-stepping and Resonance-Free Thermostat

Ab initio molecular dynamics (AIMD) based on density functional theory (DFT) has become a workhorse for studying the structure, dynamics, and reactions in condensed matter systems. Currently, AIMD simulations are primarily carried out at the level of generalized gradient approximation (GGA), which is at the second rung of DFT functionals in terms of accuracy. Hybrid DFT functionals, which form the fourth rung in the accuracy ladder, are not commonly used in AIMD simulations as the computational cost involved is 100 times or higher. To facilitate AIMD simulations with hybrid functionals, we propose here an approach using multiple time stepping with adaptively compressed exchange operator and resonance-free thermostat, that could speed up the calculations by ∼30 times or more for systems with a few hundred of atoms. We demonstrate that by achieving this significant speed up and making the compute time of hybrid functional-based AIMD simulations at par with that of GGA functionals, we are able to study several complex condensed matter systems and model chemical reactions in solution with hybrid functionals that were earlier unthinkable to be performed.

High-Throughput Condensed-Phase Hybrid Density Functional Theory for Large-Scale Finite-Gap Systems: The SeA Approach

High-throughput electronic structure calculations (often performed using density functional theory (DFT)) play a central role in screening existing and novel materials, sampling potential energy surfaces, and generating data for machine learning applications. By including a fraction of exact exchange (EXX), hybrid functionals reduce the self-interaction error in semilocal DFT and furnish a more accurate description of the underlying electronic structure, albeit at a computational cost that often prohibits such high-throughput applications. To address this challenge, we have constructed a robust, accurate, and computationally efficient framework for high-throughput condensed-phase hybrid DFT and implemented this approach in the PWSCF module of Quantum ESPRESSO (QE). The resulting SeA approach (SeA = SCDM + exx + ACE) combines and seamlessly integrates: (i) the selected columns of the density matrix method (SCDM, a robust noniterative orbital localization scheme that sidesteps system-dependent optimization protocols), (ii) a recently extended version of exx (a black-box linear-scaling EXX algorithm that exploits sparsity between localized orbitals in real space when evaluating the action of the standard/full-rank V^xx operator), and (iii) adaptively compressed exchange (ACE, a low-rank V^xx approximation). In doing so, SeA harnesses three levels of computational savings: pair selection and domain truncation from SCDM + exx (which only considers spatially overlapping orbitals on orbital-pair-specific and system-size-independent domains) and low-rank V^xx approximation from ACE (which reduces the number of calls to SCDM + exx during the self-consistent field (SCF) procedure). Across a diverse set of 200 nonequilibrium (H₂O)₆₄ configurations (with densities spanning 0.4-1.7 g/cm³), SeA provides a 1-2 order-of-magnitude speedup in the overall time-to-solution, i.e., ≈8-26× compared to the convolution-based PWSCF(ACE) implementation in QE and ≈78-247× compared to the conventional PWSCF(Full) approach, and yields energies, ionic forces, and other properties with high fidelity. As a proof-of-principle high-throughput application, we trained a deep neural network (DNN) potential for ambient liquid water at the hybrid DFT level using SeA via an actively learned data set with ≈8,700 (H₂O)₆₄ configurations. Using an out-of-sample set of (H₂O)₅₁₂ configurations (at nonambient conditions), we confirmed the accuracy of this SeA-trained potential and showcased the capabilities of SeA by computing the ground-truth ionic forces in this challenging system containing >1,500 atoms.

Hybrid Functional and Plane Waves based Ab Initio Molecular Dynamics Study of the Aqueous Fe2+ /Fe3+ Redox Reaction

Kohn-Sham density functional theory and plane wave basis set based ab initio molecular dynamics (AIMD) simulation is a powerful tool for studying complex reactions in solutions, such as electron transfer (ET) reactions involving Fe²⁺ /Fe³⁺ ions in water. In most cases, such simulations are performed using density functionals at the level of Generalized Gradient Approximation (GGA). The challenge in modelling ET reactions is the poor quality of GGA functionals in predicting properties of such open-shell systems due to the inevitable self-interaction error (SIE). While hybrid functionals can minimize SIE, standard plane-wave based AIMD at that level of theory is typically 150 times slower than GGA for systems containing ∼100 atoms. Among several approaches reported to speed-up AIMD simulations with hybrid functionals, the noise-stabilized MD (NSMD) procedure, together with the use of localized orbitals to compute the required exchange integrals, is an attractive option. In this work, we demonstrate the application of the NSMD approach for studying the Fe²⁺ /Fe³⁺ redox reaction in water. It is shown here that long AIMD trajectories at the level of hybrid density functionals can be obtained using this approach. Redox properties of the aqueous Fe²⁺ /Fe³⁺ system computed from these simulations are compared with the available experimental data for validation.

Improving the scaling and performance of multiple time stepping-based molecular dynamics with hybrid density functionals

Density functionals at the level of the generalized gradient approximation (GGA) and a plane-wave basis set are widely used today to perform ab initio molecular dynamics (AIMD) simulations. Going up in the ladder of accuracy of density functionals from GGA (second rung) to hybrid density functionals (fourth rung) is much desired pertaining to the accuracy of the latter in describing structure, dynamics, and energetics of molecular and condensed matter systems. On the other hand, hybrid density functional based AIMD simulations are about two orders of magnitude slower than GGA based AIMD for systems containing ~100 atoms using ~100 compute cores. Two methods, namely MTACE and s-MTACE, based on a multiple time step integrator and adaptively compressed exchange operator formalism are able to provide a speed-up of about 7-9 in performing hybrid density functional based AIMD. In this work, we report an implementation of these methods using a task-group based parallelization within the CPMD program package, with the intention to take advantage of the large number of compute cores available on modern high-performance computing platforms. We present here the boost in performance achieved through this algorithm. This work also identifies the computational bottleneck in the s-MTACE method and proposes a way to overcome it.

Free Energy Profiles of Proton Transfer Reactions: Density Functional Benchmark from Biased Ab Initio Dynamics

By coupling an enhanced sampling algorithm with an orbital-localized variant of Car-Parrinello molecular dynamics, the so-called atomic centered density matrix propagation model, we reconstruct the free energy profiles along reaction pathways using different density functional approximations (DFAs) ranging from locals to hybrids. In particular, we compare the computed free energy barrier height of proton transfer (PT) reactions to those obtained by a more traditional static approach, based on the intrinsic reaction coordinate (IRC), for two case systems, namely malonaldehyde and formic acid dimer. The obtained results show that both the IRC profiles and the potentials of mean force, derived from biased dynamic trajectories, are very sensitive to the density functional approximation applied. More precisely, we observe that, with the notable exception of M06-L, local density functionals always strongly underestimate the reaction barrier heights. More generally, we find that also the shape of the free energy profile is very sensitive to the density functional choice, thus highlighting the effect, often neglected, that the choice of DFA has also in the case of dynamics simulations.

Boosting the conformational sampling by combining replica exchange with solute tempering and well-sliced metadynamics

Methods that combine collective variable (CV) based enhanced sampling and global tempering approaches are used in speeding-up the conformational sampling and free energy calculation of large and soft systems with a plethora of energy minima. In this paper, a new method of this kind is proposed in which the well-sliced metadynamics approach (WSMTD) is united with replica exchange with solute tempering (REST2) method. WSMTD employs a divide-and-conquer strategy wherein high-dimensional slices of a free energy surface are independently sampled and combined. The method enables one to accomplish a controlled exploration of the CV-space with a restraining bias as in umbrella sampling, and enhance-sampling of one or more orthogonal CVs using a metadynamics like bias. The new hybrid method proposed here enables boosting the sampling of more slow degrees of freedom in WSMTD simulations, without the need to specify associated CVs, through a replica exchange scheme within the framework of REST2. The high-dimensional slices of the probability distributions of CVs computed from the united WSMTD and REST2 simulations are subsequently combined using the weighted histogram analysis method to obtain the free energy surface. We show that the new method proposed here is accurate, improves the conformational sampling, and achieves quick convergence in free energy estimates. We demonstrate this by computing the conformational free energy landscapes of solvated alanine tripeptide and Trp-cage mini protein in explicit water.

Enabling Large-Scale Condensed-Phase Hybrid Density Functional Theory-Based Ab Initio Molecular Dynamics II: Extensions to the Isobaric-Isoenthalpic and Isobaric-Isothermal Ensembles

In the previous paper of this series [Ko, H.-Y. et al. J. Chem. Theory Comput. 2020, 16, 3757-3785], we presented a theoretical and algorithmic framework based on a localized representation of the occupied space that exploits the inherent sparsity in the real-space evaluation of the exact exchange (EXX) interaction in finite-gap systems. This was accompanied by a detailed description of exx, a massively parallel hybrid message-passing interface MPI/OpenMP implementation of this approach in Quantum ESPRESSO (QE) that enables linear scaling hybrid density functional theory (DFT)-based ab initio molecular dynamics (AIMD) in the microcanonical/canonical (NVE/NVT) ensembles of condensed-phase systems containing 500-1000 atoms (in fixed orthorhombic cells) with a wall time cost comparable to semi-local DFT. In this work, we extend the current capabilities of exx to enable hybrid DFT-based AIMD simulations of large-scale condensed-phase systems with general and fluctuating cells in the isobaric-isoenthalpic/isobaric-isothermal (NpH/NpT) ensembles. The theoretical extensions to this approach include an analytical derivation of the EXX contribution to the stress tensor for systems in general simulation cells with a computational complexity that scales linearly with system size. The corresponding algorithmic extensions to exx include optimized routines that (i) handle both static and fluctuating simulation cells with non-orthogonal lattice symmetries, (ii) solve Poisson's equation in general/non-orthogonal cells via an automated selection of the auxiliary grid directions in the Natan-Kronik representation of the discrete Laplacian operator, and (iii) evaluate the EXX contribution to the stress tensor. Using this approach, we perform a case study on a variety of condensed-phase systems (including liquid water, a benzene molecular crystal polymorph, and semi-conducting crystalline silicon) and demonstrate that the EXX contributions to the energy and stress tensor simultaneously converge with an appropriate choice of exx parameters. This is followed by a critical assessment of the computational performance of the extended exx module across several different high-performance computing architectures via case studies on (i) the computational complexity due to lattice symmetry during NpT simulations of three different ice polymorphs (i.e., ice Ih, II, and III) and (ii) the strong/weak parallel scaling during large-scale NpT simulations of liquid water. We demonstrate that the robust and highly scalable implementation of this approach in the extended exx module is capable of evaluating the EXX contribution to the stress tensor with negligible cost (<1%) as well as all other EXX-related quantities needed during NpT simulations of liquid water (with a very tight 150 Ry planewave cutoff) in ≈5.2 s ((H₂O)₁₂₈) and ≈6.8 s ((H₂O)₂₅₆) per AIMD step. As such, the extended exx module presented in this work brings us another step closer to routinely performing hybrid DFT-based AIMD simulations of sufficient duration for large-scale condensed-phase systems across a wide range of thermodynamic conditions.