8.2 Molecular dynamics at constant energy 203This form is most convenient because it is very stable with respect to errors due to
finite precision arithmetic, and it does not require additional calculations in order
to find the velocities. It should be noted that all formulations have essentially the
same memory requirements. It may seem that, as this algorithm needstwoforces
the second step, we need two arrays for these, one containingF(t)and the other
F(t+h). However, the following form of the algorithm is exactly equivalent and
avoids the need for two force arrays:
v ̃(t)=v(t)+hF(t)/2, (8.10a)
r(t+h)=r(t)+hv ̃(t), (8.10b)
v(t+h)=v ̃(t)+hF(t+h)/2. (8.10c)The new forceF(t+h)is calculated between the second and third step.
The force acting on particleiresults from the interaction forces between this
particle and all the other particles in the system – usually pair-wise interactions
are used. The calculation of the forces therefore takes a relatively long time as
this requiresO(N^2 )steps. A problem in the evaluation of the force arises from
the assumption of periodic boundary conditions. These imply that the system is
surrounded by an infinite number of copies with exactly the same configuration
as in Figure 8.1. A particle therefore interacts not only with each partnerjin the
system cell we are considering but also with the images of particlejin all the copies
of the system. This means that in principle an infinite number of interactions has
to be summed over. In many cases, the force decays rapidly with distance, and
in that case remote particle copies will not contribute significantly to the force. If
the force between the particles can safely be neglected beyond separations of half
the linear system size, the force evaluation can be carried out efficiently by taking
into account, for each particle in the system, only the interactions with the nearest
copy of each of the remaining particles (seeFigure 8.1): each infinite sum over all
the copies is replaced by a single term! This is theminimum image convention.In
formula, for a cubic system cell the minimum image convention reads
rijmin=min
n
|ri−rj+nμLμ| (8.11)with the same notation as inEq. (8.3), but where the components ofnμassume
the values 0,±1, provided all the particles are kept within the system cell, by
translating them back if they leave this cell. The potential is no longer analytic in
this convention, but discontinuities will obviously be unimportant if the potential
is small beyond half the linear system size.
Often it is possible to cut the interactions off at a distancercut-offsmaller than half
the linear system size without introducing significant errors. In that case the forces
do not have to be calculated for all pairs. However, all pairs must be considered to