Basic equations and algorithms

The position of a set of particles q after a time step displacement $\Delta t$ can be obtained by a Taylor expansion

$${\bf q}_{i+1} = {\bf q}_i + \frac{d{\bf q}}{dt}\Delta t + \frac{1... ...frac{d^2{\bf q}}{dt^2}\Delta t^2 + \frac{1}{6}\frac{d^3{\bf q}}{dt^3}\Delta t^3$$ (2.101)

where in the second term in the right-hand of equation 1.101 appears the velocity ${\bf v}=\frac{d{\bf q}}{dt}$ and in the third term the acceleration ${\bf a}=\frac{d^2{\bf q}}{dt^2}=-\frac{1}{m}\frac{dV}{d{\bf q}}$

Developing the same expansion at a previous step and added to equation 1.101, with both equations truncated at third order, give place to the Verlet algorithm which is the basis of the current molecular dynamics simulations techniques

$${\bf q}_{i+1}=(2{\bf q}_i-{\bf q}_{i-1}) + {\bf a}_i(\Delta t)^2$$ (2.102)

The acceleration is taken from the derivatives of the potential energy and at the initialization point ${\bf q}_0$ the previous positions can be estimated as ${\bf q}_{i-1}={\bf q}_0 - {\bf v}_0\Delta t$. The Verlet algorithm offers several numerical problems. In addition, the fact that the velocity does not appear explictly is a problem when generating ensembles at constant temperature.

There are several improvements to the Verlet algorithm. The Leap-Frog algorithm includes the velocity in its equations

$${\bf q}_{i+1}={\bf q}_i + {\bf v}_{i+\frac{1}{2}} \Delta t$$ (2.103)

where the velocity is

$${\bf v}_{i+\frac{1}{2}}={\bf v}_{i-\frac{1}{2}} + {\bf a}_i\Delta t$$ (2.104)

velocity Verlet or the higher order predictor-corrector are other alternative integration methods. Discussion about their adequacy, the numerical stability, the energy conservation and the time-reversible character must be found in the extensive bibliography [15,50,13].