Velocity Verlet Algorithm

The goal of this section is to introduce the velocity Verlet algorithm for molecular dynamics (MD) simulations and to show how to implement it to compute time-correlation functions (e.g., $\langle v(t) v(0) \rangle$) and transport coefficients (e.g., the diffusion coefficient $D$).

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

$$\langle v(0) v(t) \rangle = Tr\{\rho \hspace{.1cm} (v(0) v(t))\}.$$ (502)

The computation of $\langle v(0) v(t) \rangle$ thus requires sampling initial conditions according to the ensemble distribution described by $\rho$ and for each initial condition compute the value of the particle velocity $v(t)$ at time $t$.

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 $v(t+\tau)$ and position $x(t+\tau)$ at time t+$\tau$ as follows:

$$\begin{array}{ll}x(t+\tau) &= x(t) + v(t) \tau + \frac{1}{2}......(t+\tau) &= v(t) + \frac{f(t)+f(t+\tau)}{2 m} \tau,\end{array}$$ (503)

where $\tau$ is a small time increment, $m$ is the particle mass and $f(t)$ is the total force acting on the particle at time t. Given the initial conditions $x(0)$ and $v(0)$ one can compute $v(t)$ and $x(t)$ simply by applying Eqs.(503) successively n times, with $n=t/\tau$.

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 $t_j=j*\tau$ (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

$$\langle v(0) v(t) \rangle = \frac{1}{T} \int_0^T dt' v(t') v(t'+t),$$ (504)

where $T$ is the time of the measurement of the diffusion constant $D$, 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).