Franz J. Vesely > CompPhys Tutorial > Selected Applications > Evaluation of Simulation Experiments

 < >
 Part III: Ch. 6

## 6.4 Evaluation of Simulation Experiments

In addition to basic thermodynamic observables such as internal energy (eq. 6.1) and pressure (equ. 6.2) many more details about the structure and the dynamics of statistical-mechanical systems may be accessed. Let us consider two such quantities, the pair correlation function $g(r)$ and the velocity autocorrelation function $C(t)$.

In your Lennard-Jones MD and MC programs, include a procedure to calculate averages of the total potential energy and the virial as defined by

$W\equiv\sum \limits_{i}\vec{K}_{i}\cdot \vec{r}_{i} = -\frac{\textstyle 1}{\textstyle 2} \sum \limits_{i} \sum \limits_{j} \vec{K}_{ij}\cdot \vec{r}_{ij}$

From these compute the internal energy and the pressure. Compare with results from literature, e.g. [ MCDONALD 74, VERLET 67]. Allow for deviations in the $5-10\%$ range, as we have omitted a correction for the finite sample size (`cutoff correction').

Subsections

### 6.4.1 Pair Correlation Function

Let

$a(r;\Gamma_{c}) = \sum \limits_{i}\delta(r_{i}- r)$

An average of this quantity represents the relative frequency, or probability densitiy of some particle being situated near $r$. In other words, $\langle a(r)\rangle=\rho(r)$ is simply the mean fluid density at position $r$:

$\rho(r)=\langle \sum \limits_{i}\delta(r_{i}-r)\rangle$

In a fluid we usually have $\rho(r)=const$; only in the presence of external fields or near surfaces $\rho(r)$ varies in a non-trivial manner.

Let us proceed to the "pair correlation function" (PCF)

$g(r)=\frac{\textstyle V}{\textstyle 4\pi r^{2}N^{2}} \langle\sum \limits_{i} \sum \limits_{j\neq i} \delta (r-r_{ij}) \rangle$

This is the conditional probability density of finding a particle at $r$, given that there is a particle at the coordinate origin. Thus $g(r)$ provides a measure of local spatial ordering in a fluid.

To determine $g(r)$, proceed like this:
• Divide the range of $r$-values (at most $[0;L/2]$, where $L$ is the side length of the basic simulation cell) into $50-200$ intervals of length $\Delta r$.
• A given configuration $\{r_{1},\dots r_{N}\}$ is scanned to determine, for each pair $(i,j)$, a "distance channel" number

$k = \mbox{int} \left( \frac{\textstyle r_{ij}}{\textstyle \Delta r}\right)$

• In a histogram table $g(k)$ the corresponding value is then incremented by $1$. This procedure is repeated every, say, $50$ MD steps (or $50N$ MC steps).
• At the end of the simulation run the histogram is normalized according to 6.9.
The typical shape of the PCF at liquid densities is depicted in Fig. 6.6.

Figure 6.6: Pair correlation function of the Lennard-Jones liquid

Significance of $g(r)$ in fluid physics:
• The average of any quantity that depends on the pair potential $u(r)$ may be expressed as an integral over $g(r)$.

$p=\frac{\textstyle NkT}{\textstyle V}-\frac{\textstyle N^{2}}{\textstyle 6V^{2}} \int\limits_{V}r\frac{\textstyle du}{\textstyle dr} g(r) dr$

• Theory: analytical approximations to $g(r)$ for a given pair potential $u(r)$.[KOHLER 72], [HANSEN 86]
• Experiment: the Fourier transform of $g(r)$, the "scattering law"

$S(k)+ \frac{\textstyle N}{\textstyle V} \int \limits_{V} [g(r)-1] e^{\textstyle i k\cdot r} dr$

is just the relative intensity of neutron or X-ray scattering at a scattering angle $\theta$ with

$k \equiv \frac{\textstyle 4 \pi}{\textstyle \lambda} \sin \frac{\textstyle \theta}{\textstyle 2}$

Augment your Lennard-Jones MD (or MC) program by a routine that computes the pair correlation function $g(r)$ according to the procedure given above; remember to apply the nearest image convention when computing the pair distances. As the subroutine ENERGY already contains a loop over all particle pairs $(i,j)$, it is best to increment the $g(r)$ histogram within that loop.

Plot the PCF and see whether it resembles the one given in Figure 6.6.

### 6.4.2 Autocorrelation Functions

An elementary example of temporal correlations of the form

$C_{a}(t) \equiv \langle a(0) a(t) \rangle$

is the velocity autocorrelation in fluids

$C(t)\equiv \langle v_{i}(0) \cdot v_{i}(t) \rangle$

Simple kinetic theory, which takes into account only binary collisions, predicts $C(t) \propto e^{-\lambda t}$. At fluid densities a different behavior is to be expected. Nevertheless the first results on $C(t)$ obtained by Alder [
ALDER 67] provided some surprises. It turned out that at intermediate fluid densities and long times $C(t) \propto t^{-3/2}$ instead of showing an exponential decay.

This has profound consequences. The diffusion constant $D$ of a liquid is given by

$D = \frac{\textstyle 1}{\textstyle 3} \int \limits_{0}^{\infty}C(t) dt$

Due to the long time tail of $C(t)$ the simulation result for $D$ is about $30$ percent larger than its kinetic estimate.

The reason for the long time tail in $C(t)$ was later explained as a collective dynamical effect: part of the momentum of a particle is stored in a microscopic vortex that dies off very slowly. [DORFMAN 72]

Figure 6.7: Velocity autocorrelation function of the Lennard-Jones fluid

Procedure for calculating autocorrelation functions $\langle a(0) a(t) \rangle$:

• At regular intervals of $20-100$ time steps, mark starting values $\{a(t_{0,m}),m,\dots M\}$. Since in the further process only the preceding $M \approx 10-20$ starting values are required, it is best to store them in registers that are cyclically overwritten.
• At each time $t_{n}$, compute the $M$ products

$z_{m} = a(t_{n}) \cdot a(t_{0,m}) , \;\;\; m, \dots M$

and relate them to the (discrete) time displacements $\Delta t_{m} \equiv t_{n}-t_{0,m}$; a particular $\Delta t_{m}$ defines a channel number

$k = \Delta t_{m} / \Delta t$

indicating the particular histogram channel to be incremented by $z_{m}$. To simplify the final normalization it is recommended to count the number of times each channel $k$ is incremented.

• Run your MD program for $2000$ time steps and store the velocity vector of a certain particle (say, no. 1) at each time step. Write and test a program that evaluates the autocorrelation function of this vector.