Multiple time step integrators in ab initio molecular dynamics

Assessing the persistence of chalcogen bonds in solution with neural network potentials

Non-covalent bonding patterns are commonly harvested as a design principle in the field of catalysis, supramolecular chemistry and functional materials to name a few. Yet, their computational description generally neglects finite temperature and environment effects, which promote competing interactions and alter their static gas-phase properties. Recently, neural network potentials (NNPs) trained on Density Functional Theory (DFT) data have become increasingly popular to simulate molecular phenomena in condensed phase with an accuracy comparable to ab initio methods. To date, most applications have centered on solid-state materials or fairly simple molecules made of a limited number of elements. Herein, we focus on the persistence and strength of chalcogen bonds involving a benzotelluradiazole in condensed phase. While the tellurium-containing heteroaromatic molecules are known to exhibit pronounced interactions with anions and lone pairs of different atoms, the relevance of competing intermolecular interactions, notably with the solvent, is complicated to monitor experimentally but also challenging to model at an accurate electronic structure level. Here, we train direct and baselined NNPs to reproduce hybrid DFT energies and forces in order to identify what are the most prevalent non-covalent interactions occurring in a solute-Cl$^-$-THF mixture. The simulations in explicit solvent highlight competition with chalcogen bonds formed with the solvent and the short-range directionality of the interaction with direct consequences for the molecular properties in the solution. The comparison with other potentials (e.g., AMOEBA, direct NNP and continuum solvent model) also demonstrates that baselined NNPs offer a reliable picture of the non-covalent interaction interplay occurring in solution.

Recent Advances in First-Principles Based Molecular Dynamics

First-principles molecular dynamics (FPMD) and its quantum mechanical-molecular mechanical (QM/MM) extensions are powerful tools to follow the real-time dynamics of a broad variety of systems in their ground as well as electronically excited states. The continued advances in computational power have enabled simulations of QM regions of larger sizes for more extended time scales. In addition, development of the parallel algorithms has boosted the performance of QM/MM methods even on existing computer architectures. In the case of density functional-based FPMD, systems of several hundreds to thousands of atoms can now be customarily simulated for tens to hundreds of picoseconds. In spite of this progress, the time scale limitations remain severe, especially when high-rung exchange-correlation functionals or high-level wave function based quantum mechanical methods are used. To ameliorate this, a large number of enhanced sampling methods have been introduced but most of the approaches that have been developed to increase the efficiency of FPMD based simulations sacrifice the real-time dynamics in favor of enhancing sampling. Here, we present some recent advances in boosting the efficiency of FPMD based simulations while keeping the full dynamic information. These include a highly efficient recent implementation of FPMD-based QM/MM simulations that not only enables fully flexible combinations of different electronic structure methods and force fields via a highly efficient communication library, it also fully exploits parallelism for both quantum and classical descriptions. The second type of acceleration methods we discuss is a large family of specially devised multiple-time-step algorithms that make use of suitable breakups of the total nuclear forces into fast components that can be calculated via lower level methods and slowly varying correction forces evaluated with a high-level method at long time intervals. The computational gain of this scheme mostly depends on the cost difference between the two methods and advantageous combinations can yield large speedups without compromising the accuracy of the high-level method. And finally, the third class of FPMD acceleration methods presented here are machine learning models to accelerated FPMD and their powerful combinations with multiple-time-step techniques. The combination of all the approaches enables substantial speedups of FPMD simulations of several orders of magnitude while fully preserving the real-time dynamics and accuracy.

Evolution of microscopic heterogeneity and dynamics in choline chloride-based deep eutectic solvents

Deep eutectic solvents (DESs) are an emerging class of non-aqueous solvents that are potentially scalable, easy to prepare and functionalize for many applications ranging from biomass processing to energy storage technologies. Predictive understanding of the fundamental correlations between local structure and macroscopic properties is needed to exploit the large design space and tunability of DESs for specific applications. Here, we employ a range of computational and experimental techniques that span length-scales from molecular to macroscopic and timescales from picoseconds to seconds to study the evolution of structure and dynamics in model DESs, namely Glyceline and Ethaline, starting from the parent compounds. We show that systematic addition of choline chloride leads to microscopic heterogeneities that alter the primary structural relaxation in glycerol and ethylene glycol and result in new dynamic modes that are strongly correlated to the macroscopic properties of the DES formed.

Tailoring the macroscopic properties of deep eutectic solvents requires knowing how these depend on the local structure and microscopic dynamics. The authors, with computational and experimental tools spanning a wide range of space- and timescales, shed light into the relationship between micro and macroscopic properties in glyceline and ethaline.

Characterizing and Contrasting Structural Proton Transport Mechanisms in Azole Hydrogen Bond Networks Using Ab Initio Molecular Dynamics

Imidazole and 1,2,3-triazole are promising hydrogen-bonded heterocycles that conduct protons via a structural mechanism and whose derivatives are present in systems ranging from biological proton channels to proton exchange membrane fuel cells. Here, we leverage multiple time-stepping to perform ab initio molecular dynamics of imidazole and 1,2,3-triazole at the nanosecond time scale. We show that despite the close structural similarities of these compounds, their proton diffusion constants vary by over an order of magnitude. Our simulations reveal the reasons for these differences in diffusion constants, which range from the degree of hydrogen-bonded chain linearity to the effect of the central nitrogen atom in 1,2,3-triazole on proton transport. In particular, we uncover evidence of two "blocking" mechanisms in 1,2,3-triazole, where covalent and hydrogen bonds formed by the central nitrogen atom limit the mobility of protons. Our simulations thus provide insights into the origins of the experimentally observed 10-fold difference in proton conductivity.

Re-integration with anchor points algorithm for ab initio molecular dynamics

A new integration scheme for ab initio molecular dynamics (MD) is proposed in this work for efficient propagation using large time steps (e.g., 2.0 fs or a larger time step with one ab initio evaluation of gradients for the dynamics point and one additional evaluation for the anchor point per dynamics step). This algorithm is called re-integration with anchor points (REAP) integrator. The REAP integrator starts from a quadratic potential energy surface based on the updated Hessian to propagate the system to the halfway of the MD step that is called the anchor point. Then, an approximate dynamics position for this step is obtained by the propagation based on an interpolated surface using the anchor point and the previous MD point. The approximate dynamics step can be further improved by the re-integration steps, i.e., integration based on the interpolated surface using the calculated energies, gradients, and updated Hessians of the previous step, the anchor point, and the approximate current step. A trajectory only needs one analytical Hessian calculation at the initial geometry, and thereafter, only calculations of gradients are required. This integrator can be considered either as a generalization of Hessian-based predictor-corrector integration with substantial improvement of accuracy and efficiency or as a dynamics on interpolated surfaces that are built on the fly. An automatic correction scheme is implemented by comparing the interpolated energies and gradients to the actual ones to ensure the quality of the interpolations at a certain level. The tests in this work show that the REAP method can increase computational efficiency by more than one order of magnitude than that of the velocity Verlet integrator and more than twice that of Hessian-based predictor-corrector integration.

Nuclear Quantum Effect and Its Temperature Dependence in Liquid Water from Random Phase Approximation via Artificial Neural Network

We report structural and dynamical properties of liquid water described by the random phase approximation (RPA) correlation together with the exact exchange energy (EXX) within density functional theory. By utilizing thermostated ring polymer molecular dynamics, we examine the nuclear quantum effects and their temperature dependence. We circumvent the computational limitation of performing direct first-principles molecular dynamics simulation at this high level of electronic structure theory by adapting an artificial neural network model. We show that the EXX+RPA level of theory accurately describes liquid water in terms of both dynamical and structural properties.

Elucidating the Proton Transport Pathways in Liquid Imidazole with First-Principles Molecular Dynamics

Imidazole is a promising anhydrous proton conductor with a high conductivity comparable to that of water at a similar temperature relative to its melting point. Previous theoretical studies of the mechanism of proton transport in imidazole have relied either on empirical models or on ab initio trajectories that have been too short to draw significant conclusions. Here, we present the results of multiple time-step ab initio molecular dynamics simulations of an excess proton in liquid imidazole reaching 1 ns in total simulation time. We find that the proton transport is dominated by structural diffusion, with the diffusion constant of the proton defect being ∼8 times higher than that of self-diffusion of the imidazole molecules. By using correlation function analysis, we decompose the mechanism for proton transport into a series of first-order processes and show that the proton transport mechanism occurs over three distinct time and length scales. Although the mechanism at intermediate times is dominated by hopping along pseudo-one-dimensional chains, at longer times the overall rate of diffusion is limited by the re-formation of these chains. These results provide a more complete picture of the traditional idealized Grotthuss structural diffusion mechanism.

Efficient computation of free energy surfaces of chemical reactions using ab initio molecular dynamics with hybrid functionals and plane waves

Ab initio molecular dynamics (AIMD) simulations employing density functional theory (DFT) and plane waves are routinely carried out using density functionals at the level of generalized gradient approximation (GGA). AIMD simulations employing hybrid density functionals are of great interest as it offers a more accurate description of structural and dynamic properties than the GGA functionals. However, the computational cost for carrying out calculations using hybrid functionals and plane wave basis set is at least two orders of magnitude higher than that using GGA functionals. Recently, we proposed a strategy that combined the adaptively compressed exchange operator formulation and the multiple time step integration scheme to reduce the computational cost by an order of magnitude [J. Chem. Phys. 151, 151102 (2019)]. In this work, we demonstrate the application of this method to study chemical reactions, in particular, formamide hydrolysis in an alkaline aqueous medium. By actuating our implementation with the well-sliced metadynamics scheme, we can compute the two-dimensional free energy surface of this reaction at the level of hybrid-DFT. This work also investigates the accuracy of the PBE0 (hybrid) and the PBE (GGA) functionals in predicting the free energetics of this chemical reaction.

TeraChem Cloud: A High-Performance Computing Service for Scalable Distributed GPU-Accelerated Electronic Structure Calculations

The encapsulation and commoditization of electronic structure arise naturally as interoperability, and the use of nontraditional compute resources (e.g., new hardware accelerators, cloud computing) remains important for the computational chemistry community. We present TeraChem Cloud, a high-performance computing service (HPCS) that offers on-demand electronic structure calculations on both traditional HPC clusters and cloud-based hardware. The framework is designed using off-the-shelf web technologies and containerization to be extremely scalable and portable. Within the HPCS model, users can quickly develop new methods and algorithms in an interactive environment on their laptop while allowing TeraChem Cloud to distribute ab initio calculations across all available resources. This approach greatly increases the accessibility of hardware accelerators such as graphics processing units (GPUs) and flexibility for the development of new methods as additional electronic structure packages are integrated into the framework as alternative backends. Cost-performance analysis indicates that traditional nodes are the most cost-effective long-term solution, but commercial cloud providers offer cutting-edge hardware with competitive rates for short-term large-scale calculations. We demonstrate the power of the TeraChem Cloud framework by carrying out several showcase calculations, including the generation of 300,000 density functional theory energy and gradient evaluations on medium-sized organic molecules and reproducing 300 fs of nonadiabatic dynamics on the B800-B850 antenna complex in LH2, with the latter demonstration using over 50 Tesla V100 GPUs in a commercial cloud environment in 8 h for approximately $1250.

Using a monomer potential energy surface to perform approximate path integral molecular dynamics simulation of ab initio water at near-zero added cost

It is now established that nuclear quantum motion plays an important role in determining water's hydrogen bonding, structure, and dynamics. Such effects are important to include in density functional theory (DFT) based molecular dynamics simulation of water. The standard way of treating nuclear quantum effects, path integral molecular dynamics (PIMD), multiplies the number of energy/force calculations by the number of beads required. In this work we introduce a method whereby PIMD can be incorporated into a DFT simulation with little extra cost and little loss in accuracy. The method is based on the many body expansion of the energy and has the benefit of including a monomer level correction to the DFT energy. Our method calculates intramolecular forces using the highly accurate monomer potential energy surface developed by Partridge-Schwenke, which is cheap to evaluate. Intermolecular forces and energies are calculated with DFT only once per timestep using the centroid positions. We show how our method may be used in conjunction with a multiple time step algorithm for an additional speedup and how it relates to ring polymer contraction and other schemes that have been introduced recently to speed up PIMD simulations. We show that our method, which we call "monomer PIMD", correctly captures changes in the structure of water found in a full PIMD simulation but at much lower computational cost.

A Versatile Multiple Time Step Scheme for Efficient ab Initio Molecular Dynamics Simulations

We present here our implementation of a time-reversible, multiple time step (MTS) method for full QM and hybrid QM/MM Born-Oppenheimer molecular dynamics simulations. The method relies on a fully flexible combination of electronic structure methods, from density functional theory to wave function-based quantum chemistry methods, to evaluate the nuclear forces in the reference and in the correction steps. The possibility of combining different electronic structure methods is based on the observation that exchange and correlation terms only contribute to low frequency modes of nuclear forces. We show how a pair of low/high level electronic structure methods that individually would lead to very different system properties can be efficiently combined in the reference and correction steps of this MTS scheme. The current MTS implementation makes it possible to perform highly accurate ab initio molecular dynamics simulations at reduced computational cost. Stable and accurate trajectories were obtained with time steps of several femtoseconds, similar to and even exceeding the ones usually adopted in classical molecular dynamics, in particular when using a generalized Langevin stochastic thermostat. Compared to the standard Velocity Verlet integration, the present MTS scheme allows for a 5- to 6-fold overall speedup, at an unaltered level of accuracy.

Computational Methods for Ab Initio Molecular Dynamics

Ab initio molecular dynamics is an irreplaceable technique for the realistic simulation of complex molecular systems and processes from first principles. This paper proposes a comprehensive and self-contained review of ab initio molecular dynamics from a computational perspective and from first principles. Quantum mechanics is presented from a molecular dynamics perspective. Various approximations and formulations are proposed, including the Ehrenfest, Born–Oppenheimer, and Hartree–Fock molecular dynamics. Subsequently, the Kohn–Sham formulation of molecular dynamics is introduced as well as the afferent concept of density functional. As a result, Car–Parrinello molecular dynamics is discussed, together with its extension to isothermal and isobaric processes. Car–Parrinello molecular dynamics is then reformulated in terms of path integrals. Finally, some implementation issues are analysed, namely, the pseudopotential, the orbital functional basis, and hybrid molecular dynamics.

Computational elucidation, mutational and hot spot-based designing of potential inhibitors against human acid-sensing ion channels (hASIC-1a) to treat various physiological conditions

Acid-sensing ion channels are ligand/proton-gated ion channels belonging to the family of the degenerin/epithelial Na⁺ channel (DEG/ENaC). They function as a sodium-selective pore for Ca²⁺ entry into neuronal cells during pathological conditions. The blocking of this channel has therapeutic importance, because at basal physiological pH (7.2), it is in a closed state and under a more acidic condition, and the ASIC1a ion channel is activated. To investigate the different states of the hASIC1a channel based on mutational analysis, structure-based virtual screening and molecular dynamics simulation studies. The system showed stability after 30 ns (after 1500 frame), and it was stabilized to an average value around 2.2Å. During the simulation, the ion channel residues in persistent contact with toxin PcTx1 were D237, E238, D347, D351, E219 and E355. These residues are important physiologically for the activation of the channel. From in silico alanine scanning, the significant hotspots obtained in hASIC1 are E344, P347, F352, D351, E355 and E219. From the sitemap analysis, it was evident that the sitemap found one of the active sites at the PcTx1 binding site with a site score of 1.086 and a D-score of 1.035 for hASIC1. We obtained a few promising hits and final potential hits from the virtual screening in hASIC1 that made interactions with the residues in the acidic pocket (E344, P347, F352, D351, E355 and E219). Based on these studies, the hits and scaffolds of potential therapeutic interest against various pathological conditions are associated with hASIC1a for future studies.

Quantum Dynamics and Spectroscopy of Ab Initio Liquid Water: The Interplay of Nuclear and Electronic Quantum Effects

Understanding the reactivity and spectroscopy of aqueous solutions at the atomistic level is crucial for the elucidation and design of chemical processes. However, the simulation of these systems requires addressing the formidable challenges of treating the quantum nature of both the electrons and nuclei. Exploiting our recently developed methods that provide acceleration by up to 2 orders of magnitude, we combine path integral simulations with on-the-fly evaluation of the electronic structure at the hybrid density functional theory level to capture the interplay between nuclear quantum effects and the electronic surface. Here we show that this combination provides accurate structure and dynamics, including the full infrared and Raman spectra of liquid water. This allows us to demonstrate and explain the failings of lower-level density functionals for dynamics and vibrational spectroscopy when the nuclei are treated quantum mechanically. These insights thus provide a foundation for the reliable investigation of spectroscopy and reactivity in aqueous environments.

Perspective: Quantum mechanical methods in biochemistry and biophysics

In this perspective article, I discuss several research topics relevant to quantum mechanical (QM) methods in biophysical and biochemical applications. Due to the immense complexity of biological problems, the key is to develop methods that are able to strike the proper balance of computational efficiency and accuracy for the problem of interest. Therefore, in addition to the development of novel ab initio and density functional theory based QM methods for the study of reactive events that involve complex motifs such as transition metal clusters in metalloenzymes, it is equally important to develop inexpensive QM methods and advanced classical or quantal force fields to describe different physicochemical properties of biomolecules and their behaviors in complex environments. Maintaining a solid connection of these more approximate methods with rigorous QM methods is essential to their transferability and robustness. Comparison to diverse experimental observables helps validate computational models and mechanistic hypotheses as well as driving further development of computational methodologies.

A Stochastic, Resonance-Free Multiple Time-Step Algorithm for Polarizable Models That Permits Very Large Time Steps

Molecular dynamics remains one of the most widely used computational tools in the theoretical molecular sciences to sample an equilibrium ensemble distribution and/or to study the dynamical properties of a system. The efficiency of a molecular dynamics calculation is limited by the size of the time step that can be employed, which is dictated by the highest frequencies in the system. However, many properties of interest are connected to low-frequency, long time-scale phenomena, requiring many small time steps to capture. This ubiquitous problem can be ameliorated by employing multiple time-step algorithms, which assign different time steps to forces acting on different time scales. In such a scheme, fast forces are evaluated more frequently than slow forces, and as the former are often computationally much cheaper to evaluate, the savings can be significant. Standard multiple time-step approaches are limited, however, by resonance phenomena, wherein motion on the fastest time scales limits the step sizes that can be chosen for the slower time scales. In atomistic models of biomolecular systems, for example, the largest time step is typically limited to around 5 fs. Previously, we introduced an isokinetic extended phase-space algorithm (Minary et al. Phys. Rev. Lett. 2004, 93, 150201) and its stochastic analog (Leimkuhler et al. Mol. Phys. 2013, 111, 3579) that eliminate resonance phenomena through a set of kinetic energy constraints. In simulations of a fixed-charge flexible model of liquid water, for example, the time step that could be assigned to the slow forces approached 100 fs. In this paper, we develop a stochastic isokinetic algorithm for multiple time-step molecular dynamics calculations using a polarizable model based on fluctuating dipoles. The scheme developed here employs two sets of induced dipole moments, specifically, those associated with short-range interactions and those associated with a full set of interactions. The scheme is demonstrated on the polarizable AMOEBA water model. As was seen with fixed-charge models, we are able to obtain large time steps exceeding 100 fs, allowing calculations to be performed 10 to 20 times faster than standard thermostated molecular dynamics.

Multiple Time-Step Dual-Hamiltonian Hybrid Molecular Dynamics - Monte Carlo Canonical Propagation Algorithm

A multiple time-step integrator based on a dual Hamiltonian and a hybrid method combining molecular dynamics (MD) and Monte Carlo (MC) is proposed to sample systems in the canonical ensemble. The Dual Hamiltonian Multiple Time-Step (DHMTS) algorithm is based on two similar Hamiltonians: a computationally expensive one that serves as a reference and a computationally inexpensive one to which the workload is shifted. The central assumption is that the difference between the two Hamiltonians is slowly varying. Earlier work has shown that such dual Hamiltonian multiple time-step schemes effectively precondition nonlinear differential equations for dynamics by reformulating them into a recursive root finding problem that can be solved by propagating a correction term through an internal loop, analogous to RESPA. Of special interest in the present context, a hybrid MD-MC version of the DHMTS algorithm is introduced to enforce detailed balance via a Metropolis acceptance criterion and ensure consistency with the Boltzmann distribution. The Metropolis criterion suppresses the discretization errors normally associated with the propagation according to the computationally inexpensive Hamiltonian, treating the discretization error as an external work. Illustrative tests are carried out to demonstrate the effectiveness of the method.

Multiple-Time Step Ab Initio Molecular Dynamics Based on Two-Electron Integral Screening

A multiple-timestep ab initio molecular dynamics scheme based on varying the two-electron integral screening method used in Hartree-Fock or density functional theory calculations is presented. Although screening is motivated by numerical considerations, it is also related to separations in the length- and timescales characterizing forces in a molecular system: Loose thresholds are sufficient to describe fast motions over short distances, while tight thresholds may be employed for larger length scales and longer times, leading to a practical acceleration of ab initio molecular dynamics simulations. Standard screening approaches can lead, however, to significant discontinuities in (and inconsistencies between) the energy and gradient when the screening threshold is loose, making them inappropriate for use in dynamics. To remedy this problem, a consistent window-screening method that smooths these discontinuities is devised. Further algorithmic improvements reuse electronic-structure information within the dynamics step and enhance efficiency relative to a naı̈ve multiple-timestepping protocol. The resulting scheme is shown to realize meaningful reductions in the cost of Hartree-Fock and B3LYP simulations of a moderately large system, the protonated sarcosine/glycine dipeptide embedded in a 19-water cluster.

Acceleration of Ab Initio QM/MM Calculations under Periodic Boundary Conditions by Multiscale and Multiple Time Step Approaches

Development of multiscale ab initio quantum mechanical and molecular mechanical (AI-QM/MM) method for periodic boundary molecular dynamics (MD) simulations and their acceleration by multiple time step approach are described. The developed method achieves accuracy and efficiency by integrating the AI-QM/MM level of theory and the previously developed semiempirical (SE) QM/MM-Ewald sum method [J. Chem. Theory Comput. 2005, 1, 2] extended to the smooth particle-mesh Ewald (PME) summation method. In the developed methods, the total energy of the simulated system is evaluated at the SE-QM/MM-PME level of theory to include long-range QM/MM electrostatic interactions, which is then corrected on the fly using the AI-QM/MM level of theory within the real space cutoff. The resulting energy expression enables decomposition of total forces applied to each atom into forces determined at the low-level SE-QM/MM method and correction forces at the AI-QM/MM level, to integrate the system using the reversible reference system propagator algorithm. The resulting method achieves a substantial speed-up of the entire calculation by minimizing the number of time-consuming energy and gradient evaluations at the AI-QM/MM level. Test calculations show that the developed multiple time step AI-QM/MM method yields MD trajectories and potential of mean force profiles comparable to single time step QM/MM results. The developed method, together with message passing interface (MPI) parallelization, accelerates the present AI-QM/MM MD simulations about 30-fold relative to the speed of single-core AI-QM/MM simulations for the molecular systems tested in the present work, making the method less than one order slower than the SE-QM/MM methods under periodic boundary conditions.