236 Molecular dynamics simulations
equations of motion the norm ofnˆis not rigorously conserved; it can suffer from
numerical errors which may growing steadily with time. We shall now see how this
can be avoided.
Let us write down the leap-frog algorithm for the equation of motion (8.110)
fornˆ:
p(t+h/ 2 )=p(t−h/ 2 )+h[−ω^2 nˆ(t)+N(t)×nˆ(t)/I] (8.111a)
nˆ(t+h)=nˆ(t)+hp(t+h/ 2 ). (8.111b)Hereprepresents the time-derivative ofnˆat timest=(n+ 1 / 2 )h. The problem
with this algorithm is that it depends onω^2 and this depends in turn on the time
derivativen ̇ˆ. A convenient way of findingω^2 is to use
p(t−h/ 2 )=p(t)−
h
2(−ω^2 nˆ+N×nˆ/I)+O(h^2 ), (8.112)so that, usingnˆ(t)·p(t)=0, we obtain
2 p(t−h/ 2 )·nˆ(t)=hω^2 +O(h^2 ). (8.113)Calling the left hand side of this equationλ, we have [2, 47]
λ= 2 p(t−h/ 2 )·nˆ (8.114a)
p(t+h/ 2 )=p(t−h/ 2 )+h[nˆ(t)×N(t)/I−λnˆ(t)] (8.114b)
nˆ(t+h)=nˆ(t)+hp(t+h/ 2 ). (8.114c)The continuum equations of motion guaranteed conservation of the norm of the unit
vectornˆ. The leap-frog algorithm will enforce this normalisation only up to an error
ofh^3. It is therefore sensible to normalisenˆafter every time step – the parameter
λcan then be viewed as the Lagrange multiplier associated with the constraint
|nˆ|^2 =1. In the next section we shall discuss a simpler method for simulating
liquid nitrogen, using constraints in a different way.
For general molecules, we have an extra degree of freedom: the angle of rotation
around a molecular axis. This is the third Euler angle, which is usually denoted as
γ. The straightforward procedure for solving the equations of motion is to calculate
the principal angular velocityωωωin terms of the time derivatives of the Euler angles.
The Euler equation of motion gives the rate of change inωωωin terms of the torque. The
time derivatives of the Euler angles can then be found again fromωωω, and these can
be used to calculate the new atomic positions. There is however a problem when the
Euler angleθ=0, as in that case the transformation fromωωωto the time derivatives
of the Euler angles becomes singular. Evans has discussed this problem and has
presented methods to avoid the instability resulting from this singularity[48]. The
most efficient one is to use the quaternion representation, in which the orientation
of the molecule is defined in terms of a four-dimensional unit vector rather than
three Euler angles. This method was implemented by Evans and Murad[49].