214 Molecular dynamics simulations
(8.27)and(8.28)does indeed scale ash^2. In the following we determine momenta
according toEq. (8.27). In the leap-frog version the momentum estimator is
p(t)=[p(t+h/ 2 )+p(t−h/ 2 )]/ 2 +O(h^2 ). (8.31)
The results for the various energy estimators can be obtained by solving the
harmonic oscillator in the Verlet algorithm analytically. The ‘Verlet harmonic
oscillator’ reads
x(t+h)= 2 x(t)−x(t−h)−h^2 x(t). (8.32)
If we substitutex(t)=exp(iωt)into the last equation, we obtain
cos(ωh)= 1 −h^2 / 2 (8.33)and this defines a frequencyωdiffering by an amount of orderh^2 from the angular
frequencyω=1 of the exact solution. The difference between the numerical and
the exact solution will therefore show a slow beat.
A striking property of the energy determined from the Verlet/leap-frog solution
is that it does not show any drift in the total energy (in exact arithmetic). This stabil-
ity follows directly from the fact that the Verlet algorithm is time-reversible, which
excludes steady increase or decrease of the energy for periodic motion. In a molecu-
lar dynamics simulation, however, the integration time, which is the duration of the
simulation, is much smaller than the period of the system, which is thePoincaré
time, that is the time after which the system returns to its starting configuration.
The error in the energy might therefore grow steadily during the simulation. It turns
out, however, that the deviation of the energy remains bounded in this case also, as
the Verlet algorithm possesses an additional symmetry, calledsymplecticity. Sym-
plecticity will be described in detail in Section 8.4.2. Here we briefly describe what
the consequences of symplecticity are for an integration algorithm. Symplecticity
gives rise to conserved quantities, and in particular, it can be shown that a discrete
analogue of the total energy is rigorously conserved (in exact arithmetic) [17]. It
turns out that this discrete energy deviates from the continuum energy at most an
amount of orderhk, for some positive integerk. Therefore, the energy cannot drift
away arbitrarily and it follows that the noise remains bounded.
To illustrate this point we return to the harmonic oscillator. In this particular
case we can actually determine the conserved discrete energy. In the leap-frog
formulation:
p(t+h/ 2 )=p(t−h/ 2 )−hx(t); (8.34a)
x(t+h)=x(t)+hp(t+h/ 2 ), (8.34b)it is equal to[18]
HD=^12 [p(t−h/ 2 )^2 +x(t)^2 −hp(t−h/ 2 )x(t)]. (8.35)