Computational Physics

208 Molecular dynamics simulations

function beyondrcut-off:

〈U〉=〈U〉cut-off+ 2 π

N(N− 1 )

V

∫∞

rcut-off

r^2 drU(r)g(r) (8.18)

where〈···〉cut-offis the average restricted to pairs with separation smaller than
rcut-off. Of course, we can determine the correlation function forrup to half the
linear system size only because of periodic boundary conditions. Verlet[12]has
used the Percus–Yevick approximation to extrapolategbeyond this range. Ofteng
is simply approximated by its asymptotic valueg(r)≡1 for larger.
Similarly, the virial equation is corrected for the potential tail:
P
nkBT

= 1 −

1

3 NkBT

〈

∑

i

∑

j>i

rij

∂U(R)

∂rij

〉

cut-off

−

2 πN

3 kBTV

∫∞

rcut-off

r^3

∂U(r)

∂r

g(r)dr,

(8.19)
whereg(r)can also be replaced by 1.
The specific heat can be calculated from Lebowitz’s formula, seeEq. (7.37).

8.3 A molecular dynamics simulation program for argon

In the previous section we described the structure of a MD program and here we give
some further details related to the actual implementation. The program simulates
the behaviour of argon. In 1964, Rahman[13]published a paper on the properties
of liquid argon – the first MD simulation involving particles with smoothly varying
potentials. Previous work by Alder and Wainwright[14]was on hard sphere fluids.
Rahman’s work was later refined and extended by Verlet[8]who introduced several
features that are still used, as we have seen in the previous section.
The Lennard–Jones pair potential turns out to give excellent results for argon:

U(r)= 4 ε

[(

σ

r

) 12

−

(σ

r

) 6 ]

. (8.20)

The optimal values for the parametersεandσareε/kB=119.8 K andσ =
3.405 Å respectively.
In the initialisation routine, the positions of a face centred cubic lattice are gen-
erated. For anL×L×Lsystem containing 4M^3 particles, the fcc lattice constanta
isa=L/M. It may be safer to put the particles not exactly on the boundary facets
of the system because as a result of rounding errors it might not always be clear
whether they belong to the system under consideration or a neighbouring copy.
The procedure inAppendix B3for generating random numbers with a Gaussian
distribution should be used in order to generate momenta according to a Maxwell
distribution. First generate all the momenta with some arbitrary distribution width.
Then calculate the total momentumptotand subtract a momentump ̄=ptot/Nfrom