Energy drift in molecular dynamics simulations (Q2458226)

Approximate solutions of Hamiltonian ordinary differential equations, specific to molecular dynamics simulations, are typically computed by the Störmer-Verlet one step method. This and other symplectic numerical methods are known to preserve the symplectic nature of the Hamiltonian flow, but they never conserve energy along trajectories. Although the system energy may be maintained in a narrow band about the true energy for quite long periods of time, there always appears a slow drift in energy away from this range. The backward error analysis allows to keep longer the energy stability. This method however fails if another source of errors, the non-smooth Hamiltonian functions which may appear through the truncation of interaction potentials, needs to be considered. The main goal of the paper is a numerical analysis of the one-dimensional test model with a non-smooth Hamiltonian, where suitable statistical assumptions allow to make precise statistical predictions about the long-term energy drift, in terms of a geometric Brownian motion. After a numerical validation of the diffusion model, a discussion is presented of its possible utility in realistic molecular systems with the Lenard-Jones interaction potential. The general methodology is expected to carry through, albeit the long-time limit stochastic processes may not be a Brownian motion nor a continuous Markov process at all.