Exercises 253Exercises
8.1 [C] For coding the leap-frog method(Eq. (8.7))two arrays are needed, one
containing the velocities at timest=(n+ 1 / 2 )h, and one for the positions att=nh.
The same holds for the velocity-Verlet algorithm.
At first sight it might seem that the Verlet algorithm would need more memory:
arrays containing the positions at timest=nh,t=(n− 1 )handt=(n+ 1 )h.
However, the valuexi[(n− 1 )h]can be overwritten byxi[(n+ 1 )h]. Use this to code
the Verlet algorithm such that only two arrays are needed. Test it for a number of
particles moving in one dimension and subject to the harmonic oscillator potential.
8.2 The neighbour list proposed by Verlet[8]needs updating every 10–20 integration
steps and this update requires of the order ofN^2 steps for a system containingN
particles. Another bookkeeping device consists of partitioning the system into cubic
volumes and keeping track of which particles are to be found in each of these
volumes. Consider a two-dimensionalL×Lsystem for convenience. We split this
up intoP×Psquares of linear sizeL/P.Pis chosen such that the potential can be
cut off safely beyondL/P. Suppose we have for each square a list of particles within
that volume. These lists will change whenever a particle leaves a square and moves
to a neighbouring one. The force evaluation now includes only particle pairs whose
members are either in the same or in neighbouring cells.
(a) How many particles are on average to be found in one square?
(b) How many pair forces are on average taken into account in this ‘cell method’?
(c) Calculate the gain in speed with respect to the method in which all pair
interactions are taken into account, assuming that the particles are distributed
more or less homogeneously over the volume.
8.3 The first molecular dynamics simulations were carried out by Alder and Wainwright
for hard spheres[14]. The discontinuity in the potential calls for a different approach
than that used for smooth potentials. The state of the system is given by the positions
riand velocitiesvi(ilabels the particles) at some timetiwhich is usually the time of
the last collision experienced byi. We must calculate the velocity changes for the
next pair undergoing a collision.
We consider the elastic collision between two hard spheres,iandj, which are
moving with velocitiesviandvj. At timettheir positions areriandrj. After the
collision, velocities arev′iandv′jrespectively. The sphere diameter isσ.
(a) Show, using energy and momentum conservation, that the changes in velocities
of the two particles are given by
vi=vi′−vi=−vj=
rij(vij·rij)
σ^2
wherevij=vi−vjandrij=ri−rjat the collision.
For each pair of particles we need to know the time at which they will collide
(note that because of PBC each pair will indeed collide at some time unless the