Overcoming Stability Limitations in Biomolecular Dynamics:
Combining Force Splitting via Extrapolation with Langevin Dynamics in LN

E. Barth and T. Schlick

An efficient method for propagating biomolecular dynamics according to the Langevin equation is presented. Termed LN for its origin in a Langevin/Normal modes scheme, the method combines force linearization with force splitting techniques. The notion of linearizing the equations of motion is central to normal mode analysis, but in our context we update the harmonic model frequently (e.g., every 1 or 3 fs). The force splitting technique forms the basis of current multiple-timestepping (MTS) schemes employed for Newtonian dynamics: by updating the slowly-varying forces less frequently than the fast components, computational savings are realized without significant sacrifice in resolution accuracy. However, since the most competitive MTS versions today (r-Respa and Impulse Verlet) are formulated to be symplectic and time-reversible integrators, they are limited by resonance artifacts to an outer timestep related to half the period of the fastest motion, or around 4--5 fs, when applied to biomolecular systems. This upper bound on the frequency of updating the slow forces --- still short relative to the characteristic motion --- limits typical speedup to factors of 4--5 in comparison to explicit simulations at a timestep of 0.5 fs (the innermost timestep used).

To achieve more significant speedup factors, LN foregoes symplecticness, merging the slow and fast forces via extrapolation rather than ``impulses''. The Langevin heat bath prevents systematic energy drifts, and the resonance behavior related to the period of the fast motion is thereby lessened significantly.

We show that LN achieves excellent agreement with small-timestep solutions of the Langevin equation in terms of thermodynamics (energy means and variances), geometry, as well as dynamics (spectral densities) for two proteins in vacuum and a large water system. Significantly, the frequency of updating the slow forces extends to 48 fs or more, resulting in speedup factors well exceeding 10. The implementation of LN in any program that employs force-splitting computations is straightforward, with only partial second-derivative information required, and sparse Hessian/vector multiplication routines. The linearization part of LN could even be replaced by direct evaluation of the fast components. The application of LN to biomolecular dynamics is well suited for configurational sampling, thermodynamic, and structural questions.