Computational Physics

8.2 Molecular dynamics at constant energy 207
Continue simulation and determine physical quantities:Integration of the
equations of motion proceeds as described above. In this part of the simulation,
the actual determination of the static and dynamic physical quantities takes place.
We determine the expectation value of a static physical quantity as a time average
according to

A=

1

n−n 0

∑n

ν>n 0

Aν. (8.16)

The indicesνlabel thentime steps of the numerical integration, and the firstn 0
steps have been carried out during the equilibration. For determination of errors in
the measured physical quantities, see the discussion in Section 7.4.
Difficulties in the determination of physical quantities may arise when the para-
meters are such that the system is close to a first or second order phase transition
(see the previous chapter): in the first order case, the system might be ‘trapped’ in
a metastable state and in the second order case, the correlation time might diverge
for large system sizes.
In the previous chapter we have already considered some of the quantities of
interest. In the case of a microcanonical simulation, we are usually interested in the
temperature and pressure. Determination of these quantities enables us to determine
theequation of state, a relation between pressure and temperature, and the system
parameters – particle number, volume and energy (NVE). This relation is hard to
establish analytically, although various approximate analytical techniques for this
purpose exist: cluster expansions, Percus–Yevick approximation, etc. [11].
The pair correlation function is useful not only for studying the details of the
system but also to obtain accurate values for the macrosopic quantities such as the
potential energy and pressure, as we shall see below. The correlation function is
determined by keeping a histogram which contains for every interval[ir,(i+ 1 )
r]the number of pairsn(r)with separation within that range. The list can be
updated when the pair list for the force evaluation is updated. The correlation
function is found in terms ofn(r)as

g(r)=

2 V

N(N− 1 )

[

〈n(r)〉

4 πr^2 r

]

. (8.17)

Similar expressions can be found for time-dependent correlation functions – see
Refs.[2]and[11].
If the force has been cut off during the simulation, the calculation of average
values involving the potentialUrequires some care. Consider for example the
potential energy itself. This is calculated at each step taking only the pairs with
separation within the minimum cut-off distance into account; taking all pairs into
account would imply losing the efficiency gained by cutting off the potential. The
neglect of the tail of the potential can be corrected for by using the pair correlation