Time-propagation algorithms in molecular dynamics

In molecular dynamics simulations, the ionic positions ${\mathbf{𝐫}}_i(t)$ and velocities ${\mathbf{𝐯}}_i(t)$ are monitored as functions of time $t$. This time dependence is obtained by integrating Newton's equations of motion. When integrating the equations of motions it is important to use symplectic algorithms which conserve the phase-space volume. To solve the equations of motion under symplectic conditions, various integration algorithms have been developed. The time dependence of a particle can be expressed in a Taylor expansion

${\mathbf{𝐫}}_i(t+\mathrm{\Delta }t)={\mathbf{𝐫}}_i(t)+{\mathbf{𝐯}}_i(t)\mathrm{\Delta }t+\frac{{\mathbf{𝐅}}_i}{2m}(t)\mathrm{\Delta }{t}^2+\frac{{\partial }^3{\mathbf{𝐫}}_i(t)}{\partial t^3}\mathrm{\Delta }{t}^3+{𝒪}(\mathrm{\Delta }{t}^4)$

A backward propagation in time by a time step $\mathrm{\Delta }t$ can be obtained in a similar way

${\mathbf{𝐫}}_i(t-\mathrm{\Delta }t)={\mathbf{𝐫}}_i(t)-{\mathbf{𝐯}}_i(t)\mathrm{\Delta }t+\frac{{\mathbf{𝐅}}_i}{2m}(t)\mathrm{\Delta }{t}^2-\frac{{\partial }^3{\mathbf{𝐫}}_i(t)}{\partial t^3}\mathrm{\Delta }{t}^3+{𝒪}(\mathrm{\Delta }{t}^4)$

Adding these two equation gives and rearrangement gives the Verlet algorithm

${\mathbf{𝐫}}_i(t+\mathrm{\Delta }t)=2{\mathbf{𝐫}}_i(t)-{\mathbf{𝐫}}_i(t-\mathrm{\Delta }t)+\frac{{\mathbf{𝐅}}_i}{2m}(t)\mathrm{\Delta }{t}^2+{𝒪}(\mathrm{\Delta }{t}^3)$

The Verlet algorithm can be rearranged to the Velocity-Verlet algorithm by inserting ${\mathbf{𝐯}}_i(t)=\frac{{\mathbf{𝐫}}_i(t)-{\mathbf{𝐫}}_i(t-\mathrm{\Delta }t)}{\mathrm{\Delta }t}$

${\mathbf{𝐫}}_i(t+\mathrm{\Delta }t)={\mathbf{𝐫}}_i(t)+{\mathbf{𝐯}}_i(t)\mathrm{\Delta }t+\frac{{\mathbf{𝐅}}_i}{2m}(t)\mathrm{\Delta }{t}^2.$

Velocity-Verlet integration scheme

The Velocity-Verlet algorithm can be decomposed into the following steps:

1. ${\mathbf{𝐯}}_i(t+\frac{1}{2}\mathrm{\Delta }t)={\mathbf{𝐯}}_i(t)+\frac{{\mathbf{𝐅}}_i(t)}{2m_i}\mathrm{\Delta }t$
2. ${\mathbf{𝐫}}_i(t+\mathrm{\Delta }t)={\mathbf{𝐫}}_i(t)+{\mathbf{𝐯}}_i(t+\frac{1}{2}\mathrm{\Delta }t)\mathrm{\Delta }t$
3. compute forces ${\mathbf{𝐅}}_i(t)$ from density functional theory or machine learning
4. ${\mathbf{𝐯}}_i(t+\mathrm{\Delta }t)={\mathbf{𝐯}}_i(t+\frac{1}{2}\mathrm{\Delta }t)+\frac{{\mathbf{𝐅}}_i(t+\mathrm{\Delta }t)}{2m_i}\mathrm{\Delta }t$

From these equations it can be seen that the velocity and the position vectors are synchronous in time.

Leap-Frog integration scheme

Another form of the Verlet algorithm can be written in the form of the Leap-Frog algorithm. The Leap-Frog algorithm consists of the following steps:

1. compute forces ${\mathbf{𝐅}}_i(t)$ from density functional theory or machine learning
2. ${\mathbf{𝐯}}_i(t+\frac{1}{2}\mathrm{\Delta }t)={\mathbf{𝐯}}_i(t-\frac{1}{2}\mathrm{\Delta }t)+\frac{{\mathbf{𝐅}}_i(t)}{m_i}\mathrm{\Delta }t$
3. ${\mathbf{𝐫}}_i(t+\mathrm{\Delta }t)={\mathbf{𝐫}}_i(t)+{\mathbf{𝐯}}_i(t+\frac{1}{2}\mathrm{\Delta }t)\mathrm{\Delta }t$

In this form the velocity and the position vectors are asynchronous in time.

Thermostats and used integrators

MDALGO	thermostat	integration algorithm
0	Nosé-Hoover	Velocity-Verlet
1	Andersen	Leap-Frog
2	Nosé-Hoover	Leap-Frog
3	Langevin	Velocity-Verlet
4	Nosé-Hoover chain	Velocity-Verlet
5	CSVR	Leap-Frog
5	Multiple Andersen	Leap-Frog

Related tags and articles

IBRION, MDALGO, Thermostats