c CCLRC
Section 3.2
The RATTLE VV2, is also a two stage algorithm. In the first stage, the VV algorithm (VV2)
calculates the velocities of the atoms in the system assuming a complete absence of the rigid bond
forces (since forces have just been recalculated afresh after VV1. The relative velocity of atom i to
j (or vice versa) constituting the rigid bond ij is not perpendicular to the bond (i.e. does not have
a non-zero component along the bond) as required by the definition of rigidity. In the second stage
the deviation from zero of the scalar product d
ij
· (v
j
- v
i
) is used retrospectively to compute the
constraint force needed to keep the bond rigid over the length of the timestep t. It is relatively
simple to show that the constraint force has the form:
B
ij
µ
ij
t
d
ij
· (v
j
- v
i
)
d
2
ij
d
ij
.
(3.15)
The velocity corrections can therefore be written as
v
corr
i
= t
B
ij
m
i
=
µ
ij
m
i
d
ij
· (v
j
- v
i
)
d
2
ij
d
ij
.
(3.16)
For a system of simple diatomic molecules, computation of the constraint force will, in principle,
allow the correct atomic positions to be calculated in one pass. However, in the general polyatomic
case this correction is merely an interim adjustment, not only because the above formula is ap-
proximate, but the successive correction of other bonds in a molecule has the effect of perturbing
previously corrected bonds. Either part of the RATTLE algorithm is therefore iterative, with the
correction cycle being repeated for all bonds until: each has converged to the correct length, within
a given tolerance for RATTLE VV1 (SHAKE) and the relative bond velocities are perpendicular
to their respective bonds within a given tolerance for RATTLE VV2 (RATTLE). The tolerance
may be of the order 10
-4
°
A to 10
-8
°
A depending on the precision desired.
The SHAKE procedure may be summarised as follows:
1. All atoms in the system are moved using the LFV algorithm, assuming an absence of rigid
bonds (constraint forces). (This is stage 1 of the SHAKE algorithm.)
2. The deviation in each bondlength is used to calculate the corresponding constraint force
) that (retrospectively) `corrects' the bond length.
3. After the correction (
) has been applied to all bonds, every bondlength is checked. If the
largest deviation found exceeds the desired tolerance, the correction calculation is repeated.
4. Steps 2 and 3 are repeated until all bondlengths satisfy the convergence criterion (this iteration
constitutes stage 2 of the SHAKE algorithm).
The RATTLE procedures may be summarised as follows:
1. RATTLE stage 1:
(a) All atoms in the system are moved using the VV algorithm, assuming an absence of
rigid bonds (constraint forces). (This is stage 1 of the RATTLE VV1 algorithm.)
(b) The deviation in each bondlength is used to calculate the corresponding constraint force
) that (retrospectively) `corrects' the bond length.
(c) After the correction (
) has been applied to all bonds, every bondlength is checked.
If the largest deviation found exceeds the desired tolerance, the correction calculation is
repeated.
46