1549380323-Statistical Mechanics Theory and Molecular Simulation

Constraints 107

Eqn. (3.9.14) has a simple solution

δλ ̃(1)l =−

σ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))

. (3.9.15)

Eqn. (3.9.15) could be used, for example, to obtainδ ̃λ(1) 1 followed immediately by an

update of allr(1)i involved in thek= 1 constraint to obtain positionsr(2)i for this con-

straint. Given the updated position, eqn. (3.9.15) is used to obtainδλ ̃(2) 1 immediately

followed by an update of allr
(2)
i involved in thek= 1 constraint, and we iterate until
thek= 1 constraint is satisfied. We then proceed to thek= 2 constraint and iterate
until it is satisfied, which will cause a slight violation of thek= 1 constraint. Note
that satisfying thekth constraint via this procedure causes alll < kconstraints to
be slightly violated. Thus, after cycling through all of the constraints in this manner,
the procedure must be repeated until the full set of constraintsis converged to within
a given tolerance. Once the converged multipliers are obtained, andthe coordinates
fully updated, the velocities must be updated as well according to

vi(∆t/2) =vi(0) +

∆t

2 mi

Fi(0) +

1

mi∆t

∑

k

λ ̃k∇iσk(0). (3.9.16)

After eqn. (3.9.16) is applied, we can proceed to the next step of updating the
velocities, which requires that the condition ̇σk(r 1 ,...,rN) = 0 be satisfied. Once the
new forces are obtained from the updated positions, the final velocities are written as

vi(∆t) =vi(∆t/2) +

∆t

2 mi

Fi(∆t) +

∆t

2 mi

∑

k

μk∇iσk(∆t)

=v′i+

1

mi

∑

k

μ ̃k∇iσk(∆t), (3.9.17)

whereμkhas been used to denote the multipliers for the velocity step to indicate
that they are different from those used for the position step, and ̃μk= (∆t/2)μk. The
multipliersμkare now obtained by enforcing the condition

∑N

i=1

∇iσk(∆t)·vi(∆t) = 0 (3.9.18)

on the velocities. Substituting in forvi(∆t), we obtain a set ofNclinear equations

∑N

i=1

∇iσk(∆t)·

(

v′i+

1

mi

∑

l

μ ̃l∇iσl(∆t)

)

= 0 (3.9.19)

for the multipliers ̃μl. These can be solved straightforwardly by matrix inversion or,
for large systems, iteratively by satisfying the condition for each constraint in turn and