212 Molecular dynamics simulations
a steady increase or decrease, and thenoise, fluctuations on top of the drift. Drift is
obviously most undesirable. In microcanonical MD we want to sample the points in
phase space with a given energy; these points form a hypersurface in phase space –
the so-calledenergy surface. If the system drifts away steadily from this plane it is
obviously not in equilibrium.
It is very important to distinguish in all these cases between two sources of
error: those resulting from the numerical integration method as opposed to those
resulting from finite precision arithmetic, inherent to computers. For example, we
shall see below that the Verlet algorithm is not susceptible to energy drift in exact
arithmetic. Drift will however occur in practice as a result of finite precision of
computer arithmetic, and although different formulations of the Verlet algorithm
have different susceptibility to this kind of drift, this depends also on the particular
way in which numbers are rounded off in the computer.
Recently, there has been much interest insymplectic integrators. After con-
sidering the Verlet algorithm in some detail, we shall describe the concept of
symplecticity^2 and its relevance to numerical integration methods.
8.4.1 The Verlet algorithm revisited
Properties of the Verlet algorithmIn this section we treat the Verlet algorithm
x(t+h)= 2 x(t)−x(t−h)+h^2 F[x(t)] (8.26)in more detail with emphasis on issues which are relevant to MD. A derivation of
this algorithm can be found in Appendix A7.1. The error per integration step is
of the orderh^4. Note that we take the mass of the particle(s) involved equal to 1.
Unless stated otherwise, we analyse the one-dimensional single-particle version of
the algorithm. The momenta are usually determined as
p(t)=[x(t+h)−x(t−h)]/( 2 h)+O(h^2 ). (8.27)Note that there is no need for a more accurate formula, as the accumulated error in
the positions after many steps is also of orderh^2. We shall check this below, using
also a more accurate expression for the momenta [16]:
p(t)=[x(t+h)−x(t−h)]/( 2 h)−
h
12{F[x(t+h)]−F[x(t−h)]} +O(h^3 ).
(8.28)This form can be derived by subtracting the Taylor expansions forx(t+h)and
x(t−h)aboutt, and approximating dF[x(t)]/dtby{F[x(t+h)]−F[x(t−h)]}/h.
(^2) Some authors use the term ‘symplecticness’ instead of ‘symplecticity’.