Multiple Time Step Algorithm Research Articles

Accelerating Ab Initio Quantum Mechanical and Molecular Mechanical (QM/MM) Molecular Dynamics Simulations with Multiple Time Step Integration and a Recalibrated Semiempirical QM/MM Hamiltonian.

Molecular dynamics (MD) simulations employing ab initio quantum mechanical and molecular mechanical (ai-QM/MM) potentials are considered to be the state of the art, but the high computational cost associated with the ai-QM calculations remains a theoretical challenge for their routine application. Here, we present a modified protocol of the multiple time step (MTS) method for accelerating ai-QM/MM MD simulations of condensed-phase reactions. Within a previous MTS protocol [Nam J. Chem. Theory Comput. 2014, 10, 4175], reference forces are evaluated using a low-level (semiempirical QM/MM) Hamiltonian and employed at inner time steps to propagate the nuclear motions. Correction forces, which arise from the force differences between high-level (ai-QM/MM) and low-level Hamiltonians, are applied at outer time steps, where the MTS algorithm allows the time-reversible integration of the correction forces. To increase the outer step size, which is bound by the highest-frequency component in the correction forces, the semiempirical QM Hamiltonian is recalibrated in this work to minimize the magnitude of the correction forces. The remaining high-frequency modes, which are mainly bond stretches involving hydrogen atoms, are then removed from the correction forces. When combined with a Langevin or SIN(R) thermostat, the modified MTS-QM/MM scheme remains robust with an up to 8 (with Langevin) or 10 fs (with SIN(R)) outer time step (with 1 fs inner time steps) for the chorismate mutase system. This leads to an over 5-fold speedup over standard ai-QM/MM simulations, without sacrificing the accuracy in the predicted free energy profile of the reaction.

A multiple time step algorithm for trajectory surface hopping simulations.

A multiple time step (MTS) algorithm for trajectory surface hopping molecular dynamics has been developed, implemented, and tested. The MTS scheme is an extension of the ab initio implementation for Born-Oppenheimer molecular dynamics presented in the work of Liberatore et al. [J. Chem. Theory Comput. 14, 2834 (2018)]. In particular, the MTS algorithm has been modified to enable the simulation of non-adiabatic processes with the trajectory surface hopping (TSH) method and Tully's fewest switches algorithm. The specificities of the implementation lie in the combination of Landau-Zener and Tully's transition probabilities during the inner MTS time steps. The new MTS-TSH method is applied successfully to the photorelaxation of protonated formaldimine, showing that the important characteristics of the process are recovered by the MTS algorithm. A computational speed-up between 1.5 and 3 has been obtained compared to standard TSH simulations, which is close to the ideal values that could be obtained with the computational settings considered.

Wavefunction-Based Electrostatic-Embedding QM/MM Using CFOUR through MiMiC.

We present an interface of the wavefunction-based quantum chemical software CFOUR to the multiscale modeling framework MiMiC. Electrostatic embedding of the quantum mechanical (QM) part is achieved by analytic evaluation of one-electron integrals in CFOUR, while the rest of the QM/molecular mechanical (MM) operations are treated according to the previous MiMiC-based QM/MM implementation. Long-range electrostatic interactions are treated by a multipole expansion of the potential from the QM electron density to reduce the computational cost without loss of accuracy. Testing on model water/water systems, we verified that the CFOUR interface to MiMiC is robust, guaranteeing fast convergence of the self-consistent field cycles and optimal conservation of the energy during the integration of the equations of motion. Finally, we verified that the CFOUR interface to MiMiC is compatible with the use of a QM/QM multiple time-step algorithm, which effectively reduces the cost of ab initio MD (AIMD) or QM/MM-MD simulations using higher level wavefunction-based approaches compared to cheaper density functional theory-based ones. The new wavefunction-based AIMD and QM/MM-MD implementations were tested and validated for a large number of wavefunction approaches, including Hartree-Fock and post-Hartree-Fock methods like Møller-Plesset, coupled-cluster, and complete active space self-consistent field.

Artificial Neural Networks as Propagators in Quantum Dynamics.

The utilization of artificial neural networks (ANNs) provides strategies for accelerating molecular simulations. Herein, ANNs are implemented as propagators of the time-dependent Schrödinger equation to simulate the quantum dynamics of systems with time-dependent potentials. These ANN propagators are trained to map nonstationary wavepackets from a given time to a future time within the discrete variable representation. Each propagator is trained for a specified time step, and iterative application of the propagator enables the propagation of wavepackets over long time scales. Such ANN propagators are developed and applied to one- and two-dimensional proton transfer systems, which exhibit nuclear quantum effects such as hydrogen tunneling. These ANN propagators are trained for either a specific time-independent potential or general potentials that can be time-dependent. Hierarchical, multiple time step algorithms enable parallelization, and the extension to higher dimensions is straightforward. This strategy is applicable to quantum dynamical simulations of diverse chemical and biological processes.

Numerical simulations of liquid-solid flows with free surface by coupling IMPS and DEM

Liquid-solid two-phase flow with free surface is a significant issue in nature and engineering fields. In present study, an in-house solver MPSDEM-SJTU based on a fully Lagrangian coupled method is developed to solve the problems of liquid-solid two-phase flows. The improved moving particle semi-implicit (IMPS) method is applied for the simulation of incompressible viscous liquid flows while the discrete element method (DEM) is used to model the interaction among the solid particles. The coupling strategy is based on the local averaging technique. Multiple time-step algorithm is adopted for those two methods to balance the stability and the efficiency. The free surface detection and the solid-fluid neighbor search are modified. A simulation of multi-balls collapse is conducted to verify the DEM program of the present solver. The distributions of solid balls at different instants are in good agreement with the experimental results. Then, this solver is applied for two typical problems of two-phase flows, including the two-phase dam-break and the liquid-solid flows in a rotating cylindrical tank. The flow front, flow patterns and the shape of solid beds all agree well with experimental observations.

Hamiltonian based resonance-free approach for enabling very large time steps in multiple time-scale molecular dynamics

Extended phase-space isokinetic methods in their deterministic [Minary et al., Phys. Rev. Lett. 93, 150201 (2004)] and stochastic forms [Leimkuhler et al., Mol. Phys. 111, 3579 (2013)] have proved tremendously successful in allowing multiple time-scale molecular dynamics simulations to be performed with very large time steps. These methods work by coupling the physical degrees of freedom to a set of Nosé-Hoover chain or Nosé-Hoover Langevin thermostats via an isokinetic constraint, which has the effect of avoiding resonance artifacts that plague multiple time-step algorithms. In this paper, we introduce a new resonance-free approach that achieves the same gains in time step but without the imposition of isokinetic constraints or the introduction of extended phase-space variables. Rather, we modify the physical Hamiltonian that effects the same regulation of resonances achieved by the isokinetic constraints. In so doing, we show that sampling errors can be controlled and performance improvements are possible within a simpler Hamiltonian framework. The method is demonstrated in simulations of the structure of liquid water and, in conjunction with enhanced sampling, in generation of the Ramachandran free-energy surface of the solvated alanine dipeptide.

On the use of multiple-time-step algorithms to save computing effort in molecular dynamics simulations of proteins.

Computer simulation of proteins in aqueous solution at the atomic level of resolution is still limited in time span and system size due to limited computing power available and thus employs a variety of time‐saving techniques that trade some accuracy against computational effort. Examples of such time‐saving techniques are the application of constraints to particular degrees of freedom or the use of a multiple‐time‐step (MTS) algorithm distinguishing between particular forces when integrating Newton's equations of motion. The application of two types of MTS algorithms to bond‐stretching forces versus the remaining forces in molecular dynamics (MD) simulations of a protein in aqueous solution or of liquid water is investigated and the results in terms of total energy conservation and the influence on various other properties are compared to those of MD simulations of the same systems using bond‐length, and for water bond‐angle, constraints. At comparable computational effort, the use of bond‐length constraints in proteins leads to better energy conservation and less distorted properties than the two MTS algorithms investigated.

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations.

Background:Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts. Methods:The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution. Results:In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces for weak coupling and Nosé-Hoover (chain) thermostats. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions or a stochastic dynamics thermostat, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs. Conclusion:TR schemes are not advisable to use in combination with weak coupling or Nosé-Hoover (chain) thermostats due to the occurrence of significant numerical artifacts at interfaces.

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations.

Background:Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts. Methods:The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution. Results:In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces for weak coupling and Nosé-Hoover (chain) thermostats. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions or a stochastic dynamics thermostat, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs. Conclusion:TR schemes are not advisable to use in combination with weak coupling or Nosé-Hoover (chain) thermostats due to the occurrence of significant numerical artifacts at interfaces.

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

Background: Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts. Methods: The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution. Results: In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces for weak coupling and Nosé-Hoover (chain) thermostats. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions or a stochastic dynamics thermostat, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs. Conclusion: TR schemes are not advisable to use in combination with weak coupling or Nosé-Hoover (chain) thermostats due to the occurrence of significant numerical artifacts at interfaces.

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

Background:Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts.Methods:The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution.Results:In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces. The presence of the observed artefacts was found to be independent of the interface shape and the thermostatting method used. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs.Conclusion:A SR scheme avoids numerical artefacts and offers an interesting alternative to TR RESPA with respect to performance optimization.

Molecular dynamics based enhanced sampling of collective variables with very large time steps.

Enhanced sampling techniques that target a set of collective variables and that use molecular dynamics as the driving engine have seen widespread application in the computational molecular sciences as a means to explore the free-energy landscapes of complex systems. The use of molecular dynamics as the fundamental driver of the sampling requires the introduction of a time step whose magnitude is limited by the fastest motions in a system. While standard multiple time-stepping methods allow larger time steps to be employed for the slower and computationally more expensive forces, the maximum achievable increase in time step is limited by resonance phenomena, which inextricably couple fast and slow motions. Recently, we introduced deterministic and stochastic resonance-free multiple time step algorithms for molecular dynamics that solve this resonance problem and allow ten- to twenty-fold gains in the large time step compared to standard multiple time step algorithms [P. Minary et al., Phys. Rev. Lett. 93, 150201 (2004); B. Leimkuhler et al., Mol. Phys. 111, 3579-3594 (2013)]. These methods are based on the imposition of isokinetic constraints that couple the physical system to Nosé-Hoover chains or Nosé-Hoover Langevin schemes. In this paper, we show how to adapt these methods for collective variable-based enhanced sampling techniques, specifically adiabatic free-energy dynamics/temperature-accelerated molecular dynamics, unified free-energy dynamics, and by extension, metadynamics, thus allowing simulations employing these methods to employ similarly very large time steps. The combination of resonance-free multiple time step integrators with free-energy-based enhanced sampling significantly improves the efficiency of conformational exploration.

Evaluation of atomic pressure in the multiple time-step integration algorithm.

In molecular dynamics (MD) calculations, reduction in calculation time per MD loop is essential. A multiple time-step (MTS) integration algorithm, the RESPA (Tuckerman and Berne, J. Chem. Phys. 1992, 97, 1990-2001), enables reductions in calculation time by decreasing the frequency of time-consuming long-range interaction calculations. However, the RESPA MTS algorithm involves uncertainties in evaluating the atomic interaction-based pressure (i.e., atomic pressure) of systems with and without holonomic constraints. It is not clear which intermediate forces and constraint forces in the MTS integration procedure should be used to calculate the atomic pressure. In this article, we propose a series of equations to evaluate the atomic pressure in the RESPA MTS integration procedure on the basis of its equivalence to the Velocity-Verlet integration procedure with a single time step (STS). The equations guarantee time-reversibility even for the system with holonomic constrants. Furthermore, we generalize the equations to both (i) arbitrary number of inner time steps and (ii) arbitrary number of force components (RESPA levels). The atomic pressure calculated by our equations with the MTS integration shows excellent agreement with the reference value with the STS, whereas pressures calculated using the conventional ad hoc equations deviated from it. Our equations can be extended straightforwardly to the MTS integration algorithm for the isothermal NVT and isothermal-isobaric NPT ensembles. © 2017 Wiley Periodicals, Inc.

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.

Efficient multiple time-stepping algorithms of higher order

Multiple time-stepping (MTS) algorithms allow to efficiently integrate large systems of ordinary differential equations, where a few stiff terms restrict the timestep of an otherwise non-stiff system. In this work, we discuss a flexible class of MTS techniques, based on multistep methods. Our approach contains several popular methods as special cases and it allows for the easy construction of novel and efficient higher-order MTS schemes. In addition, we demonstrate how to adapt the stability contour of the non-stiff time-integration to the physical system at hand. This allows significantly larger timesteps when compared to previously known multistep MTS approaches. As an example, we derive novel predictor–corrector (PCMTS) schemes specifically optimized for the time-integration of damped wave equations on locally refined meshes. In a set of numerical experiments, we demonstrate the performance of our scheme on discontinuous Galerkin time-domain (DGTD) simulations of Maxwell's equations.

A multiple time stepping algorithm for efficient multiscale modeling of platelets flowing in blood plasma

We developed a multiple time-stepping (MTS) algorithm for multiscale modeling of the dynamics of platelets flowing in viscous blood plasma. This MTS algorithm improves considerably the computational efficiency without significant loss of accuracy. This study of the dynamic properties of flowing platelets employs a combination of the dissipative particle dynamics (DPD) and the coarse-grained molecular dynamics (CGMD) methods to describe the dynamic microstructures of deformable platelets in response to extracellular flow-induced stresses. The disparate spatial scales between the two methods are handled by a hybrid force field interface. However, the disparity in temporal scales between the DPD and CGMD that requires time stepping at microseconds and nanoseconds respectively, represents a computational challenge that may become prohibitive. Classical MTS algorithms manage to improve computing efficiency by multi-stepping within DPD or CGMD for up to one order of magnitude of scale differential. In order to handle 3–4 orders of magnitude disparity in the temporal scales between DPD and CGMD, we introduce a new MTS scheme hybridizing DPD and CGMD by utilizing four different time stepping sizes. We advance the fluid system at the largest time step, the fluid-platelet interface at a middle timestep size, and the nonbonded and bonded potentials of the platelet structural system at two smallest timestep sizes. Additionally, we introduce parameters to study the relationship of accuracy versus computational complexities. The numerical experiments demonstrated 3000x reduction in computing time over standard MTS methods for solving the multiscale model. This MTS algorithm establishes a computationally feasible approach for solving a particle-based system at multiple scales for performing efficient multiscale simulations.

A multiple time step scheme for multiresolved models of Macromolecules

In hybrid particle models where coarse-grained beads and atoms are used simultaneously, two clearly separate time scales are mixed. If such models are used in molecular dynamics simulations, a multiple time step (MTS) scheme can therefore be used. In this manuscript, we propose a simple MTS algorithm which approximates for a specific number of integration steps the slow coarse-grained bead-bead interactions with a Taylor series approximation while the atom-atom ones are integrated every time step. The procedure is applied to a previously developed hybrid model of a melt of atactic polystyrene (di Pasquale, Marchisio, and Carbone, J. Chem. Phys. 2012, 137, 164111). The results show that structure, local dynamics, and free diffusion of the model are not altered by the application of the integration scheme which can confidently be used to simulate multiresolved models of polymer melts.

Stochastic, resonance-free multiple time-step algorithm for molecular dynamics with very large time steps

Molecular dynamics is one of the most commonly used approaches for studying the dynamics and statistical distributions of physical, chemical, and biological systems using atomistic or coarse-grained models. It is often the case, however, that the interparticle forces drive motion on many time scales, and the efficiency of a calculation is limited by the choice of time step, which must be sufficiently small that the fastest force components are accurately integrated. Multiple time-stepping algorithms partially alleviate this inefficiency by assigning to each time scale an appropriately chosen step-size. As the fast forces are often computationally cheaper to evaluate than the slow forces, this results in a significant gain in efficiency. However, such approaches are limited by resonance phenomena, wherein motion on the fastest time scales limits the step sizes associated with slower time scales. In atomistic models of biomolecular systems, for example, resonances limit the largest time step to around 5-6 fs. Stochastic processes promote mixing and ergodicity in dynamical systems and reduce the impact of resonant modes. In this paper, we introduce a set of stochastic isokinetic equations of motion that are shown to be rigorously ergodic, largely free of resonances, and can be integrated using a multiple time-stepping algorithm which is easily implemented in existing molecular dynamics codes. The technique is applied to a simple, illustrative problem and then to a more realistic system, namely, a flexible water model. Using this approach outer time steps as large as 100 fs are shown to be possible.

Dual Effect of Phosphatidyl (4,5)-Bisphosphate PIP2 on Shaker K+ Channels

Dual Effect of Phosphatidyl (4,5)-Bisphosphate PIP2 on Shaker K+ Channels