Calculating time correlation functions 513can be computed rather easily within a molecular dynamics simulation (Berne and
Harp, 1970a; Berne and Harp, 1970b; Berne, 1971). In particular, we will discuss three
commonly employed approaches for generating these functions.
13.4.1 The direct method
The most straightforward and rigorous approach for computing time correlation func-
tions from molecular dynamics calculations is based on a direct interpretation of eqn.
(13.2.30). In this scheme, a set of configurations is sampled from the equilibrium dis-
tributionf 0 (H(x)). This can be achieved either by a molecular dynamics or a Monte
Carlo simulation. If molecular dynamics is used to carry out the sampling, appropriate
thermostatting and/or barostatting techniques should be employed to generate the
desired equilibrium distribution. For Monte Carlo, the Metropolis or hybrid Monte
Carlo methods of Chapter 7 (see Sections 7.3 and 7.4) can be used for a canonical
distribution; these can be supplemented, if necessary, with volumesampling to gen-
erate an isothermal-isobaric distribution. An advantage of molecular dynamics is that
well-equilibrated coordinatesandvelocities are generated and can be retained for sub-
sequent calculation of the time correlation function. Given an adequate sampling of
f 0 (H(x)), each of the configurations obtained is used as an initial condition to Hamil-
ton’s equations in order to generate a dynamical trajectory of a pre-specified length.
The time correlation function is then computed by performing the average required
by eqn. (13.2.30) over all trajectories at each time step.
Suppose we have sampledKconfigurations fromf 0 (H(x)), which we will denote
x(λ),λ= 1,...,K. Note that x stands for a complete phase space vector ofdNcoordi-
nates anddNmomenta inddimensions. If each x(λ)is used to generate a trajectory
ofMsteps, then we will have, overall,Ktrajectories consisting of configurations at
discrete time pointst= 0,∆t,2∆t,....,M∆t, denoted{x( 0 λ),x(∆λt),...,x(Mλ)∆t}, where x( 0 λ)
is the sampled configuration x(λ). Then, the autocorrelation functionCAB(t) between
two propertiesAandBwith corresponding phase space functionsa(x) andb(x) at
timet=n∆tis given by
CAB(n∆t) =1
K
∑K
λ=1a(x(λ))b(x(nλ∆)t). (13.4.1)Thus, by lettingnrun from 0 toM, the time correlation function atMtime points
will be generated. The calculational procedure is illustrated in Fig. 13.5(a). The length
of each trajectory is determined by the decay or correlation time of the particular time
correlation function under consideration.
While eqn. (13.4.1) is easy to implement, its convergence usually requires a large
number of initial configurations to be sampled fromf 0 (H(x)) and subsequent genera-
tion of a large number of trajectories. The inefficiency of this method lies in the fact
that each trajectory gives just a single contribution to each time point in the correla-
tion function. Our next approach circumvents this problem using a single trajectory.
13.4.2 Time correlation functions from a single trajectory
The use of a single trajectory to compute a time correlation function relies on two
assumptions: 1) the system under study is large enough that the thermodynamic limit