Consider the task of computing the diffusion coefficient according to Eq. (496). The quantity of interest is, therefore, the equilibrium ensemble average of the velocity autocorrelation function

(502) |

The computation of thus requires sampling initial conditions according to the ensemble distribution described by and for each initial condition compute the value of the particle velocity at time .

The *velocity Verlet algorithm* (W. C. Swope, H. C. Andersen, P.
H. Berens, and K. R. Wilson, *J. Chem. Phys.* **76**, 637 (1982))
computes the particle velocity
and position
at time t+
as follows:

(503) |

where is a small time increment, is the particle mass and is the total force acting on the particle at time t. Given the initial conditions and one can compute and simply by applying Eqs.(503) successively n times, with .

Note that by implementing the algorithm one generates a *phase space
trajectory* (i.e., a sequence of ``snapshots'' for the particle coordinates
and velocities at all intermediate times
(with j=1,2,...,n)).

Molecular dynamics simulations thus provide the sequence of microscopic configurations through which the model system passes in time. Such detailed microscopic information allows one to compute the result of a measurement of an observable (i.e., an ensemble average) according to the time average introduced by Eq. (111) (i.e., simply by averaging the value of the observable throught the whole manifold of microscopic configurations generated during the time of the measurement). Therefore, another way of computing the ensemble average introduced by Eq. (502) is

(504) |

where is the time of the measurement of the diffusion constant , a time that is much larger than the relaxation time of the velocity autocorrelation function.

**Exercise:**

Compute the velocity autocorrelation function for a fluid of argon atoms
using the program developed for computing the radial distribution function.
(Hint: substitute the Metropolis Monte Carlo algorithm by the velocity
Verlet algorithm).