Computational Physics

572 Appendix A

Verlet algorithm

The Verlet algorithm is a simple method for integrating second order differential
equations of the form

x ̈(t)=F[x(t),t]. (A.39)

It is used very often in molecular dynamics simulations because of its simplicity
and stability (see Chapter 8). It had been mentioned already at the end of the
eighteenthcentury[ 8 , 9 ]. It has a fixed time discretisation interval and it needs only
one evaluation of the ‘force’Fper step. The Verlet algorithm is easily derived by
adding the Taylor expansions for the coordinatexatt=±haboutt=0:

x(h)=x( 0 )+hx ̇( 0 )+

h^2

2

F[x( 0 ),0]+

h^3

6

...

x( 0 )+O(h^4 ) (A.40)

x(−h)=x( 0 )−hx ̇( 0 )+

h^2

2

F[x( 0 ),0]−

h^3

6

...

x( 0 )+O(h^4 ) (A.41)

leading to

x(h)= 2 x( 0 )−x(−h)+h^2 F[x( 0 ),0]+O(h^4 ). (A.42)

Knowing the values ofxat time 0 and−h, this algorithm predicts the value ofx(h).
We always need the last two values ofxto produce the next one. If we only have the
initial position,x( 0 ), and velocity,v( 0 ), at our disposal, we approximatex(h)by

x(h)=x( 0 )+v( 0 )h+

h^2

2

F[x( 0 ),0]+O(h^3 ). (A.43)

The algorithm is invariant under time reversal (as is the original differential equa-
tion). This means that if, after having integrated the equations for some time, time
is reversed by swapping the two most recent values ofx, exactly the same path in
phase space should be followed backward in time. In practice, there will always be
small differences as a result of finite numerical precision of the computer.
As we shall see in Problem A.3, the accumulated error in the position after a large
number of integration steps is of orderO(h^2 ). As there is no point in calculating
the velocities with higher precision than that of the positions, we can evaluate the
velocities using the simple formula

v( 0 )=

x(h)−x(−h)

2 h

(A.44)

which gives the velocity with the requiredO(h^2 )error.
Using the velocities at half-integer time steps:

v[(n+ 1 / 2 )h]={x[(n+ 1 )h]−x(nh)}/h+O(h^2 ), (A.45)