1549380323-Statistical Mechanics Theory and Molecular Simulation

100 Microcanonical ensemble

ri(t+ ∆t) +ri(t−∆t) = 2ri(t) +

∆t^2

mi

Fi(t), (3.8.4)

which, after rearrangement, becomes

ri(t+ ∆t) = 2ri(t)−ri(t−∆t) +

∆t^2

mi

Fi(t). (3.8.5)

Eqn. (3.8.5) is a numerical solver known as theVerlet algorithm(Verlet, 1967). Given a
set of initial coordinates{r 1 (0),...,rN(0)}and initial velocities{v 1 (0),...,vN(0)}, eqn.
(3.8.2) can be used to obtain a set of coordinates,{r 1 (∆t),...,rN(∆t)}, after which
eqn. (3.8.5) can be used to generate a trajectory of arbitrary length. Note that the
Verlet algorithm only generates positions. If needed, the velocitiescan be constructed
at any point in the trajectory via the centered difference formula

vi(t) =

ri(t+ ∆t)−ri(t−∆t)

2∆t

. (3.8.6)

3.8.2 The velocity Verlet algorithm

Although appealing in its simplicity, the Verlet algorithm does not explicitly evolve
the velocities, and this is somewhat inelegant, as phase space is composed of both
positions and velocities (or momenta). Here, we will derive a variant of the Verlet
integrator, known as thevelocity Verletalgorithm (Swopeet al., 1982), that explicitly
evolves positions and velocities. Consider, again, the expansion of the coordinates up
to second order in ∆t:

ri(t+ ∆t)≈ri(t) + ∆tvi(t) +

∆t^2

2 mi

Fi(t). (3.8.7)

Interestingly, we could alsostartfromri(t+ ∆t) andvi(t+ ∆t), computeFi(t+ ∆t)
and evolve backwards in time tori(t) according to

ri(t) =ri(t+ ∆t)−∆tvi(t+ ∆t) +

∆t^2

2 mi

Fi(t+ ∆t). (3.8.8)

Substituting eqn. (3.8.7) forri(t+ ∆t) into eqn. (3.8.8) and solving forvi(t+ ∆t)
yields

vi(t+ ∆t) =vi(t) +

∆t

2 mi

[Fi(t) +Fi(t+ ∆t)]. (3.8.9)

Thus, the velocity Verlet algorithm uses both eqns. (3.8.7) and (3.8.9) to evolve the
positions and velocities simultaneously. The Verlet and velocity Verletalgorithms sat-
isfy two properties that are crucial for the long-time stability of numerical solvers.
The first is time-reversibility, which means that if we take as initial conditionsr 1 (t+
∆t),...,rN(t+ ∆t),v 1 (t+ ∆t),...,vN(t+ ∆t) and step backward in time using a time
step−∆t, we will arrive at the stater 1 (t),...,rN(t),v 1 (t),...,vN(t). Time-reversibility