202 Molecular dynamics simulations
wherer(t)is the position of the particle at timet=nh(his the time step;nis an
integer). From now on we choose units such thatm=1. The error per time step is of
orderh^4 and a worst case estimate for the error over a fixed time interval containing
many time steps is of orderh^2 (see Problem A3). To start up the algorithm we need
the positions of the particles at two subsequent time steps. As we have only the
initial (t=0) positions and velocity at our disposal, the positions att=hare
calculated as
r(h)=r( 0 )+hv( 0 )+h^2
2
F[r(t= 0 )] (m≡ 1 ), (8.5)with an error of orderh^3.
During the integration, the velocities can be calculated as
v(t)=r(t+h)−r(t−h)
2 h+O(h^2 ). (8.6)When using periodic boundary conditions in the simulation, we must check for each
particle whether it has left the simulation cell in the last integration step. If this is
the case, the particle is translated back over a lattice vectorLμto keep it inside
the cell (we shall see below that this procedure facilitates the common procedure
for evaluating the forces with periodic boundary conditions). The velocity must
obviously be determined before such a translation.
There exist two alternative formulations of the Verlet algorithm, which are exactly
equivalent to it in exact arithmetic but which are less susceptible to errors resulting
from finite numerical precision in the computer than the original version. The first
of these, theleap-frogform, introduces the velocities at time steps precisely in
between those at which the positions are evaluated:
v(t+h/ 2 )=v(t−h/ 2 )+hF[r(t)], (8.7a)
r(t+h)=r(t)+hv(t+h/ 2 ). (8.7b)These steps are then repeated over and over. Note that they must always be applied
in the given order: the second step usesv(t+h/ 2 )which is calculated in the first
step.
Another form is the so-called velocity-Verlet algorithm[7]which is also more
stable than the original Verlet form and which, via the definition
v(t)=r(t+h)−r(t−h)
2 h(8.8)
evaluates velocities and positions at the same time instances:
r(t+h)=r(t)+hv(t)+h^2 F(t)/2, (8.9a)
v(t+h)=v(t)+h[F(t+h)+F(t)]/2. (8.9b)