Reference System Propagator Algorithm Research Articles

Group-based evaluation of temperature and pressure for molecular dynamics simulation with a large time step

Recently, we proposed novel temperature and pressure evaluations in molecular dynamics (MD) simulations to preserve the accuracy up to the third order of a time step, δt [J. Jung, C. Kobayashi, and Y. Sugita, J. Chem. Theory Comput. 15, 84-94 (2019); J. Jung, C. Kobayashi, and Y. Sugita, J. Chem. Phys. 148, 164109 (2018)]. These approaches allow us to extend δt of MD simulations under an isothermal-isobaric condition up to 5 fs with a velocity Verlet integrator. Here, we further improve the isothermal-isobaric MD integration by introducing the group-based evaluations of system temperature and pressure to our previous approach. The group-based scheme increases the accuracy even using inaccurate temperature and pressure evaluations by neglecting the high-frequency vibrational motions of hydrogen atoms. It also improves the overall performance by avoiding iterations in thermostat and barostat updates and by allowing a multiple time step integration such as r-RESPA (reversible reference system propagation algorithm) with our proposed high-precision evaluations of temperature and pressure. Now, the improved integration scheme conserves physical properties of lipid bilayer systems up to δt = 5 fs with velocity Verlet as well as δt = 3.5 fs for fast motions in r-RESPA, respectively.

Velocity jump processes: An alternative to multi-timestep methods for faster and accurate molecular dynamics simulations

We propose a new route to accelerate molecular dynamics through the use of velocity jump processes allowing for an adaptive time step specific to each atom-atom pair (two-body) interactions. We start by introducing the formalism of the new velocity jump molecular dynamics, ergodic with respect to the canonical measure. We then introduce the new BOUNCE integrator that allows for long-range forces to be evaluated at random and optimal time steps, leading to strong savings in direct space. The accuracy and computational performances of a first BOUNCE implementation dedicated to classical (non-polarizable) force fields are tested in the cases of pure direct-space droplet-like simulations and of periodic boundary conditions (PBC) simulations using Smooth Particle Mesh Ewald method. An analysis of the capability of BOUNCE to reproduce several condensed-phase properties is provided. Since electrostatics and van der Waals two-body contributions are evaluated much less often than with standard integrators using a 1 fs time step, up to a 400% direct-space acceleration is observed. Applying the reversible reference system propagator algorithms [RESPA(1)] to reciprocal-space (many-body) interactions allows BOUNCE-RESPA(1) to maintain large speedups in PBC while maintaining precision. Overall, we show that replacing the BAOAB standard Langevin integrator by the BOUNCE adaptive framework preserves a similar accuracy and leads to significant computational savings.

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.

Accurate Classical Polarization Solution with No Self-Consistent Field Iterations.

We present a new solution for classical polarization that does not require any self-consistent field iterations, the aspect of classical polarization that makes it computationally expensive. The new method builds upon our iEL/SCF Lagrangian scheme that defines a set of auxiliary induced dipoles whose original purpose was to serve as a time-reversible initial guess to the SCF solution of the set of real induced dipoles. In the new iEL/0-SCF approach the auxiliary dipoles now drive the time evolution of the real induced dipoles such that they stay close to the Born-Oppenheimer surface in order to achieve a truly SCF-less method. We show that the iEL/0-SCF exhibits no loss of simulation accuracy when analyzed across bulk water, low to high concentration salt solutions, and small solutes to large proteins in water. In addition, iEL/0-SCF offers significant computational savings over more expensive SCF calculations based on traditional 1 fs time step integration using symplectic integrators and is as fast as reversible reference system propagator algorithms with an outer 2 fs time step.

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.

Molecular Simulation Studies on the Elongation of Gold Nanowires in Benzenedithiol

The bonding geometry at the metal−molecule interface plays an important role in determining the conductance behavior of metal−molecule−metal junctions. This bonding geometry has to be determined a priori in quantum mechanical current−voltage (I−V) calculations. To identify the detailed metal−molecule bonding configurations, we perform classical molecular simulations by combining grand canonical Monte Carlo (GCMC) sampling with molecular dynamics (MD) to explore the dynamic elongations of gold nanowires in the presence of benzenedithiol (BDT) molecules. A specific multitime-scale double reversible reference system propagator algorithm (double-RESPA) has been designed for the metal−organic complex in MD simulations to improve the simulation efficiency. We investigate the variations of bonding sites and bonding angles of BDT molecules on a stretched Au nanowire at a constant chemical potential. The density of BDT and the number of bonded and nonbonded BDT molecules in the simulation box is monitored during t...

A reversible algorithm for Nosé molecular dynamics simulations. Equilibrium properties of liquid methane

The Trotter factorisation of the Liouville propagator is applied to derive reversible reference system propagator algorithms (RESPA). This strategy is used to treat constant temperature molecular dynamics simulations of liquid methane based on a Nosé heat bath. It is shown how the RESPA method can be used to compute thermodynamics, structural and transport properties of liquid methane. The results reported in this paper are also compared to constraint molecular dynamics simulations, quantum methods and experimental data.

Multiple time step molecular dynamics simulation for interaction between dislocations and grain boundaries

A multiple time step algorithm, called reversible reference system propagator algorithm, is introduced for the long time molecular dynamics simulation. In contrast to the conventional algorithms, the multiple time method has better convergence, stability and efficiency. The method is validated by simulating free relaxation and the hypervelocity impact of nano-clusters. The time efficiency of the multiple time step method enables us to investigate the long time interaction between lattice dislocations and low-angle grain boundaries.

Structure and Dynamics of a Benzenedithiol Monolayer on a Au(111) Surface

We use the universal force field (UFF) developed by Rappe et al. (Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024) and the specific classical potentials developed from ab initio calculations for Au−benzenedithiol (BDT) molecule interaction to perform molecular dynamics (MD) simulations of a BDT monolayer on an extended Au(111) surface. The simulation system consists of 100 BDT molecules and three rigid Au layers in a simulation box that is rhombic in the plane of the Au surface. A multiple time scale algorithm, the double-reversible reference system propagator algorithm (double RESPA) based on the Nose−Hoover dynamics scheme, and the Ewald summation with a boundary correction term for the treatment of long-range electrostatic interactions in a 2-D slab have been incorporated into the simulation technique. We investigate the local bonding properties of Au−BDT contacts and molecular orientation distributions of BDT molecules. These results show t...

Molecular Dynamics Simulations of Template-Assisted Nucleation: Alcohol Monolayers at the Air−Water Interface and Ice Formation

Perhaps the best characterized template nucleation systems are the alcohol monolayer assemblies at the air−water interface. These monolayers, consisting of aliphatic alcohols [CnH2n+1OH, n = 16−31) have been the subject of hexagonal ice (Ih) induction experiments over a range of temperatures. To learn more about the molecular basis of alcohol monolayer template−induced ice nucleation phenomena, we performed canonical molecular dynamics (MD) simulations (1 ns) of pure C29, C30, and C31 alcohol monolayer−water droplet systems, using the reversible reference system propagation algorithm (r-RESPA). We find that the MD-determined monolayer physical parameters (at 278.15 K) exhibit reasonable agreement with experimentally determined values for crystalline C30 and C31 monolayers at 278.15 K. More importantly, as the simulation temperatures approach monolayer-specific freezing points, the water layer immediately below the monolayer surfaces adopted “ice-like” lattice parameters and hexagonal or c-centered rectang...

Efficient multiple time step method for use with Ewald and particle mesh Ewald for large biomolecular systems

The particle–particle particle–mesh (P3M) method for calculating long-range electrostatic forces in molecular simulations is modified and combined with the reversible reference system propagator algorithm (RESPA) for treating the multiple time scale problems in the molecular dynamics of complex systems with multiple time scales and long-range forces. The resulting particle–particle particle–mesh Ewald RESPA (P3ME/RESPA) method provides a fast and accurate representation of the long-range electrostatic interactions for biomolecular systems such as protein solutions. The method presented here uses a different breakup of the electrostatic forces than was used by other authors when they combined the Particle Mesh Ewald method with RESPA. The usual breakup is inefficient because it treats the reciprocal space forces in an outer loop even though they contain a part that changes rapidly in time. This does not allow use of a large time step for the outer loop. Here, we capture the short-range contributions in the reciprocal space forces and include them in the inner loop, thereby allowing for larger outer loop time steps and thus for a much more efficient RESPA implementation. The new approach has been applied to both regular Ewald and P3ME. The timings of Ewald/RESPA and P3ME/RESPA are compared in detail with the previous approach for protein water solutions as a function of number of atoms in the system, and significant speedups are reported.

Time-Trimming Tricks for Dynamic Simulations: Splitting Force Updates to Reduce Computational Work

Time-Trimming Tricks for Dynamic Simulations: Splitting Force Updates to Reduce Computational Work

Langevin stabilization of molecular dynamics

In this paper we show the possibility of using very mild stochastic damping to stabilize long time step integrators for Newtonian molecular dynamics. More specifically, stable and accurate integrations are obtained for damping coefficients that are only a few percent of the natural decay rate of processes of interest, such as the velocity autocorrelation function. Two new multiple time stepping integrators, Langevin Molly (LM) and Brünger–Brooks–Karplus–Molly (BBK–M), are introduced in this paper. Both use the mollified impulse method for the Newtonian term. LM uses a discretization of the Langevin equation that is exact for the constant force, and BBK–M uses the popular Brünger–Brooks–Karplus integrator (BBK). These integrators, along with an extrapolative method called LN, are evaluated across a wide range of damping coefficient values. When large damping coefficients are used, as one would for the implicit modeling of solvent molecules, the method LN is superior, with LM closely following. However, with mild damping of 0.2 ps−1, LM produces the best results, allowing long time steps of 14 fs in simulations containing explicitly modeled flexible water. With BBK–M and the same damping coefficient, time steps of 12 fs are possible for the same system. Similar results are obtained for a solvated protein–DNA simulation of estrogen receptor ER with estrogen response element ERE. A parallel version of BBK–M runs nearly three times faster than the Verlet-I/r-RESPA (reversible reference system propagator algorithm) when using the largest stable time step on each one, and it also parallelizes well. The computation of diffusion coefficients for flexible water and ER/ERE shows that when mild damping of up to 0.2 ps−1 is used the dynamics are not significantly distorted.

Ewald summation and multiple time step methods for molecular dynamics simulation of biological molecules

Methods by which to determine conditions for a molecular dynamics (MD) simulation of biological molecules were investigated. Derivation of the optimal parameters of the Ewald summation was described so as to give same precision to the real space, the reciprocal space summations and the van der Waals interaction. Later, the procedure by which to determine the condition of the multiple time step method by RESPA (REference System Propagator Algorithm; Tuckerman et al., 1992, J. Chem. Phys., 97, 1990) was described as exemplified by MD simulations of a solvated β-sheet peptide. The conservation of the total energy in a microcanonical ensemble was measured to investigate the stability of the simulation conditions. The most feasible respective combinations of the time steps were: 0.25 fs for bond, angle and torsion interactions; 2 fs for van der Waals interaction and Ewald real-space summation; and 4 fs for Ewald reciprocal-space summation. Though it retained an acceptable accuracy, this condition accelerated the simulation ten-fold compared to that in which a simple velocity-Verlet method with a time step of 0.25 fs was used. The update of the correction term due to excluded neighbors was then investigated. Better results were obtained when the correction was updated with the real-space than when it was updated with the reciprocal-space summation. Finally, an MD simulation as long as 50 ps performed under the optimal Ewald and RESPA parameters was thus determined. The trajectory showed a good stability, indicating the feasibility of the parameters.

A molecular dynamics study of a short-chain polyethylene melt.: I. Steady-state shear

Utilizing a united atom potential model and reversible reference system propagator algorithm (rRESPA) multi-timestep dynamics, we have performed equilibrium and nonequilibrium molecular dynamics simulations of a monodisperse C 100H 202 polyethylene melt at 448 K and 0.75 g/cm 3. We report a variety of properties calculated at equilibrium including rotational relaxation time and self-diffusion coefficient as well as shear-enhanced diffusion and rheological properties calculated under steady-state shearing conditions. Shear thinning is observed in the viscosity and normal stress coefficients over the range of strain rates studied. A minimum in the hydrostatic pressure is observed at an intermediate strain rate that is associated with a minimum in the intermolecular Lennard–Jones potential energy as well as transitions in the strain-rate-dependent behavior of several other viscous and structural properties of the system. The shear field also imposes significant alignment of the chains with the flow direction, approaching a limiting angle of approximately 3° at high strain rate. In addition, the self-diffusion coefficients (calculated in terms of the unconvected positions according to the Cummings–Wang formalism) are markedly enhanced under shear compared to the equilibrium state (up to two orders of magnitude at the highest shear rate studied).

Computationally efficient canonical molecular dynamics simulations by using a multiple time‐step integrator algorithm combined with the particle mesh Ewald method and with the fast multipole method

An efficient implementation of the canonical molecular dynamics simulation using the reversible reference system propagator algorithm (r-RESPA) combined with the particle mesh Ewald method (PMEM) and with the macroscopic expansion of the fast multipole method (MEFMM) was examined. The performance of the calculations was evaluated for systems with 3000, 9999, 30,000, 60,000, and 99,840 particles. For a given accuracy, the optimal conditions for minimizing the CPU time for the implementation of the Ewald method, the PMEM, and the MEFMM were first analyzed. Using the optimal conditions, we evaluated the performance and the reliability of the integrated methods. For all the systems examined, the r-RESPA with the PMEM was about twice as fast as the r-RESPA with the MEFMM. The difference arose from the difference in the numerical complexities of the fast Fourier transform in the PMEM and from the transformation of the multipole moments into the coefficients of the local field expansion in the MEFMM. Compared with conventional methods, such as the velocity-verlet algorithm with the Ewald method, significant speedups were obtained by the integrated methods; the speedup of the calculation was a function of system size, and was a factor of 100 for a system with 3000 particles and increased to a factor of 700 for a system with 99,840 particles. These integrated calculations are, therefore, promising for realizing large-scale molecular dynamics simulations for complex systems. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 201–217, 2000

Longer time steps for molecular dynamics

Simulations of the dynamics of biomolecules have been greatly accelerated by the use of multiple time-stepping methods, such as the Verlet-I/r-RESPA (reversible reference system propagator algorithms) method, which is based on approximating “slow” forces as widely separated impulses. Indeed, numerical experiments have shown that time steps of 4 fs are possible for these slow forces but unfortunately have also shown that a long time step of 5 fs results in a dramatic energy drift. To overcome this instability, a symplectic modification of the impulsive Verlet-I/r-RESPA method has been proposed, called the mollified impulse method. The idea is that one modifies the slow part of the potential energy so that it is evaluated at “time averaged” values of the positions, and one uses the gradient of this modified potential for the slow part of the force. By filtering out excitations to the fastest motions, these averagings allow the use of longer time steps than does the impulse method. We introduce a new mollified method, Equilibrium, that avoids instability in a more effective manner than previous averaging mollified methods. Our experiments show that Equilibrium with a time step of 6 fs is as stable as the impulsive Verlet-I/r-RESPA method with a time step of 4 fs. We show that it may be necessary to include the effect of nonbonded forces in the averaging to make yet longer time steps possible. We also show that the slight modification of the potential has little effect on accuracy. For this purpose we compare self-diffusion coefficients and radial distribution functions against the Leapfrog method with a short time step (0.5 fs).