1549380323-Statistical Mechanics Theory and Molecular Simulation

106 Microcanonical ensemble

Substituting in forri(∆t), we obtain a set ofNcnonlinear equations for theNcun-
known multipliersλ ̃ 1 ,..., ̃λNc:

σl

(

r′ 1 +

1

m 1

∑

k

̃λk∇ 1 σk(0),...,r′N+^1

mN

∑

k

̃λk∇Nσk(0)

)

= 0. (3.9.10)

Unless the constraints are of a particularly simple form, eqns. (3.9.10) will need to
be solved iteratively. A simple procedure for doing this, known as theSHAKE algo-
rithm (Ryckaertet al., 1977), proceeds as follows. First, if a good initial guess of the

solution,{λ ̃
(1)
k }, is available (for example, the multipliers from the previous molecular
dynamics time step), then the coordinates can be updated according to

r(1)i =r′+

1

mi

∑

k

̃λ(1)

k ∇iσk(0). (3.9.11)

The exact solution for the multipliers is now written as ̃λk= ̃λ(1)k +δ ̃λ(1)k , andri(∆t) =

r(1)i + (1/mi)

∑

kδ ̃λ

(1)

k ∇iσk(0), so that eqn. (3.9.10) becomes

σl

(

r(1) 1 +

1

m 1

∑

k

δλ ̃k∇ 1 σk(0),...,r(1)N +

1

mN

∑

k

δλ ̃(1)k ∇Nσk(0)

)

= 0. (3.9.12)

Next, eqn. (3.9.12) is expanded to first order in a Taylor series aboutδ ̃λ
(1)
k = 0:

σl(r(1) 1 ,...,r(1)N) +

∑N

i=1

∑Nc

k=1

1

mi

∇iσl(r(1) 1 ,...,r(1)N)·∇iσk(r 1 (0),...,rN(0))δλ ̃(1)k ≈ 0. (3.9.13)

Eqn. (3.9.13) is a matrix equation for the changesδ ̃λ(1)k in the multipliers. If the
dimensionality of this equation is not too large, then it can be inverteddirectly to
yield the full set ofδ ̃λ(1)k simultaneously. This procedure is known as matrix-SHAKE
or M-SHAKE (Kraeutleret al., 2001). Because eqn. (3.9.12) was approximated by

linearization, however, adding the correction

∑

kδ ̃λ

(1)

k ∇iσk(0) tor

(1)
i does not yield
a fully convergedri(∆t). We, therefore, definer(2)i =r(1)i + (1/mi)

∑

kδ
̃λ(1)
k ∇iσk(0)
and writeri(∆t) =r(2)i +(1/mi)

∑

kδ
λ ̃(2)
k ∇iσk(0) and use eqn. (3.9.13) with the “(1)”
superscript replaced by “(2)” for another iteration. The procedure is repeated until
the constraint conditions are satisfied to a given small tolerance.
If the dimensionality of eqn. (3.9.13) is high due to a large number of constraints,
then a further time-saving approximation can be made. We replace the full matrix
Alk=

∑N

i=1(1/mi)∇iσl(r

(1)

1 ,...,r

(1)

N)·∇iσk(r^1 (0),...,rN(0))δ ̃λ

(1)
k by its diagonal ele-
ments only, leading to

σl(r(1) 1 ,...,r(1)N) +

∑N

i=1

1

mi

∇iσl(r(1) 1 ,...,r(1)N)·∇iσl(r 1 (0),...,rN(0))δ ̃λ(1)l ≈ 0. (3.9.14)