Computational Physics

240 Molecular dynamics simulations

by chemical bonds can be described in terms of stretching, bending and torsion, as
described in Section 8.6.2. In general, partial constraints cannot be treated using
the methods given previously. Trying to use the constraints to reduce the equa-
tions to a smaller set and formulating equations for the rigid fragments in terms of
quaternions is quite complicated. Ryckaertet al.[51–54] devised a simple and effi-
cient iterative method for treating arbitrary constraints which is now still the most
important method for MD with polyatomic molecules. Analogous to the method
of constraints for rigid molecules, the rigidity of the fragments can be imposed by
constraints, which are all expressed in Cartesian coordinates for simplicity. The
Lagrange multipliers are determined after each integration step by substituting the
new positions into the constraint equations.
The algorithm, called SHAKE, is formulated within the framework of the Verlet
algorithm. The forces experienced by the particles consist of physical and of con-
straint forces. The constraints are given byσk(R)=0, wherek=1,...,M;M
is the number of constraints. We denote the physical force on particleibyFiand
the constraint force is

∑M

k= 1 λk∇iσk(R), whereλkis the Lagrange multiplier which
is to be determined. At time stept=nhwe have at our disposal the positions at
timestandt−h. These positions satisfy the constraint equationsσk(R)=0to
numerical precision. The aim is to find the positions at timet+h, satisfying the
constraint equation. First we calculate the new positionsr ̃i(t+h)without taking
the constraints into account:

̃ri(t+h)= 2 ri(t)−ri(t−h)+h^2 Fi[ri(t)]. (8.125)

The final positionsri(t+h)can be written as

ri(t+h)=r ̃i(t+h)−

∑M

k= 1

λk∇iσk[R(t)]. (8.126)

Theλkare found in an iterative procedure. We number the iterations by an index
l. In each iteration, a loop over the constraintskis carried out, and in each step of
this loop, the Lagrange parameterλkand all the particle positions are updated. The
positions are updated according to

rinew=riold−h^2 λ(kl)∇iσk[R(t)]. (8.127)

The parameterλ(kl)is found from a first order expansion ofσk(R), requiring that
this vanishes:

σk[Rnew]≈σk[Rold]−h^2 λ(kl)

∑

i

∇iσk[Rold]∇iσk[R(t)]=0, (8.128)

leading to

λ(kl)=

σk[Rold]

h^2

{∑

i∇iσk[R

old]∇iσk[R(t)]}. (8.129)