c CCLRC
Section 3.3
(d) Steps (b) and (c) are repeated until all bondlengths satisfy the convergence criterion
(this iteration constitutes stage 2 of the RATTLE VV1 algorithm).
2. Forces calculated afresh.
3. RATTLE stage 2:
(a) All atom velocities are updated to a full step, assuming an absence of rigid bonds. (This
is stage 1 of the RATTLE VV2 algorithm.)
(b) The deviation of d
ij
· (v
j
- d
i
) in each bond is used to calculate the corresponding
constraint force that (retrospectively) `corrects' the bond velocities.
(c) After the correction (
) has been applied to all bonds, every bond velocity is checked
against the above condition. If the largest deviation found exceeds the desired tolerance,
the correction calculation is repeated.
(d) Steps (b) and (c) are repeated until all bonds satisfy the convergence criterion (this
iteration constitutes stage 2 of the RATTLE VV2 algorithm).
The parallel version of the RATTLE algorithm, as implemented in DL POLY 3 , is derived from
the RD SHAKE algorithm [
] although its implementation in the Domain Decomposition frame-
work requires no global merging operations and is consequently significantly more efficient. The
routine constraints shake is called to apply corrections to the atomic positions and the routine
constraints rattle to apply corrections to the atomic velocities of constrained particles.
It should be noted that the fully converged constraint forces G
ij
make a contribution to the system
virial and the stress tensor.
The contribution to be added to the atomic virial (for each constrained bond) is
W = -d
ij
· G
ij
.
(3.17)
The contribution to be added to the atomic stress tensor (for each constrained bond) is given by
= d
ij
G
ij
,
(3.18)
where and indicate the x, y, z components. The atomic stress tensor derived from the pair
forces is symmetric.
3.3
Potential of Mean Force (PMF) Constraints and the Evalua-
tion of Free Energy
A generalization of bond constraints can be made to constrain a system to some point along a
reaction coordinate. A simple example of such a reaction coordinate would be the distance between
two ions in solution. If a number of simulations are conducted with the system constrained to
different points along the reaction coordinate then the mean constraint force may be plotted as a
function of reaction coordinate and the function integrated to obtain the free energy for the overall
process [
]. The PMF constraint force, virial and contributions to the stress tensor are obtained in
a manner analagous to that for a bond constraint (see previous section). The only difference is that
the constraint is now applied between the centres of two groups which need not be atoms alone.
DL POLY 3 reports the PMF constraint virial, W
P M F
, for each simulation. Users can convert this
to the PMF constraint force from
G
P M F
=
W
P M F
d
P M F
,
(3.19)
47