1549380323-Statistical Mechanics Theory and Molecular Simulation

Constraints 105

3.9.1 The SHAKE and RATTLE algorithms

For time-independent holonomic constraints, the Euler–Lagrangeequations of motion,
eqns. (1.9.11) and (1.9.12), in Cartesian coordinates are

d

dt

(

∂L

∂r ̇i

)

−

∂L

∂ri

=

∑Nc

k=1

λkaki

∑N

i=1

aki·r ̇i= 0, (3.9.3)

where
aki=∇iσk(r 1 ,...,rN). (3.9.4)

Note that these are equivalent to

mi ̈ri=Fi+

∑Nc

k=1

λk∇iσk

d

dt

σk(r 1 ,...,rN) = 0. (3.9.5)

The constraint problem amounts to integrating eqn. (3.9.3) subject to the conditions
thatσk(r 1 ,...,rN) = 0 and ̇σk(r 1 ,...,rN,r ̇ 1 ,...,r ̇N) =

∑

i∇iσk(r^1 ,...,rN)·r ̇i= 0. We
wish to develop a numerical scheme in which the constraint conditionsare satisfied
exactly as part of the integration algorithm.
Starting from the velocity Verlet approach, for example, we begin with the position
update, which, when holonomic constraints are imposed, reads

ri(∆t) =ri(0) + ∆tvi(0) +

∆t^2

2 mi

Fi(0) +

∆t^2

2 mi

∑

k

λk∇iσk(0), (3.9.6)

whereσk(0)≡σk(r 1 (0),...,rN(0)). In order to ensure that the constraint is satisfied
exactly at time ∆t, we impose the constraint condition directly on the numerically
obtained positionsri(∆t) and determine, on the fly, the multipliers needed to enforce
the constraint. Let us define

r′i=ri(0) + ∆tvi(0) +

∆t^2

2 mi

Fi(0) (3.9.7)

so that
ri(∆t) =r′i+

∑

k

̃λk∇iσk(0), (3.9.8)

whereλ ̃k = (∆t^2 /2)λk. Then, for each constraint conditionσl(r 1 ,...,rN) = 0, we
impose
σl(r 1 (∆t),...,rN(∆t)) = 0, l= 1,...,Nc. (3.9.9)