The mollified method

We have devised the mollified impulse method (MOLLY), a family of integrators [35] 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 [49]. 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 [44,46]. We have reported time steps 50% longer using Equilibrium to simulate water [49], with a real computing speedup of about 38%. Speedups of 350% are possible using mild Langevin damping while still computing correct dynamical properties [46].

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 [20,40,75,84]. After searching the literature, we implemented an efficient geometrical method to detect and update H-bonds, cf. [52,67,77,90].

We report here for the first time that H-bond MOLLY is stable for time steps of 8fs for flexible water. This represents an asymptotic twofold speedup over Verlet-I/r-RESPA. Figure 1 compares a version of MOLLY that does not include the H-bond terms, and H-bond MOLLY.

Figure: Results for 1600 fs of simulation for water at about 300 K using Verlet-I/r-RESPA, BSpline/HBond-BSpline mollified integrators using LongAverage as averaging B-Splines, all using $(\Delta t, \delta t) = (8.0; 1.0)$. The fluctuation of the total energy is averaged out by showing only the average value of every $4$ data points (spanning $32$ fs). The block-averaged total energy is then shifted to let them distinguish from each other. The integration is seen that both Verlet-I/r-RESPA and LongAverage BSpline MOLLY are unstable at $\Delta t = 8$ fs, but LongAverage HBond-BSpline MOLLY is stable at $8$ fs.

These methods still exhibit mild nonlinear instabilities and the splitting of time scales is not done very well, and thus the more powerful approach presented next is needed. These limitations are addressed in the next section.