Correcting mesh-based force calculations to conserve both energy and momentum in molecular dynamics simulations

The high cost of evaluating forces in molecular dynamics makes it necessary to use approximations. The most effective approximations for nonbonded 2-body interactions, such as particle–particle particle–mesh, particle–mesh Ewald, and multilevel summation, split the potential into a short-range part which is evaluated directly and a long-range part which is interpolated from a mesh. This results in a force that, though still conservative, is not translation-invariant, resulting in drift in the linear momentum. A common remedy is to apply a uniform correction to all forces to conserve linear momentum, but this results in a force that is not conservative. Described here is a mass-weighted correction, based on the simple idea of constraining the center of mass, which does yield a conservative force. Consider a system of N particles with masses mi and positions ri = ri(t), i = 1, 2, . . . , N, which evolve in time t as dictated by Newton’s second law of motion mir→¨i=F→i(R,R˙,t), (1) where R denotes r1, r2,…, rN,. Often linear momentum Σimi ri is conserved, which is equivalent to ΣiFi = 0. In particular, a translation-invariant potential energy function U(r→1+d→,…,r→N+d→)=U(r→1,…,r→N) yields forces Fi = −∇iU(R) that sum up to zero. For a computer simulation the cost of evaluating forces can be very high. In particular, 2-body nonbonded interactions require 𝒪(N2) operations unless approximations are used. The most efficient approximations are obtained with methods such as particle–particle particle–mesh [1], particle–mesh Ewald (PME) [2], and multilevel summation [3–5], which split the potential into a short-range part calculated directly and a long-range part interpolated from a mesh. The use of a mesh yields an approximation Ua(R) to the long-range part of the potential energy which is not translation-invariant. This leads not only to a violation of Newton’s third law but causes particles to exert forces on themselves [6]! (Methods based on a hierarchical clustering of interactions, such as the fast multipole method [7], do not have this problem; however, they compute forces that are not conservative [8] and appear to be otherwise less efficient for molecular dynamics [9,4], though these studies are based on older versions of the methods.) Failure to conserve linear momentum is at best a nuisance, since it means that even if the center of mass is initially at rest, it will start moving as the simulation progresses. A common remedy [2,10,11] is to subtract out the average net force (1/N) ∑kF→ka from each force F→ia. However, the resulting force is not conservative even if the original approximate force is conservative. For a conservative force, the Jacobian matrix of Fi with respect to rj must equal the transpose of the Jacobian matrix of Fj with respect to ri. However, such is not the case for the adjusted force F→ia(R)−(1/N)∑kF→ka(R), whose Jacobian matrix with respect to position rj is given by ∂∂r→j(F→ia−1N∑kF→ka)=−ℋ↔ija+1N∑kℋ↔kja, Where ℋ↔ija is the (i, j)th 3 × 3 block of the Hessian of the potential energy function Ua. The negative effect of this simple force correction on conservation of energy is observed in [12, p. 30]. A system of 216 rigid SPC waters is simulated for 10 ns in a cube of side length 18.7 A. PME is used with a 18 × 18 × 18 grid, a direct sum cutoff of 8 A, and the Verlet integrator with a 2 fs time step. The use of the force correction produces a drift of −0.14 kcal/mol amidst fluctuations of ±0.05 kcal/mol. Without the force correction the drift is not discernible. A mass-weighted correction to the force is conservative, however. Assume that the initial linear momentum is zero: ∑imir→˙i(0)=0 (which is attained by subtracting a constant velocity from each velocity). Denote the displacement of the center of mass by g→(R)=1mtot∑kmk(r→k−r→k(0)), where mtot = Σkmk, and conservation of linear momentum is equivalent to the holonomic constraint g→(R)=0. This can be combined with the equations of motion (1) by adding constraint forces mir→¨i=F→i(…)+mimtotλ→, where λ is a set of Lagrange multipliers and its coefficient is obtained from the Jacobian matrix of g(R) with respect to ri. Eliminating λ (by twice differentiating the constraints and substituting in the equations of motion) leads to mir→¨i=F→i(…)−mimtot∑kF→k(…). A direct proof that the corrected force is conservative is obtained by defining a corrected potential energy function Uc(R)=Ua(r→1−g→(R),…,r→N−g→(R)). This is translation-invariant, and the use of a potential leads to conservative forces F→ic(r→1,…,r→N)=−∇iUc(r→1,…,r→N)=F→ia(r→1−g→(R),…,r→N−g→(R))−mimtot∑kF→ka(r→1−g→(R),…,r→N−g→(R)). (2) The assumption of zero initial linear momentum together with conservation of linear momentum implies g(R(t)) ≡ 0. Practical numerical integrators also exactly conserve linear momentum, so the corrected force simplifies to F→ic(R)=F→ia(R)−mimtot∑kF→ka(R). Extension of this technique to nonzero initial linear momentum is possible using Eq. (2), but the calculations become more complicated and time consuming. If the term g(R) is omitted, the force will not be conservative. To test the value of the mass-weighted correction, a system of 826 TIP3 water molecules in a cubic box of length 29.08 A was simulated for 10 ns. Three simulations were performed with the same (equilibrated) initial positions and velocities: one with no force correction, one with a uniform correction, and one with a mass-weighted correction. The simulation was performed with NAMD [13] using PME with cubic B-spline interpolation from a 24 × 24 × 24 grid and the Verlet integrator with a time step of 2 fs. A tolerance of 10−6 was used for the direct part of the Ewald sum, which dictates the choice of the parameter in the Ewald splitting. A switching distance of 7 A and a cutoff of 8 A was used for Lennard-Jones interactions with an XPLOR switching function (which is continuously differentiable). Rigid water constraints were solved analytically using SETTLE [14]. Figure 1 confirms that the force-corrected simulations conserve momentum whereas the uncorrected simulation terminates with a speed of 0.26 A/ps and a displacement of 1180 A for the center of mass. (The reason that the plot shows momentum not to be exactly conserved for the two force-corrected simulations is that NAMD, for parallel performance reasons, outputs the linear momentum before applying the force correction.) The values plotted in this graph, and the one that follows, are block-averaged over 50-ps intervals. Figure 2 shows total energy as a function of time. The simulations without a force correction and with a mass-weighted correction each experience a little energy drift. Theoretical and empirical evidence [15] indicates that some drift occurs for the Verlet integrator even with a conservative force and a modest stepsize. However, the simulation with the uniform correction experiences an energy drift that is an order of magnitude greater. Similar results were also obtained for two other 10-ns simulations. Hence, a simulation with the weighted correction can be run an order of magnitude longer before temperature change becomes problematic. (The Verlet integrator is symplectic for a conservative force but merely preserves volume in phase space for a nonconservative force.) Fig. 1 Speed of center of mass (A/ps) vs. time (ns) for a periodic cube of 826 TIP3P water molecules with a 2 fs time step. The solid line is for uncorrected forces, the dashed line for a mass-weighted correction, and the dotted line for a uniform correction ... Fig. 2 Time-averaged total energy (kcal/mol) vs. time (ns) for a periodic cube of 826 TIP3P water molecules with a 2 fs time step. The solid line is for uncorrected forces, the dashed line for a mass-weighted correction, and the dotted line for a uniform correction ... Implementation of a force correction on a parallel distributed-memory computer requires a reduction operation followed by a broadcast. Fortunately, the mass-weighted correction is latency tolerant—unlike the uniform correction. Practical integrators such as the Verlet scheme and its extension to constraints [16] and to multiple time stepping [17,18] incorporate force evaluations as impulses and exactly integrate these impulses. Hence force corrections (mi/mtot)λn introduced at time t = n Δ t change the positions by an amount (t – n Δ t) Δtλn/mtot for t = n Δ t, which is the same for all particles. More-over, this remains true as long as only translation-invariant forces are applied. Hence the correction to positions and velocities need not be applied until just before the next evaluation of forces that are not translation-invariant. In practice, these forces are long-range nonbonded electrostatic forces, evaluated using the largest time step. Such a technique is implemented in version 2.6 of NAMD [19].