Efficient Multiscale MOLLY Integrators for Biomolecular Simulations

Molecular dynamics for a classical unconstrained simulation requires the solution of Newton's equations of motion:

$$\mathbf{M}\frac{\displaystyle {\mathrm d}^2}{\displaystyle {\mathrm d}t^2}{X}(t) = -\nabla U(X(t))$$ (1)

where M is a diagonal matrix of atomic masses, X(t) is the collective atomic position vector, and U is the potential energy.

The Verlet-I/r-RESPA MTS impulse method splits the force into different components whose dynamics correspond to different time scales, which are then represented as appropriately weighted impulses, with weights determined by consistency. Henceforward the time dependence will not be included in the position or velocity vector. The impulse method moves the system in time according to

$$\mathbf{M}\frac{\mathrm{d^{2}}}{\mathrm{d}t\mathrm{^{2}}}X=-\... ...ox{\boldmath$\delta$ }(t-n\Delta t)\nabla U^{\mathrm{slow}}(X)$$ (2)

where the partitioning of U into U^fast and U^slow is chosen such that an appropriate time step $\Delta t$ for the slow part of the force is larger than a time step $\delta t$for the fast part.

Verlet-I/r-RESPA was proposed but not implemented by the authors of [58] and [59] and independently discovered by the authors of [148], who also demonstrated its usefulness. It permits an increase to 4fs in the length of the longest time step $\Delta t.$ When the method was introduced, it was predicted that there would occur resonances that might induce instability if the frequency of the slow force impulse coincides with a normal mode frequency of the system [20,59]. Resonance produces an oscillation in the positions whose amplitude increases with time. More surprisingly, there is also resonance for long time steps just smaller than half the period of the fastest normal mode [14,54,138]. There is also empirical evidence that time steps of 5fs or greater are not possible with this method [21].

MOLLY is a family of integrators [54] that counteracts the instabilities present in the MTS Verlet-I/r-RESPA integrator. This is accomplished by perturbing the potential using time averaged positions. It involves two distinct steps:
time averaging:

$$U^{{\rm slow}}(X)\rightarrow U^{{\rm slow}}({\mathcal{A}}(X)),$$ (3)

with the force defined as a gradient of this averaged potential, a process called
mollification:

$$-\nabla U^{{\rm slow}}(X)\rightarrow -\mathbf{{\mathcal A}}_{X}(X)^{{\rm T}}\nabla U^{ {\rm slow}}({\mathcal A}(X)),$$ (4)

where $\mathbf{{\mathcal A}}_{X}(X)$ is a matrix. This perturbation is supposed to compensate for finite $\Delta t$artifacts. The force used by MOLLY is the gradient of the perturbed potential [77]. MOLLY can be seen as a filter that eliminates components of the slow force impulse in the directions of the fast forces, and thus improves the stability of Verlet-I/r-RESPA. Different averaging functions give rise to MOLLY integrators with different stability and accuracy properties. We have developed a time averaging MOLLY method, Equilibrium*, that completely eliminates the components of $-\nabla U^{ {\rm slow}}$ in the directions of the fast forces [74,75,76]. We have reported time steps of 6fs using Equilibrium* to simulate water [77]. For a linear problem solved using Equilibrium*, the slow force time step can be chosen independently of $U^{\rm fast}$, and hence this is a true multiscale integrator, unlike Verlet-I/r-RESPA and other MOLLY integrators. This theoretical result is somewhat limited in practice, for the following reasons. For efficiency, the time averaging function is approximated by taking into account only part of the fast forces--this induces a sparse matrix in the mollification procedure. Also, unlike the model problems, MD is nonlinear.

Thus, in practical implementations of MOLLY we need to determine what forces to include in the time averaging. Besides using those that obviously are most important for stability, we have used the fact that the systems particularly sensitive to instability are those solvated in water, and that water is a hydrogen bonded system. A hydrogen bond (H-bond) is a strong long lasting nonbonded interaction that exists in several crystals and proteins. H-bonds are semi-localized in their range but may form networks. The presence of H-bonds accounts for many important properties of liquid water, proteins, DNA, and their interactions [37,66,120,139]. After searching the literature, we implemented an efficient geometrical method to detect and update H-bonds, cf. [80,109,125,147]. Simulations of water have shown that the H-bond MOLLY is superior to other MOLLY integrators, and that H-bond interactions are the most important terms to include in the time averaging [8]. We propose to extend this algorithm to handle macromolecules in an efficient way, especially a parallel method for dynamic H-bond detection (H-bonds break and form continuously in MD simulations), and to handle the sparse matrix structures resulting from H-bond updating in an efficient manner. This is in line with my experience in sparse matrix decomposition methods, cf. []. We have also exploited the special structure of the mollification filter matrix of Equilibrium* for water to make the computation faster. This is being extended to general systems []. Finally, we are applying nonlinear analysis to our methods to try to improve their stability [,,].