Use of higher-order integrators in molecular dynamics?

Question

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.

Answer 1

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.