next up previous contents
Next: Predictor-corrector algorithm Up: Time integration algorithm Previous: Time integration algorithm

The Verlet algorithm

  In molecular dynamics, the most commonly used time integration algorithm is probably the so-called Verlet algorithm [13]. The basic idea is to write two third-order Taylor expansions for the positions ${\bf r} (t)$, one forward and one backward in time. Calling $\bf v$ the velocities, $\bf a$ the accelerations, and $\bf b$ the third derivatives of ${\bf r}$ with respect to t, one has:

\begin{displaymath}
{\bf r} (t+\Delta t) = {\bf r} (t) + {\bf v} (t) \Delta t
 +...
 ...(t) \Delta t^2 + (1/6) {\bf b} (t) \Delta t^3
 + O(\Delta t^4) \end{displaymath}

\begin{displaymath}
{\bf r} (t-\Delta t) = {\bf r} (t) - {\bf v} (t) \Delta t
 +...
 ...(t) \Delta t^2 - (1/6) {\bf b} (t) \Delta t^3
 + O(\Delta t^4) \end{displaymath}

Adding the two expressions gives
\begin{displaymath}
{\bf r} (t+\Delta t) = 2{\bf r} (t) - {\bf r} (t-\Delta t)
 + {\bf a} (t) \Delta t^2 + O(\Delta t^4) \end{displaymath} (7)
This is the basic form of the Verlet algorithm. Since we are integrating Newton's equations, ${\bf a} (t)$ is just the force divided by the mass, and the force is in turn a function of the positions ${\bf r} (t)$:
\begin{displaymath}
{\bf a} (t) = - (1/m) {\bf\nabla} V\left( {\bf r}(t) \right) \end{displaymath} (8)

As one can immediately see, the truncation error of the algorithm when evolving the system by $\Delta t$ is of the order of $\Delta t^4$, even if third derivatives do not appear explicitly. This algorithm is at the same time simple to implement, accurate and stable, explaining its large popularity among molecular dynamics simulators.

A problem with this version of the Verlet algorithm is that velocities are not directly generated. While they are not needed for the time evolution, their knowledge is sometimes necessary. Moreover, they are required to compute the kinetic energy K, whose evaluation is necessary to test the conservation of the total energy E=K+V. This is one of the most important tests to verify that a MD simulation is proceeding correctly. One could compute the velocities from the positions by using

\begin{displaymath}
{\bf v} (t) = \frac { {\bf r}(t+\Delta t) - {\bf r}(t-\Delta t) }
 { 2 \Delta t } . \end{displaymath}

However, the error associated to this expression is of order $\Delta t^2$ rather than $\Delta t^4$.

To overcome this difficulty, some variants of the Verlet algorithm have been developed. They give rise to exactly the same trajectory, and differ in what variables are stored in memory and at what times. The leap-frog algorithm, not reported here, is one of such variants [20] where velocities are handled somewhat better.

An even better implementation of the same basic algorithm is the so-called velocity Verlet scheme, where positions, velocities and accelerations at time $t+\Delta t$ are obtained from the same quantities at time t in the following way:

\begin{displaymath}
{\bf r} (t + \Delta t) = {\bf r} (t) + {\bf v} (t) \Delta t +
 (1/2) {\bf a} (t) \Delta t^2 \end{displaymath}

\begin{displaymath}
{\bf v} (t + \Delta t/2) = {\bf v} (t) + (1/2) {\bf a} (t) \Delta t \end{displaymath}

\begin{displaymath}
{\bf a} (t + \Delta t) = - (1/m) {\bf\nabla} 
 V \left( {\bf r}(t+\Delta t) \right) \end{displaymath}

\begin{displaymath}
{\bf v} (t + \Delta t) = {\bf v} (t + \Delta t/2) +
 (1/2) {\bf a} (t + \Delta t) \Delta t \end{displaymath}

Note how we need 9N memory locations to save the 3N positions, velocities and accelerations, but we never need to have simultaneously stored the values at two different times for any one of these quantities.


next up previous contents
Next: Predictor-corrector algorithm Up: Time integration algorithm Previous: Time integration algorithm
Furio Ercolessi
9/10/1997