The Verlet algorithm

Page courtesy of Prof.Furio Ercolessi at SISSA

  In molecular dynamics, the most commonly used time integration algorithm is probably the so-called Verlet algorithm (L. Verlet, Phys. Rev. 159, 98 (1967); Phys. Rev. 165, 201 (1967). ) The basic idea is to write two third-order Taylor expansions ${\bf r} (t)$, one forward and one backward in time.

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)$:

$${\bf a} (t) = - (1/m) {\bf\nabla} V\left( {\bf r}(t) \right)$$

The simplest form, using the positions, is explained here. However, we will not use this form in CompMeth.

(see e.g. The Feynman Lectures on Physics, Vol. 1, Addison-Wesley, 1963, Chapter 9 ``Newton's Laws of Dynamics'')

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:

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

$${\bf v} (t + \Delta t/2) = {\bf v} (t) + (1/2) {\bf a} (t) \Delta t$$

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

$${\bf v} (t + \Delta t) = {\bf v} (t + \Delta t/2) + (1/2) {\bf a} (t + \Delta t) \Delta t$$

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.

Choosing the timestep dt is essential to success. It should be at least an order of magnitude less than the "typical" times of the system (as defined by phonon frequencies or ratio of velocity to acceleration). Too large, and errors will accrue in the integration. Too small, and errors will occur from rounding in the computation.

Here are simple codes which use velocity Verlet to integrate a falling particle and Verlet with xmgrace to illustrate a bouncing particle

Your code will require a loop over time, with the velocity and