8.7 Long-range interactions 241Each step will therefore shift the positions more closely to the point where they all
satisfy the constraint. The iterative process is stopped when all the constraints are
smaller (in absolute value) than some small positive number.
The algorithm can be summarised as follows:
CalculateR ̃(t+h)using (8.125);
SetRoldequal toR ̃(t+h);
REPEAT
Calculateλ(kl)from(8.129);
FORk=1TOMDO
SetRoldequal toRnew
UpdateRoldtoRnewusing(8.127);
END FOR
UNTIL Constraints are satisfied.
The SHAKE algorithm turns out to be quite efficient: for a system of 48 atoms with
112 constraints, typically 25 iterations are necessary in order to achieve convergence
of the constraints within a relative accuracy of 5× 10 −^7 [52].
8.7 Long-range interactions
Coulombic and gravitational many-particle systems are of great interest because
they describe plasmas, electrolytic solutions, and celestial mechanical systems. The
interaction is described by a pair-potential which in three dimensions is proportional
to 1/r– in two dimensions it is lnr. The long range character of this potential poses
problems. First of all, it is not clear whether the potential can be cut off beyond some
finite range. One might hope that for a charge-neutral Coulomb system, screening
effects could justify this procedure. Unfortunately, for most systems of interest, the
screening length exceeds half the linear system size that can be achieved in practice,
so we cannot rely on this screening effect to justify cutting off the potential, as this
would essentially alter the form of the screening charge cloud. Also, when using
the minimum image convention with periodic boundary conditions, equally charged
particles tend to occupy opposite ends of a half diagonal of the system unit cell in
order to minimise their interaction energy, thus introducing unphysical anisotropies.
Therefore, we cannot cut off the potential and all pairs of interacting particles must
be taken into account when calculating the forces.
Connected with this is an essential difference in the treatment of periodic or
nonperiodic system cells. In the latter case, we simply use the 1/rpotential (or
lnrin two dimensions), but in the periodic case we must face the problem that in
general the sum over the image charges in the periodic replicas does not converge.
This can be remedied by subtracting an offset from the potential (note that adding