# Use of higher-order integrators in molecular dynamics?

The usual verlet integration scheme for propagating molecular dynamics according to Newton's equations of motion looks like the following: $$x(t+\Delta t)=2x(t)-x(t-\Delta t)+a(t)\Delta t^2$$

This is accurate up to $$O(\Delta t^2)$$. There are, however, higher-order integrators which make a smaller discretization error and depend upon the third time-derivative of $$x(t)$$, the so-called jerk. I am wondering if the use of these higher-order integrators have been used at all, especially for relatively simple potentials where one could presumably calculate the jerk term easily?

Also, how would one actually go about calculating the third derivative of position with respect to time? For instance, in practice we actually calculate $$a(t)=-\nabla V(x(t))/m$$. That is, we calculate the gradient of the potential, taking advantage of Newton's second law. Is there any analogous relationship for the third-derivative of $$x(t)$$? Is it possibly related to the Hessian? Or would we be stuck using some kind of finite-difference scheme to get the third-derivative, which would presumably be very slow and thus invalidate of using a more accurate integrator, hence allowing us to take larger time steps.

• Note that the jerk term disappears in your first equation since odd powers cancel out. I'm not sure about higher-order terms (e.g. where we use $a(t)\Delta t^2+c(t)\Delta t^4+\mathcal O(\Delta t^6)$) since $\Delta t$ is expected to be quite small anyway. It might potentially make sense if the problem at hand requires combining a significant number of uncertainties though. Nov 17, 2021 at 9:32
• @TheSimpliFire True. I overlooked that. I guess that's probably why people don't use higher-order integrators for MD. It results in needing to calculate the fourth-derivative with respect to time. Nonetheless, I'd be curious if people have tried this at all. Nov 17, 2021 at 18:14
• higher order methods are often more accurate over short times, but not over long times. Frenkel & Smit as well as Allen & Tildesley cover integrators fairly well. It is better to shadow the real trajectory, on average, well over a long time, than be very accurate for a short while but start diverging over a long time Nov 24, 2021 at 18:56

Higher-order integrators are used, but usually the way they perform the calculation is not through directly calculating higher order derivatives, but essentially through multiple force calculations. The Wikipedia page on symplectic integrators gives some information on this, including 3rd order and 4th order examples. I've personally used the 4th order integrator in semiclassical simulations, where we were concerned about high accuracy to preserve quantum phase interference effects (although I'm not sure it was actually necessary).

In condensed phase systems, we're often simulating at a constant temperature, frequently with a stochastic thermostat, so the error in question should be how well we preserve the ensemble, not the energy itself. The "BAOAB" Langevin splitting has been shown to have a coincidental cancellation of errors that makes it a higher order integrator at the cost of a lower order scheme.

Non-symplectic predictor-corrector schemes used to be used more frequently, but it's generally agreed that symplectic methods are superior (thanks to some guaranteed properties coming from them being canonical transformations). I've used 4th-order Adam-Bashforth-Moulton predictor-correctors, as well as a 6th order Gear predictor corrector (modified to use adaptive time steps; needed when I hit a system where fixed time steps didn't work). The Gear predictor-corrector is presented in an appendix of [EDIT: the first edition of] Allen and Tildesley's "Computer Simulation of Liquids," which for a long time was one of the best books to learn MD (it's a little out of date now, being as it was published in 1987, [EDIT: and it looks like the new edition quite reasonably removed this appendix]). The Gear PC tries to actually estimate the values of the derivatives (whereas ABM uses a higher discretization in time). Some of my own old notes on ABM integrator and 4th order symplectic integrators can be found under the heading "Qual Prep" here (warning -- you'll have to decipher my handwriting).

Overall, although higher-order integrators have certainly been used, I think the overall sense in the community is that the cost-benefit ratio just isn't worth it.

Also, how would one actually go about calculating the third derivative of position with respect to time?

This should be feasible just by applying the chain rule. However, I think this would indeed require the whole Hessian matrix, which rules it out for large-scale MD.

• There's a second edition of Allen and Tildesley - my copy was published in 2017. I haven't been through it systematically, but there certainly has been a lot of new and updated material added. Nov 19, 2021 at 9:14
• Thanks for the info, @IanBush! From the contents pages available on Amazon's "Look inside", it appears the second edition removed the appendix I was referring to -- which makes sense, since the Gear predictor-corrector is not widely used in modern MD. I updated the answer accordingly. Nov 19, 2021 at 23:17
• You have hinted at it by mentioning "some guaranteed properties" of symplectic integrators, but it might be worth spelling out for beginners: Despite their high order, Gear methods have terrible energy conservation because they are not time reversible. Symplectic methods like Verlet conserve energy almost exactly. Apr 19 at 8:06