Efficient molecular dynamics using geodesic integration and solvent–solute splitting

We present an approach to Langevin dynamics in the presence of holonomic constraints based on decomposition of the system into components representing geodesic flow, constrained impulse and constrained diffusion. We show that a particular ordering of the components results in an integrator that is an order of magnitude more accurate for configurational averages than existing alternatives. Moreover, by combining the geodesic integration method with a solvent–solute force splitting, we demonstrate that stepsizes of at least 8 fs can be used for solvated biomolecules with high sampling accuracy and without substantially altering diffusion rates, approximately increasing by a factor of two the efficiency of molecular dynamics sampling for such systems. The methods described in this article are easily implemented using the standard apparatus of modern simulation codes.


Introduction
Molecular dynamics (MD) is a widely used and powerful tool for studying molecular systems with extensive applications to the simulation of macromolecules both for fundamental biology/biochemistry [1,2] and, increasingly, for medical applications [3].In molecular dynamics simulation of complex systems, the most important restriction (ignoring issues of force field quality) is the size of the timestep that can be used to accurately compute trajectories.The goal of simulation is typically to unlock behaviors that occur on timescales of microceconds or more, for example state-to-state protein conformational dynamics [4] such as partial folds of proteins [5], but the use of typical molecular dynamics methods limits the timestep to a few femtoseconds.
In recent decades, with the explosion in the use of molecular dynamics, there has appeared a wide variety of integration algorithms for molecular dynamics simulation, with the typical goal being to sample the canonical ensemble defined by constant temperature, particle number and system volume.The two principal types of tools that are in widespread use for increasing the timestep in molecular sampling are multiple timestepping [6,7] and the use of holonomic constraints [8][9][10].Caution is needed when applying the former in molecular dynamics since (i) many schemes exhibit significant resonance artifacts [11,12] which substantially limit their viability for practical applications, and (ii) methods to tame the resonances may end up severely corrupting the calculation of dynamical properties such as transition and diffusion rates.
The method of constraints consists of introducing algebraic relations, for example to freeze selected bond lengths [9,10] and/or bond angles [13].It is commonly stated that the use of constraints can enable an increase in the simulation timestep for organic molecules in detailed solvent to between two and four femtoseconds, a substantial improvement on the 1-2fs typically used for fully flexible models.It is worth mentioning that a doubling of the timestep may mean halving the computational budget for a laboratory that relies largely on such MD simulations, so the cost impacts (or, conversely, the accessible range of simulation) due to increasing the timestep are considerable.The general wisdom regarding the stepsize threshold for biomolecular dynamics, with constrained stochastic dynamics, still places the limit at around 4fs [14].Efforts to go beyond this have largely involved more complicated computations (e.g.Hessian matrices and/or partial normal modes, and schemes for mollification of resonance artefacts) [15][16][17][18], often limiting their efficiency (or scalability on parallel architectures) or their suitability for implementation in major molecular modelling software.Typical methods are paired with special techniques to compute the long range electrostatic forces, such as particle mesh Ewald summation.Some extreme multiple timestepping proposals can use "outer" stepsizes (in a multiple timestepping context) of more than 100fs [19][20][21], but without any consideration of the dynamical approximation.
There is no general principle for constructing integrators for constrained stochastic systems such as Langevin dynamics.Extending an idea of Leimkuhler and Patrick [22] for deterministic dynamics on a Riemannian manifold, we propose to view the splittings of a constrained molecular dynamics problem as a collection of separate constrained problems and to solve or approximate each constrained flow.In practice, the constraints are implemented using one or more steps of the SHAKE method [8] (or its RATTLE variant [23]).As particular examples, we obtain geodesic versions of a popular splitting due to Bussi and Parinello [24] and our own "BAOAB" integrator [25,26].We find that the accuracy and stability advantages of the geodesic BAOAB scheme are so dramatic (across a wide range of stepsizes) as to render many popular alternatives uncompetitive.In the case of a biomolecular model in detailed solvent, the use of very large stepsizes leads to instabilities due to unconstrained solute modes.In this case, we show that combining the geodesic BAOAB integrator with a multiple timestepping scheme based on solvent-solute force splitting, the stepsize for biomolecules can be increased to 8 or even 9fs without significantly compromising the resolution of long timescale dynamics or sampling accuracy and with only modest additional cost per timestep.This approach is conceptually simple and straightforward to implement using existing, off-the-peg, integration components present in most MD packages (and without the need for Hessian matrices, knowledge of specific collective variables, normal mode analyses or other system-dependent analysis).

Numerical methods for constrained Langevin dynamics
The analysis of MD stability barriers has largely focussed on deterministic molecular dynamics, although some studies have emphasized the use of Langevin dynamics [15,17,24,27].The incorporation of stochastic perturbations of the force-field, combined with friction in such a way as to preserve the fluctuation-dissipation relation, dramatically changes the nature of the molecular dynamics problem.Resonances present in the deterministic system are destroyed and the ergodic nature of stochastic dynamics improves the accuracy and reliability of numerical methods, even in the long time limit.
The mathematical formulation of constrained Langevin dynamics [28] builds on the standard theory of constrained Hamiltonian systems (which may be derived starting from the constrained Euler-Lagrange equations).We write the constrained Langevin molecular dynamics system in an overdetermined form: ) where M = diag(m 1 , m 2 , . . ., m N ), with m i the mass associated to the ith degree of freedom, F = −∇U ( q ) is the vector of total force (with U : R N → R the potential energy function), γ represents a prescribed friction coefficient (or collision rate), and η(t) is a vector-valued, stationary, zero-mean Gaussian process whose components satisfy where δ ij is the Kronecker delta and δ(•) is the Dirac delta.The {g j } are smooth functions which define the constraint relationships, and the λ i are Lagrange multipliers that are chosen to maintain the constraints, m > 0 is the number of constraints, and k B T > 0 represents the temperature scaled by Boltzmann's constant.( q, p ) represents a phase space point (position and momentum vectors), restricted to the configuration manifold M = { q |g j ( q) = 0, j = 1, 2, . . ., m} and cotangent space T * M = { p |∇g j ( q ) T M −1 p = 0, j = 1, 2, . . ., m}, respectively, We assume the number of degrees of freedom N to be large, making direct sampling methods (e.g.Monte-Carlo methods) infeasible or inefficient.The cotangency constraints (2.4) are redundant, since they are automatically satisfied along solutions which lie in the configuration manifold, but it is useful to include them here for later reference.In the case of constraints, the target canonical (NVT) probability density may be written using the Dirac delta function: [28][29][30] where β = (k B T) −1 and Z is a suitable normalization constant so that the total integral of ρ g β is one.Both the equations of motion and the probability measure can alternatively be written in an intrinsic form or in a system of local coordinates, but we favor the cartesian form for simplicity.In what follows we assume that the constraints are non-degenerate, meaning that the Jacobian matrix of the constraints is of full rank at every point, and also that the stochastic constrained system is ergodic.The best existing methods for unconstrained Langevin dynamics are based on a splitting and composition technique [24,25,[31][32][33] that builds directly on the procedures used to construct symplectic integrators for Hamiltonian dynamics [34,35] The common approach is to break the stochastic vector field into parts by writing it as an additive decomposition, with each part assumed to be solvable (in the distributional sense).For example the BAOAB method [25] uses a splitting into three parts, where "B" stands for a "kick" using the force evaluated at the current position, "A" represents a linear drift (along the direction of the current velocity), and "O" is an exact (in the sense of the associated time-evolved probability distribution initiated from the point distribution at the initial conditions) solve of the Ornstein-Uhlenbeck stochastic differential equations associated to each component: The primary issue for Langevin discretization in the case of molecular dynamics is the effective perturbation of the invariant distribution induced by the numerical method , i.e. the sampling bias due to finite timestep.Although it is difficult to calculate this perturbation, it is often possible to compare, in various asymptotic limits, the perturbations and thus to select optimal methods for approximation of equilibrium averages [25,36].In addition, it is possible to analyze special cases by direct calculation, including the harmonic model problem or perturbations thereof.Of particular interest in the biomolecular dynamics setting are results which demonstrate that stepsizes of up to 2.7fs may be used safely for fully flexible, unconstrained (detailed solvent) models without degradation of configurational sampling accuracy [26] if the BAOAB splitting is used.This observation is likely the consequence of two facts: (i) BAOAB is exact for configurational sampling in the case of harmonic oscillators and, moreover, (ii) BAOAB has much smaller errors than alternative splittings for mildly anharmonic perturbations of harmonic oscillators [37].It is typically assumed that the stepsize in a Verlet molecular dynamics simulation of a fully-flexible biological model must be restricted to about 1fs to have sufficient accuracy for chemical studies, where "accuracy" is usually taken in terms of the drift in energy in long simulations.The increased stepsize usable for the stochastic MD method is not due to a change in the stability threshold but is a direct consequence of the specific choice of splitting scheme and the change of benchmark to statistical NVT configurational sampling from energy conservation.Note also that averages of kinetic energy and momentum-dependent observables computed using the BAOAB scheme are not as accurate as position-dependent averages-the scheme is effectively optimized for configurational sampling , however we hasten to add that the method is convergent for all sensible kinetic quantities and so it is possible, for example, to calculate autocorrelation functions and thus diffusion rates, as we demonstrate in Section 7. Two of the most popular constrained dynamics schemes in the literature are SHAKE and RATTLE.The SHAKE method of Ryckaert, Ciccotti and Berendsen [8] constrains the position variables to the constraint manifold g( q ) = 0 via oblique projection typically taken to be along the direction orthogonal to the manifold at the previous timestep.By contrast the "RATTLE" method (Andersen [23]) incorporates an additional orthogonal projection of the momenta to the cotangent space at a given point.For a deterministic system, the latter consists of the following sequence of calculations: A second set of multipliers {µ n j } must be computed in order to maintain the tangency conditions.RATTLE can be viewed as half a kick in the momenta, followed by oblique projection to the configuration manifold, then half a kick in the momenta and an orthogonal projection to the cotangent space (RATTLE projection).It has been shown previously that RATTLE and SHAKE are mathematically conjugate, i.e. formally equivalent from step to step under a change of variables, and both are symplectic [38].
In the case of multiple constraints, an iterative procedure must typically be used to solve the nonlinear system to implement SHAKE projection.The method of Ciccotti, Berendsen and Ryckaert [8] uses iteration over the equations (in the manner of Gauss-Seidel iteration) performing a single Newton step on each successive constraint.An alternative is to perform a Newton iteration leading to the need to solve a system of linear equations at each step [14,39].For the RATTLE projection, in the case of multiple constraints, a similar linear system must be solved.When the constraints result in a system of decoupled planar rigid bodies, the SETTLE method [40] is an exact solution of the constraints.
For constrained Langevin dynamics, typical methods are constructed by interleaving unconstrained Langevin steps with projections, usually implemented via SHAKE or RATTLE steps [28,41].In building large timestep methods, the step-project approach could, in some cases, lead to the introduction of large intermediate errors.Although the constraints are intended to be satisfied at the end of the step, there are potential problems: (i) the step may fail as there may be no solution (or more than one solution) consistent with the projection direction, and (ii) even if the step succeeds, significant error may have been introduced by the procedure, due to the oblique projection (Figure 1).

(a) (b) (c)
Figure 1: Illustration of the limitations of SHAKE projection.In SHAKE, an unconstrained step is taken from a point q n on the constraint manifold using the applied forces resulting in q unc , then a projection back to the manifold is applied along the normal direction to the constraint manifold obtained at the start of the step.(a) The SHAKE projection may not be unique, and with a large timestep, even the "nearest" solution might arise from a different part of the constraint manifold.(b) There may be no solution at all.(c) Even when there is a unique solution, significant error may arise from the oblique projection.

The geodesic integration algorithm
We now describe an integration scheme for stochastic molecular dynamics which (a) preserves the constraints and (b) is symplectic in the deterministic case (γ = 0).Although SHAKE and RATTLE are often proposed as constrained analogues of the Verlet method, we argue that the alternative approach described in this section is a more direct and geometrically coherent constrained analogue.For the deterministic case, the geodesic integrator presented here is a special case of the Riemannian manifold symplectic integrator [22] (specifically the method described in the comment at the end of Sec 4.2 of that article).We combine the scheme with efficient constraint solvers, stochastic perturbation, and, later in this article, a further splitting of the force field, in order to develop an effective method for biomolecular modelling.The idea of the geodesic integrator is straightforward.If we split the Hamiltonian as H = H 1 + H 2 , then it is desirable to solve each part, exactly preserving the constraints and tangency constraints.Splitting Hamiltonian dynamics in the natural way into kinetic and potential terms, as in the Verlet method, and introducing Lagrange multipliers to maintain the constraints, results in a pair of systems of the form We see that, for the kinetic part (A g ), the cotangency constraint is redundant, but for the system in the part defined by the potential energy (B g ), the configuration constraint is redundant (since the positions are fixed), whereas the cotangency constraint will still ensure that the changes in the momentum vector are restricted to the cotangent space.The Lagrange multipliers have been given separate names (λ, µ) to distinguish their separate roles in satisfying configurational and cotangency conditions, respectively.
The solution for part B g is easily obtained, it can be viewed as a linear projection of the unconstrained kick imparted by the forces into the cotangent space.However, the solution of the A g part is a geodesic curve (an unforced motion on the constraint surface) and so requires the solution of a nonlinear system of ODEs.A "geodesic integrator" is thus obtained by solving the A g and B g parts separately and then composing the flows.
When the constraints have a simple structure, such as planar rigid bodies (e.g.water), one may compute their geodesic flow exactly.Indeed this is exactly the framework for which the SETTLE method [40] has been developed and one finds this algorithm implemented in many production codes.Accurate methods for propagating more general rigid bodies are also available and used in gas state simulations [43].An alternative Langevin dynamics method for constrained systems of this type (i.e. which can be decomposed as a system of independent rigid bodies coupled only through forces) has been proposed by Davidchack, Oulridge and Tretyakov [44], again taking advantage of the quaternion representation to describe the constraint manifold.However our interest is ultimately in systems of a more general type for which it is typically impossible to solve the A g part exactly.Nonetheless, from an efficiency perspective we see that the geodesic integrator only requires the numerical solution of a nonlinear flow to find the drift in the configuration variables and the step does not involve any evaluation of the force field defined by the potential energy function.
In the more general situation the geodesics must be approximated by a numerical method.There are several alternative methods for computing them.The obvious approach is to employ a series of SHAKE/RATTLE steps which preserve the constraint manifold and respect the cotangency conditions.The combined method can be viewed as a sort of multiple timestepping method in which the multiple steps are not taken with fast components of the force field but merely to resolve the geodesic flow.This approach was robust and stable in our simulations, although to push the timestep up maximally in the case of a biomolecule, we needed to employ an additional level of solvent-solute decomposition as described later in this article.Higher order integrators [45] could be used to compute geodesics, but from our experiments, we doubt that this extra accuracy would be needed in practice in the biomolecular dynamics setting.
For the solution of B g we insert the equation for the acceleration into the co-tangency constraint and solve for the Lagrange multipliers.This results in a system of differential equations of the following form: The matrix Π is a projector (satisfying Π 2 = Π) onto the cotangent space defined by where G is an m × N matrix whose rows are given by the gradients of the constraint functions, i.e.
We solve the equations for q (t) and p (t) in this step of the algorithm exactly by use of the formula which is a projected "kick." In the case of Langevin dynamics, the stochastic terms (due to the constrained Ornstein-Uhlenbeck process) result in a projected SDE system of the form [28,46] These equations have the (weak, i.e. exact in the sense of averages) solution where R(t) ∼ N (0, 1) is a vector of i.i.d.normal random numbers.The computation of the exponential is easily performed by Rodrigues' formula, since For the square root, one thus obtains If we have Π p (0) = p (0) then the evolution reduces to The detailed presentation of the g-BAOAB algorithm based on these formulas is given in the Appendix.Of course other choices of the composition sequence may be adopted.As discussed earlier in this article, we have favorable experience and supporting analysis of the g-BAOAB method in the setting of unconstrained molecular dynamics [26].We also compare the g-OBABO method in the experiments.
Since the current implementation of the geodesic integrator relies on a series of SHAKE/RATTLE steps, alternatives to vanilla SHAKE, such as Newton-based techniques [14,39], LINCS [47], etc. could in principle be readily employed to implement our method.

Analysis of the sampling error
In this section we discuss the order of accuracy of the geodesic integrator for long term (ergodic) averages.Our approach to accuracy (in the sense of long term averages) is at once more transparent and more relevant for MD than the standard method based on convergence on finite time intervals ("weak order"), see, e.g.[41], which does not, in the end, typically give an asymptotic (t → ∞) result due to the difficulty of computing global bounds for the error that hold on infinite time intervals.A rigorous mathematical foundation for Langevin splitting methods for unconstrained systems, including the development of asymptotic expansions for invariant distributions and an understanding of the effect on accuracy of different ordering of the components of composition methods, is now available [36].We summarize this approach here.The analysis proceeds from study of the Fokker-Planck operator L acting on the phase space density.Under standard assumptions, the unique steady state of the Fokker-Planck equation ∂ρ/∂t = Lρ is the Gibbs density ρ β which corresponds to thermodynamic equilibrium.The splitting of the stochastic differential equations of Langevin dynamics generates a corresponding splitting of the Fokker-Planck operator.Using the Baker-Campbell-Hausdorff expansion, one then derives an effective operator whose steady states are perturbations of ρ β .That is, one obtains the following relation, sometimes referred to as a Talay-Tubaro expansion: [48] where L 1 is the leading perturbation of the operator L and f 1 describes the corresponding perturbation of the density.Then one obtains, after using Lρ β = 0, an equation for f 1 .
To analyze the error expansion for the invariant distribution in the constrained case we use the method of parameterization, whereby the constrained system is replaced by an unconstrained system in a reduced coordinate system.For a general constraint manifold this reduction can only be carried out locally and one typically needs to assemble a collection of such parameterizations to span the entire constraint manifold.However, in molecular dynamics many constraint networks we would be interested in are tree-structured [39] which implies that a single global parameterization is possible.We assume such a case here for the purpose of presenting a simplified general analysis.
After reduction to canonical coordinates θ, p θ , which follows by similar analysis to the deterministic case [35], the stochastic constrained system is equivalent to Langevin dynamics applied to a conservative system with a nonconstant mass matrix having a Hamiltonian of the form with symmetric positive definite matrix J( θ).The Langevin equations of motion become where η(t) denotes a vector of Wiener increments in the appropriate dimensional space and ∇ represents the gradient with respect to θ.The dynamics have an associated invariant distribution ρg β , where Each of the substeps of the geodesic integrator, A g , B g and O g map directly to associated parameterized integration steps Âg , Bg , Ôg : This case is not directly covered by our previous convergence analysis [36], which always assumed a constant mass matrix.Nevertheless, it is possible in principle to carry out similar calculations based on the Talay-Tubaro expansion [48].
The Fokker-Planck (forward Kolmogorov) equation gives the evolution of a distribution ρ under the exact (weak) flow, via ρ t = exp(tL)ρ 0 , where L is the adjoint of the infinitesimal generator of the vector field.The associated evolution operator can be computed for each of the A, B and O parts: where X : Y := tr(XY T ).As we have used an additive splitting of the vector field in (4.1), the overall exact operator becomes For a splitting scheme composed by solving the corresponding vector fields in sequence, the overall evolution of a distribution is characterized by the product of the exponential of these operators.For example, for the g-BAOAB scheme with stepsize δt we have exp(L Under an ergodicity assumption, we have a unique distribution such that with computed averages of an observable f (ignoring sampling error) given as If we make the ansatz that for some smooth perturbation ρ 1 , then using the Baker-Campbell-Hausdorff expansion it is useful to rewrite by virtue of a Jacobi identity, where [X, Y ] := XY − Y X is the commutator of X and Y .Thus the equation L which implies that ρ 1 ( θ, p θ ) ≡ 0 and hence implying we obtain at least second-order agreement with the exact result for the g-BAOAB scheme.Similar analysis holds for any scheme using this splitting such that its codifying string is a palindrome.However, while we expect a second order error for general observable f , the prefactor to the δt 2 term may be very different for schemes with different compositions.In the unconstrained case, our previous work [25,36] has shown that for particular families of observables some higher order terms terms can be cancelled out completely, after averaging out over the momenta.Unfortunately, preliminary calculations and numerical experiments demonstrate that such a property is lost for general symmetric positive definite J( θ); hence the g-BAOAB scheme does not have the same superconvergence properties as in the case of constant mass matrix.Nevertheless, as numerical experiments presented below demonstrate, the error in sampling using g-BAOAB is dramatically smaller than for other geodesic splitting schemes.

Solvent-solute splitting
In our simulations of biomolecules we use constraints in a relatively conservative way, freezing only the water molecules and the bond lengths between hydrogens and heavier atoms (N,C,O).
After the introduction of constraints, we find that the dominant oscillatory mode is associated to unconstrained angle bonds in the solute molecule, which is typically a polymer or protein, thus the atoms involved in the protein model need to be simulated using smaller timesteps than the atoms in the solvent bath.This is precisely the setting which motivated the development of the RESPA (multiple timestepping) algorithm [6,7], and we turn to this device to enable a longer timestep than would otherwise be feasible.
In biomolecular simulations with detailed solvent, the number of water molecules that need to be incorporated must grow in direct proportion to the volume of the periodic simulation box, whereas even if it is in a compact state, the solute never comes close to packing the box.For an alanine dipeptide simulation (22 atoms), around 500 water molecules are typically incorporated in the simulation.For larger molecules, e.g. the bovine pancreatic trypsin inhibitor (BPTI) with 1101 atoms, the total number of atoms (14281 in one representative simulation [49]) is again much greater due to the incorporation of detailed solvent.A simulation of mouse acetylcholinesterase [50] involved 8289 solute atoms and 75615 solvent atoms.Thus the number of solvent-solvent interaction forces is much larger than the number of interactions between atoms of the solute and/or between atoms of the solute and those of the solvent.Thus the computational cost of computing the solvent-solvent interactions is typically the dominant cost of the timestep.While there are ways of reducing the long-ranged force computations, they remain by far the dominant cost in typical bio-MD simulation.
In our scheme we employ RESPA in an extremely restricted form, with just two or three "fast" solute steps to each outer step.In this formulation, with a complicated reduced model which is deterministically and stochastically coupled, and a small differential in timestep, we do not anticipate or observe the resonance artifacts that can sometimes arise in other RESPA-based simulations [11,12].In the fast system we incorporate both the internal solute-solute interactions as well as the interactions between solute and solvent, since on the scale of the long stepsize the neglect of these could introduce steric clashes.In the spirit of geodesic integration, we resolve the constrained flow in tandem with the solute kicks instead of projecting only at the end of the sequence.
Let q = ( qp, qs) denote coordinates representing the position of the protein (or more generally, the solute) and solvent degrees of freedom respectively, with a consequent division of the potential energy terms as U ( q ) = Uss( qs) + Ups( qp, qs) + Upp( qp) and a corresponding splitting of the B g vector field into B g p (including forces due to Upp and Ups) and B g s involving only the solvent-solvent terms, the assumption being that computation of B g p is much cheaper than that of B g s .The evolution operator for the scheme with stepsize δt, using Kp intermediate solute steps, becomes Similar analysis to that included above shows that we would expect this scheme to be second order.In the coordinate description, one finds that [L Bg and hence when Kp = 1 the scheme reduces to g-BAOAB.Note that while we would expect stability to improve as we increase Kp, it is not necessarily true that the error is also reduced.In the unconstrained case, the BAOAB scheme has a special structure in its leading order error term that leads to a fortunate cancellation in the observed bias.This error term is modified with any changes to the scheme, so even attempting to increase the accuracy by solving some terms more accurately can lead, somewhat surprisingly, to an increase in the bias introduced to the simulation.
The stability restriction on the size of the timestep of the g-BAOAB method is a consequence of the properties of the system being simulated.In typical cases we expect this to be a biological molecule embedded in a water bath.The spectral peaks in water due to intramolecular interactions begin around 3500cm −1 (with the fastest motions related to chemical O-H, N-H and C-H bonds).These are followed by other bond stretches (e.g. the O-C-O symmetric stretch at 2400cm −1 ), other bond stretches involving carbons, and finally angle bonds arising at around 1600cm −1 [51].The intermolecular forces yield a first peak for the out-of-plane libration at 725cm −1 .If these modes were purely harmonic, the first bond stretch would represent a stepsize restriction at about 3fs, although anharmonicities alter these thresholds in numerical methods (see Schlick et al [11], and Skeel and Srinivas [52]), and the incorporation of stochastic perturbations further complicates this picture.In practice, one does in fact find that the stepsize stability threshold for liquid water is limited to about 2.7fs, or around 10% below the harmonic estimate of 3.04fs.In our simulations, we remove the fastest bonds involving hydrogen which still leaves various solute bond stretches and angles with slightly longer periods and which normally would require stepsizes (from a harmonic estimate) of around 4.5fs for stable integration.
These observations suggest that Kp = 2 or Kp = 3 interior steps of the constrained solute molecule are all that is needed to overcome the stability restrictions due to internal unconstrained angle bonds with an overall integration stepsize of up to the vicinity of 10fs; in practice (see below) we found again a reduction due to anharmonic effects, but were still able to use stepsizes of around 8-9fs.

Application of g-BAOAB to TIP3P Water
We begin with an analysis of the performance of g-BAOAB on a simple system: a box of water.For convenience, we have implemented the algorithm within the Tinker software package [53], see the supplementary material for more details.We consider a constant temperature simulation of 216 TIP3P [54] water molecules in a 18.643Å periodic box, included as one of the Tinker benchmark simulations.We run at 300K using a 9Å cutoff for all interactions, using constraints on both angles and bond lengths making each water molecule completely rigid.We compared the following: • the Midpoint-Euler-Verlet-Midpoint-Euler (MEVME) method of Lelievre, Rousset and Stoltz [28], • the scheme implemented in the Tinker standard package (v 7.1) and referred to there as the Velocity Verlet Stochastic Dynamics method (it is in fact derived from a method of Guarneri and Still [42]), and run using the command stochastic, • geodesic versions of the Bussi-Parinello method [24], denoted g-OBABO, and the geodesic BAOAB (g-BAOAB) method.
We further compare using either one or five steps of the RATTLE algorithm for each of the A pieces in the g-OBABO and g-BAOAB schemes.We run simulations using each scheme with a timestep of between 1fs and 10fs, over a fixed time period of 5ns and with the friction parameter set to a standard value 1/ps.For each experiment, results are averaged over five independent runs.Computed results are compared to the average of five baseline simulations computed using the default Tinker integrator with a timestep of 0.5fs.In Figure 2 we show how the average total potential energy varies as we change the timestep for each scheme.
It is clear from the figure that the g-BAOAB scheme is far superior for the computation of configurational quantities such as average potential energy compared to the other schemes.When using one RATTLE step per A piece (each RATTLE iteration increased wall clock time by 4% in our simulations) the timestep for the g-BAOAB scheme can be pushed as far as 9fs for TIP3P without significant error.The observed error at 9fs is comparable to the error in standard schemes Error in average potential energy The error in average total potential energy for the TIP3P model system is plotted for each scheme in comparison to the baseline result of -2078kcal/mol.The MEVME and Tinker default schemes were unable to complete a simulation past a timestep of 7fs (and then with high error), whereas the g-BAOAB scheme gives a very small error in total potential energy up to 9fs.Using more accurate geodesics improves the stability of g-OBABO, but with relatively large error.at between 2fs and 3fs.Increasing the number of RATTLE steps does not provide a significant qualitative difference to the behavior of the g-BAOAB scheme, since the sampling error is too large to further resolve the precision of the integrator and because the baseline result is computed using the Tinker default integrator which introduces substantial bias.
Interestingly, in the g-OBABO scheme, increasing the number of RATTLE steps taken improves the stability of the scheme but has little effect on the error.One would expect that converting the drift steps in the Tinker or MEVME schemes to a geodesic update using multiple RATTLE steps would improve the stability of the schemes beyond 7fs, but likely not improve the existing error results for smaller timesteps.
Figure 3 gives the computed radial distribution function for the O-H distance in the simulation, using the Tinker package's radial analysis program.We plot the curves for simulations using a timestep of 7fs, and compare to the baseline result at 0.5fs.
Although errors are present in all schemes, it is again the g-BAOAB scheme which gives the best results.The other schemes appear to heat the system artificially, pulling the G(r) curve closer to unity in the initial peak and trough.We can verify this effect by looking at the diffusion coefficient averaged over all repeated experiments, given in Table 1.
In assessing the accuracy of diffusion rates, we are faced with several issues.First, it is well known that the TIP3P model does not accurately model the self-diffusion of water [55], with an observed rate of 5.06×10 −5 cm 2 /s, around double the realistic value.In addition, Langevin dynamics reduces the diffusion compared to microcanonical calculation by an amount dependent on the magnitude of the friction coefficient [17].Our comparisons are therefore only sensible relative to an accurate discretization of Langevin dynamics with the same friction coefficient, in this case γ =1/ps.(If the key purpose of the simulation is highly accurate intermediate time dynamics, then one should either resort to a very small friction, use microcanonical simulation or else employ an alternative "gentle" thermostat [56].) The diffusion constant is approximated using Tinker's diffuse program, which computes mean squared displacements from snapshots taken every 1ps for each trajectory.Our baseline computation gives the mean diffusion constant as 4.515×10 −5 cm 2 /s, with an error of

O-H RDF computed at 7fs
Baseline (TINKER 0.5fs)  1: Computed diffusion coefficient of water (×10 −5 cm 2 /s) 0.186×10 −5 cm 2 /s.All the schemes give consistent results for the diffusion coefficient at a 3fs timestep, and the g-BAOAB scheme gives reasonably accurate results using a 9fs timestep, even when using a single RATTLE step.However, there is evident bias introduced in the other schemes, with the increased diffusion suggesting a higher temperature being sampled than that prescribed by the thermostat.
In this example, each additional RATTLE step increased the wall clock time by approximately 4%, making the g-BAOAB scheme run using five RATTLE iterations take about 36% longer compared to the standard Tinker integrator.Thus, on the basis of efficiency measured as accuracy per unit computational work, the optimal method in this case involves just a single RATTLE step in each A g solve.With a single RATTLE step the overall wall clock time is similar for all of the popular methods tested.
We note that the case of independent rigid bodies coupled through forces, such as the TIP3P model described in this section, can also be handled using the quaternionic Langevin schemes of Davidchack, Oulridge and Tretyakov [44].One of the methods ("Langevin C") of that article uses a BAOAB-inspired scheme in quaternionic representation; the authors of that article likewise observed a maximum stepsize of around 9fs with high configurational sampling accuracy, suggesting a close relationship between our scheme and theirs in this special case.However, their method does not extend to more complicated constraint manifold structures such as those associated to a protein.

Application to a biomolecular model problem
We consider a simulation of solvated alanine dipeptide, parameterized using the Amber99 force field with explicit solvent and an 8Å cutoff for all long-range interactions.The NVT simulation was conducted at 300K in a 24Å periodic cube, with 442 TIP3P water molecules kept rigid and with all bonds to hydrogen in the alanine dipeptide molecule constrained.
Testing showed that the solvent-solvent interactions were by far the most costly computation in the simulation, taking roughly 95% of the total force-call time.Hence we expect that performing the multiple timestep version of g-BAOAB will in practice be only moderately more costly than a standard constrained Langevin scheme, given an efficient implementation.(Additional discussion of efficiency is taken up in the conclusion.)Timestep (fs) Results are computed from ten independent trajectories of fixed time 5ns, using a friction of 0.5ps −1 and using one RATTLE step for the A update in g-OBABO and g-BAOAB.In Figure 4 we plot the error in the average potential energy observed in the simulation, when compared to a baseline result computed from twenty trajectories using the Tinker thermostat at a timestep of 0.5fs.
All of the standard Langevin schemes (without multiple force calls) became unstable beyond a timestep of 4fs, however the g-BAOAB scheme demonstrates an order of magnitude improvement in the observed error in the scheme.It is notable that even when run at 4fs, the g-BAOAB scheme is still over four times more accurate than conventional schemes run at half the timestep.However, increasing the number of solute calls evidently doubles the stability of the schemes, with a modest impact on the wall clock time in an efficient implementation.
When using three protein steps we were able to run stable simulations at stepsizes of up to 10fs.As this was the limit of stability for the solvent-only case (demonstrated in Section 6), we would not expect any further stability increase beyond this point.The scheme also made larger errors in average potential energy for runs with stepsize above 8fs.The mean observed error at a stepsize of 9fs was 102kcal/mol (against a total energy of around -4240kcal/mol).
An unexpected result is that using more solute steps appears to slightly increase the reported error in average potential energy.As we are effectively using a smaller timestep for the solute dynamics it would be reasonable to expect an improvement in overall accuracy compared to the more general g-BAOAB scheme.However, this result can be explained given the framework we have developed.Using more solute steps alters the sequence of Fokker-Planck operators that characterize the evolution of measure for a particular scheme.In the unconstrained dynamics, the BAOAB scheme is a special case of the five-term palindrome sequences that can be shown to exhibit special cancellations that reduce the total observed error.Changing the method alters the resulting propagation of distribution and appears to destroy the structure inherent in the operator expansion that allows these cancellations.We also compare the rate of decay of the autocorrelation function for an indicator function of the φ dihedral angle of the alanine dipeptide protein.We use the indicator function f (φ) where and compute its autocorrelation function from sixteen 5ns trajectories, generated using the standard Tinker scheme at 1fs.This indicator approximately marks a region covering one free energy well in the φ coordinate, as is shown in Figure 5.In an effort to reduce the overall sampling error we fit the averaged curves for each scheme to a decaying exponential function (using MATLAB's fit function) and compare the resulting exponents.The normalized exponents are plotted in Figure 6.Despite the noise in computing the autocorrelation function, the g-BAOAB method (both variants) gives agreement in the integrated autocorrelation time (within sampling error) up to around 8fs.With 3 solute steps, the method can be pushed to 9fs, while at 10fs a significantly poorer result is obtained, with diffusion happening much faster than expected.This is likely due to the effect of bias introduced into the distribution lowering the free energy barriers present in the system, thus changing the timescales observed in barrier crossings.
To show this effect, we estimate the height of the free energy barrier (marked in the right plot of Figure 5) for −180 < φ < −40.We plot the results in Figure 7, with the baseline result of 1.76 ± 0.01kcal/mol denoted with a dashed line.We recover a good estimate for the g-BAOAB method rspa.royalsocietypublishing.org Proc R Soc A 0000000

Normalized decay exponent for autocorrelation of φ
Figure 6: The decay of the autocorrelation function of an indicator on the φ variable is estimated by fitting the computed curve to an exponential, and then normalized based on a baseline run.
The g-BAOAB scheme, in each of the variants, gives an accurate result between 0.9 and 1.1 until 8fs.
Timestep (fs) The estimated free energy barrier height in the φ coordinate is plotted as the stepsize is varied for each scheme.The g-BAOAB method is stable and accurate to stepsizes above 8fs, giving a substantially correct effective free energy barrier height.
with stepsizes up to around 9fs.By contrast, for example the Tinker scheme performs poorly at the maximum 4fs stepsize, with an error beyond ten percent, contributing to an increased diffusion across the barrier.
In this example the constraint topology of the solute molecule is relatively simple, consisting of a number of small decoupled, tree-structured groups of atoms.The constraints are therefore quite easy to maintain and it was only necessary to use a single RATTLE step in each A piece to maintain good stability.We expect that if more bonds (or bonds and angles) are constrained, then the resolution of the constraint geometry will become more delicate and likely necessitate a more accurate approximation of the geodesics.

Conclusions
The geodesic BAOAB method offers both accuracy and stability enhancements for the biomolecule in atomistic solvent.The accuracy improvement due to the use of the g-BAOAB method (without the solute-solvent splitting) at a 4fs timestep, the largest stepsize usable by most standard constraint algorithms, is around one order of magnitude for configurational averaging (see Figure 4).This is consistent with our experience of the corresponding unconstrained scheme [26] and indicates that the g-BAOAB method is the proper analogue of BAOAB for constrained systems.
When simulations are performed in a pre-converged setting (in which sampling error, rather than perfect sampling bias, is the main obstacle) the stability rather than the accuracy of the method becomes central.By combining g-BAOAB with a multiple timestepping procedure, the constrained bio-MD stability threshold can be reliably lifted to around 8fs. Implementing the large timesteps efficiently requires separation of solvent and solute force computations and the efficient treatment of neighbor lists.In our unoptimized Tinker implementation we found that the cost of a full timestep using our 8fs (with two solute steps) g-BAOAB method was about 20% greater than a 4fs step of the standard Tinker scheme.However, we attribute much of the difference to computational work in the updating of neighbor lists (Tinker performed this operation for the full system at each intermediate force call, regardless of which atoms actually needed updating) and the small size of the alanine dipeptide system.In larger systems and with moderate optimization of the force calculation, we would expect the overhead to be no greater than around 10%.On the other hand, careful attention to some of the figures of the last section (see, e.g., Figure 7) suggests that the stepsize for the standard algorithms should be limited to less than around 3fs whereas g-BAOAB appears to be robust at 8fs, so the net efficiency gain is likely to be more than 100%.
The scheme proposed here incorporates three separate mechanisms: (i) the more accurate approximation of geodesics, (ii) the use of a solute-solvent force decomposition in RESPA-style multiple timestepping, and (iii) the fortuitous cancellation of errors in the invariant distribution due to the specific ordering of the component parts (BAOAB).It is difficult to completely separate the roles of each aspect.For example the constraints needed to be resolved once per solute step, so they are in effect more accurately treated at the large stepsize whenever solute-solvent splitting is used.We conjecture that many other integrators could be stabilized by the geodesic integration procedure and the use of solute-solvent splitting, but this only makes sense provided the sampling bias due to discretization error can be controlled.For example, using more accurate resolution of the geodesics (5 RATTLE steps) improved the stability of the g-OBABO method for TIP3P water, allowing stepsizes to be increased substantially, but with much higher sampling bias than for g-BAOAB.
Currently there is demand for algorithms that scale to much larger molecular sizes.In the case of a biomolecule in solvent, this involves more efficient treatment of the long range electrostatic forces, typically using Particle-Mesh Ewald and domain decomposition techniques.Restrictions in the way this is implemented within the Tinker package prevented us from performing simulations using Ewald summation within the solute/solvent splitting.Ewald summation leaves the essential solute/solvent dynamics in the vicinity of the biomolecule to be handled in a detailed atomistic model (i.e.our solvent-solute and geodesic integration methods would still certainly apply to this part).Regardless of the methods used, the solvent-related interactions remain the dominant part of a typical biomolecular simulation.
A natural question raised by the current work is where the next barriers lie.Our simulations constrained only the bond stretches involving hydrogen atoms.It is likely that the stepsize could be further increased (to 10fs or higher) by incorporating additional constraints on all bond stretches or even angle bonds, although there is a risk that these may in some cases alter potentially important dynamical processes.There is nothing to prevent our technique from being combined with a coloured noise thermostat [57] and other improvements [18,58] to increase the stepsize further, albeit at the cost of additional algorithmic complexity, however we note that those methods generally encounter stability thresholds above about 12fs due the out of plane libration mode and other intermolecular interactions, so breaking this barrier would likely require significantly changes in methodology.It is unclear whether the substantial additional complexity of designing a method to bridge to 12fs would be justified by the relatively slender efficiency gains.
Langevin dynamics is not the only possible choice for thermostatted MD.A recent article of Peters, Goga and Berendsen [46] addressed some other viable canonical sampling methods, including "impulsive" methods which randomize selected velocities, as well as DPD-like thermostatting schemes, also in the presence of constraints.Another family of approaches is based on "degenerate thermostats" which have a relatively mild perturbative effect on dynamical properties [56,59].Since our methods are modular and based on easily available functionalities they can also be combined with these alternative thermostats.
The methods presented here are designed to be easy to implement in standard molecular software packages.Implementation is currently being performed by the authors in the LAMMPS software package [60].Moreover, the g-BAOAB scheme is also currently being implemented in the Molecular Integrator Software Tools package of the ExTASY project [61].Using MIST, the g-BAOAB algorithm will be available in conjunction with AMBER [62] and Gromacs [63] (with support planned for many other codes).

Figure 4 :
Figure4: The error in the observed average potential energy for solvated alanine dipeptide, computed at timesteps of between 2 and 8fs.Results are averaged over ten 5ns trajectories and are to be compared to a baseline result of -4240kcal/mol for the average total potential energy.Using additional protein force evaluations for the g-BAOAB scheme (the value of Kp used in (5.1) is given in brackets) we are able to significantly increase the timestep.

Figure 5 :
Figure 5: (Left) The free energy plot in the (φ, ψ) dihedral angle space, computed from a baseline result.Contours of free energy are plotted at 0.5, 2, 3.5 and 5 kcal/mol.(Right) The marginal distribution is plotted for the φ dihedral angle.The barrier height estimated in Figure 7 is marked with the dash line.

Figure 7 :
Figure7: The estimated free energy barrier height in the φ coordinate is plotted as the stepsize is varied for each scheme.The g-BAOAB method is stable and accurate to stepsizes above 8fs, giving a substantially correct effective free energy barrier height.
Figure 3: Computed O-H radial distribution functions G(r) for each scheme at a timestep of 7fs, compared to the baseline result computed at 0.5fs (black, dashed).The initial peak and trough are highlighed in the insets.