Computational Physics

8.2 Molecular dynamics at constant energy 205

we look at Header(IX,IY,IZ) and then move down from particle I to the following
by taking for the next particle the value Link(I). Using this in the force calculation
leads to the pseudocode:

FOR all cells with indices (IX,IY,IZ) DO

{Fill the list xt, yt and zt with the particles of the central cell}

icnt = 0;

j = Header(IX,IY,IZ);

WHILE (j>0) DO

j = link(j);

icnt = icnt + 1;

xt(icnt) = x(j); yt(icnt) = y(j); zt(icnt) = z(j);

LocNum = icnt;

END WHILE

{Now, LocNum is the number of particles in the central cell}

FOR half of the neighbouring cells DO

Find particles in the same way as central cell

and append them to the list xt, yt, zt;

END FOR

Calculate Lennard–Jones forces between all particles in the central cell;

Calculate Lennard–Jones forces between particles in central and

neighbouring cells;

END FOR

Note that we loop over onlyhalfthe number of neighbouring cells in order to avoid
double counting of particle pairs. The cell method is less efficient than the neighbour
list method as the blocks containing possible interaction candidates for each particle
substantially bigger than the spheres of the neighbour list. The advantage of the
present method lies in its suitability for parallel computing (see Chapter 16).
Cutting off the force violates energy conservation although the effect is small if
the cut-off radius is chosen suitably. To avoid this violation, the pair potentialU(r)
can be shifted so that it becomes continuous atrcut-off. The shifted potential can be
written in terms of the original one as

Ushift(r)=U(r)−U(rcut-off). (8.12)

The force is not affected by this shift; it remains discontinuous at the cut-off and
this gives rise to inaccuracies in the integration. Applying a shift in the force in
addition to the shift in the potential yields [ 9 , 10 ]

Uforce shift(r)=U(r)−U(rcut-off)−

d

dr

U(rcut-off)(r−rcut-off) (8.13)