Verlet Algorithm

In molecular dynamics, the most commonly used time integration algorithm is probably the so-called Verlet algorithm [L. Verlet, Computer experiments on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules, Physical Review 159, 98 (1967)]. 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{equation} {\bf r} (t+\Delta t) = {\bf r} (t) + {\bf v} (t) \Delta t + \frac{1}{2} {\bf a}(t) \Delta t^2 + (1/6) {\bf b} (t) \Delta t^3

O(\Delta t^4) \end{equation} $$ $$ \begin{equation} {\bf r} (t-\Delta t) = {\bf r} (t) - {\bf v} (t) \Delta t + \frac{1}{2} {\bf a}(t) \Delta t^2 - (1/6) {\bf b} (t) \Delta t^3
O(\Delta t^4) \end{equation} $$ Adding the two expressions gives $$\begin{equation} {\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{equation}
$$ 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{equation} {\bf a} (t) = - \frac{1}{m} {\bf\nabla} V\left( {\bf r}(t) \right) \end{equation}
$$ 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.

While the velocities 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 subtracting the previous expression to obtain:

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

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

The main problem with the Verlet algorithm is that it is not self starting, and the first step needs to be computed by different means. An additional problem is that the new velocity is found by computing the difference between two quantities of the same order of magnitude. When using computers which always operate with finite numerical precision, such an operation results in a loss of numerical precision and may give rise to substantial roundoff error.

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{eqnarray} {\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 \end{eqnarray}$$

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.

Here, we modify the code for particle2 implementing velocity Verlet:



In [2]:

    
class particle2(object):
    
    def __init__(self, mass=1., x=0., y=0., vx=0., vy=0.):
        self.mass = mass
        self.x = x
        self.y = y
        self.vx = vx
        self.vy = vy
       
    def euler(self, fx, fy, dt):
        self.vx = self.vx + fx*dt
        self.vy = self.vy + fy*dt
        self.x = self.x + self.vx*dt
        self.y = self.y + self.vy*dt
        
    def get_force(self):  # returns force per unit of mass (acceleration)
        GM=4*math.pi*math.pi # We use astronomical units
        r = math.sqrt(self.x*self.x+self.y*self.y)
        r3 = r * r * r
        fx = -GM*self.x/r3
        fy = -GM*self.y/r3
        return (fx,fy)
        
    def verlet(self, dt):
        (fx,fy) = self.get_force() # before I move to the new position
        self.x += self.vx*dt + 0.5*fx*dt*dt
        self.y += self.vy*dt + 0.5*fy*dt*dt
        self.vx += 0.5*fx*dt
        self.vy += 0.5*fy*dt
        (fx,fy) = self.get_force() # after I move to the new position
        self.vx += 0.5*fx*dt
        self.vy += 0.5*fy*dt

Verlet Algorithm

Challenge 2.3: