Molecular dynamics is a natural tool to explore the conformational space of molecules, particularly proteins and other biomolecules [59, p. 434]. A conformation is a particular geometry of the molecule. Stable geometries are believed to correspond to an energy minimum [26]. For example, MD simulations permit ligands and receptors to explore their configurations in space to determine the efficiency of molecular docking. Study of the dynamics itself is important since the biochemical functionality of some complex systems such as enzymes depend on dynamic processes like gating or side-chain rotations in receptor molecules [66,101]. MD is also commonly used to study thermodynamic properties of systems [1, p. 46].
Molecular dynamics (MD) solves Newton's equations of motion by evaluating
pairwise force interactions between particles, numerically integrating the
system in time, and imposing boundary conditions, cf.
[9,34,38,54,78]. This is a challenging problem not
only because the number of particles is large but also because stability
severely limits the length of the time steps relative to the total length
needed for simulations--time steps are in the order of femtoseconds
(
seconds) whereas simulations of a few microseconds (
seconds) up to one second are most desired for processes such as protein
folding. The time scales span large ranges, about 15 orders of magnitude.
Multiple time stepping (MTS) integrators such as Verlet-I/r-RESPA have been used to lengthen the time step for most of the interactions in the equations of motion, cf. [34,91]. These methods evaluate different parts of the force at different frequencies. Limitations on the step size in MTS integrators are still severe, and these are mostly due to stability rather than accuracy.
We have reported 50% speedup over Verlet-I/r-RESPA by using a mollified version that perturbs the potential eliminating sources of instability [49], and speedups of 350% using mild Langevin damping while still computing correct dynamical properties [46]. Recent analysis and experiments have revealed more clearly what is necessary to remove remaining instabilities [47,60]. Using results from this analysis we have devised stabler MTS integrators that already yield speedups from 2 to 4 over Verlet-I/r-RESPA. See Sections 2.1.1 and 3.1.
The approach of Elber and coworkers [71] has been to eliminate high
frequencies altogether to an arbitrary point. They accomplish this by formulating
MD as a boundary value problem and using implicit methods to produce a
``reaction path'' to arbitrary accuracy. This method has produced impressive
outcomes, such as recent results of a millisecond simulation of an enzyme
[93]. However, the applicability of this method is limited by the
need to have an initial and a final configuration. Our proposed method uses
symplectic semi-implicit methods to create an approximate MD simulation that
is not limited by instability. Also, our method does not completely eliminate
fast motions but incorporates their effect in an approximate manner. Speedups
of several orders of magnitude are possibleover current MD technology.
See Sections 2.1.4 and 3.2.