216 Molecular dynamics simulations
as we shall see below. Most important for molecular dynamics is the property
that the total energy fluctuates within a narrow range around the exact one. Some
comparison has been carried out between nonsymplectic phase space conserving
and symplectic integrators [25], and this gave no indication of the superiority of
symplectic integrators above merely phase-space conserving ones. As symplectic
integrators are not more expensive to use than nonsymplectic time-reversible ones,
their use is recommended as the safest option. Investigating the merits of the various
classes of integration methods for microcanonical MD is a fruitful area for future
research.
Frictional forcesLater we shall encounter extensions of the standard MD method where a fric-
tional force is acting on the particles along the direction of the velocity. The Verlet
algorithm can be generalised to include such frictional forces and we describe this
extension for the one-dimensional case which can easily be generalised to more
dimensions. The continuum equation of motion is
̈x=F(x)−γx ̇, (8.36)and expandingx(h)andx(−h)aroundt=0 in the usual way (seeAppendix A7.1)
gives
x(h)=x( 0 )+hx ̇( 0 )+h^2 [−γx ̇( 0 )+F( 0 )]/ 2 +h^3
...x
( 0 )/ 6 +O(h^4 ) (8.37a)
x(−h)=x( 0 )−hx ̇( 0 )+h^2 [−γx ̇( 0 )+F( 0 )]/ 2 −h^3...
x( 0 )/ 6 +O(h^4 ). (8.37b)Addition then leads to
x(h)= 2 x( 0 )−x(−h)+h^2 [−γx ̇( 0 )+F( 0 )]+O(h^4 ) (8.38)wherex ̇( 0 )remains to be evaluated. If we write
x ̇( 0 )=[x(h)−x(−h)]/( 2 h)+O(h^2 ), (8.39)and substitute this into( 8. 38 ),we obtain
( 1 +γh/ 2 )x(h)= 2 x( 0 )−( 1 −γh/ 2 )x(−h)+h^2 F( 0 )+O(h^4 ). (8.40)A leap-frog version of the same algorithm is
x(h)=x( 0 )+hp(h/ 2 ); (8.41a)p(h/ 2 )=( 1 −γh/ 2 )p(−h/ 2 )+hF( 0 )
1 +γh/ 2. (8.41b)
If the massmis not equal to unity, the factors 1±γh/2 are replaced by 1±γh/( 2 m).
It is often useful to simulate the system with a prescribed temperature rather
than at constant energy. InSection 8.5we shall discuss a constant-temperature MD