Radial distribution function: the structure of a Lennard-Jones liquid

The problem: how ordered is a liquid?

Liquid argon at 90 K, 6 K above its triple point at 83.8 K, holds 1.38 g/cm³, about 21 atoms per cubic nanometer. It has no lattice, and yet its atoms are not placed at random: around any one atom the others crowd at about one atomic diameter, thin out, and crowd again at two. How many neighbors does an atom have, and how far away do they sit?

The radial distribution function answers both. Count the atoms in a thin shell of radius \(r\) and thickness \(dr\) around an atom, average over atoms and time, and divide by the count that the mean density \(\rho\) would put there if the atoms ignored each other:

Where \(g = 1\) there is no structure. X-ray and neutron scattering measure \(g(r)\) through the structure factor \(S(q)\) the diffractometer records: \(S(q) - 1\) is the density \(\rho\) times the Fourier transform of \(g(r) - 1\).

The atoms interact through the Lennard-Jones potential \(V(r) = 4\varepsilon[(\sigma/r)^{12} - (\sigma/r)^6]\), the pair potential of the seven-atom cluster, with \(\sigma\) the distance where it crosses zero and \(\varepsilon\) the depth of its well. Lengths below are in \(\sigma\), energies in \(\varepsilon\), temperatures in \(\varepsilon/k_B\), and times in \(\tau = \sigma\sqrt{m/\varepsilon}\), a unit of time and not the \(r/v\) at closest approach of velocity Verlet. For argon \(\sigma = 0.3405\) nm, \(\varepsilon/k_B = 119.8\) K and \(\tau = 2.16\) ps, so 90 K is temperature 0.75 and the real liquid's density is 0.82 \(\sigma^{-3}\). The simulation below runs at 0.84, the classic state point of this liquid, 1.41 g/cm³ in argon units. Here 256 atoms move in a periodic cube, pushed by the velocity Verlet loop.

The first peak sits at 1.09 σ, 0.371 nm for argon, where the density is three times the mean. In the lower panel the total energy stays flat while the potential energy swings, the check that the run is the motion Newton's equations give. Step 6 draws it.

Setup

Everything is in reduced units, with \(k_B = m = 1\). The cell prints the box side and translates the state point back to argon.

import numpy as np
import matplotlib.pyplot as plt

N, RHO, T_TARGET = 256, 0.84, 0.75   # atoms, density / σ⁻³, temperature / (ε/k_B)
RC, DT = 2.5, 0.005                  # cutoff / σ, time step / τ
L = (N / RHO) ** (1 / 3)             # box side / σ
rng = np.random.default_rng(1)

SIGMA, EPS_K, MASS = 0.3405e-9, 119.8, 39.948 * 1.66054e-27   # argon: m, K, kg
tau = SIGMA * np.sqrt(MASS / (EPS_K * 1.380649e-23))

plt.rcParams.update({
    "figure.figsize": (7, 3.6), "figure.dpi": 110,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25,
    "font.size": 11, "lines.linewidth": 1.8,
})
INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"

print(f"box side L = {L:.3f} σ")
print(f"argon: T = {T_TARGET * EPS_K:.0f} K, density {RHO * MASS / SIGMA**3 / 1000:.2f} g/cm³, "
      f"τ = {tau * 1e12:.2f} ps")

box side L = 6.730 σ
argon: T = 90 K, density 1.41 g/cm³, τ = 2.16 ps

Step 1: Fill a periodic box with an fcc lattice and thermal velocities

Random starting positions put some pairs at a fraction of \(\sigma\), where the \(r^{-12}\) wall throws them apart at enormous speed. A crystal at the right density has no such overlaps and melts within a few \(\tau\). The face-centered cubic cell of general chemistry is four atoms, at \((0,0,0)\), \((\tfrac12,\tfrac12,0)\), \((\tfrac12,0,\tfrac12)\) and \((0,\tfrac12,\tfrac12)\) in units of the cell side; broadcasting adds them to the origins of \(n \times n \times n\) cells of side \(L/n\). So \(N\) must be \(4n^3\): 32, 108, 256, 500, 864, and so on, and the 256 atoms here fill 4 × 4 × 4 cells. Any other \(N\) stops the cell with an error instead of building a lattice of the wrong size.

The box is periodic: it is surrounded on every side by copies of itself, so the 256 atoms have no walls and no surface. An atom that leaves through one face is the copy entering through the opposite one.

Velocities are drawn from a Gaussian, the center-of-mass motion is removed, and the rest is scaled to the target. Equipartition gives \(\tfrac12 k_B T\) per degree of freedom, and fixing the total momentum removes three of the \(3N\), so \(T = \sum v^2 / (3(N-1))\):

def temperature(v):
    return np.sum(v**2) / (3 * (N - 1))

n_cells = round((N / 4) ** (1 / 3))                         # cells per side
assert 4 * n_cells**3 == N, f"N = {N} is not 4 n³"
a = L / n_cells                                              # fcc cell side
basis = np.array([[0, 0, 0], [0.5, 0.5, 0], [0.5, 0, 0.5], [0, 0.5, 0.5]])
cells = np.indices((n_cells,) * 3).reshape(3, -1).T          # cell origins
r_lattice = a * (cells[:, None, :] + basis[None, :, :]).reshape(-1, 3)

v = rng.normal(size=(N, 3))
v -= v.mean(axis=0)                                          # no drift of the whole box
v *= np.sqrt(T_TARGET / temperature(v))

print(f"positions {r_lattice.shape}, nearest neighbors at {a / np.sqrt(2):.3f} σ")
print(f"T = {temperature(v):.4f}, largest total momentum component {np.abs(v.sum(axis=0)).max():.1e}")

positions (256, 3), nearest neighbors at 1.190 σ
T = 0.7500, largest total momentum component 1.8e-14

The neighbors on the lattice sit at 1.19 σ, a little beyond the minimum of the potential at \(2^{1/6}\sigma = 1.12\,\sigma\), and the total momentum is \(1.8 \times 10^{-14}\), rounding error against single velocities near 0.9.

Step 2: Compute forces with the minimum image and a cutoff

In a large system most pairs are far apart and feel almost nothing: \(V(2.5\,\sigma) = -0.016\,\varepsilon\), under 2 % of the well depth. Simulation codes therefore ignore pairs beyond a cutoff \(r_c = 2.5\,\sigma\) and find the rest with neighbor lists, so the cost grows like \(N\) instead of \(N^2\). With 256 atoms all 32,640 distances are cheap, so this code computes them anyway, and the cutoff stays because it defines the standard truncated Lennard-Jones liquid.

Of all the periodic copies of atom \(j\), atom \(i\) feels the nearest one. Per component that is the difference minus the box side times the rounded number of box sides in it, the minimum image. Atoms at \(x = 0.2\) and \(x = 6.5\) are 0.43 σ apart across the face:

dx = 6.5 - 0.2
print(f"{dx:.2f} σ inside the box, {dx - L * np.round(dx / L):+.2f} σ across the face")

6.30 σ inside the box, -0.43 σ across the face

Keep \(r_c < L/2\) (here 2.5 < 3.37) and at most one copy of each atom lies inside the cutoff, so the minimum image is exact. The potential is shifted by \(V(r_c)\) to reach zero at the cutoff, for a reason Pitfall 1 shows. The force on atom \(i\) from atom \(j\) is \(24\varepsilon[2(\sigma/r)^{12} - (\sigma/r)^6]\,(\mathbf r_i - \mathbf r_j)/r^2\):

V_RC = 4 * (RC**-12 - RC**-6)

def forces(r):
    d = r[None, :, :] - r[:, None, :]                  # d[i, j] points from atom i to atom j
    d -= L * np.round(d / L)                           # minimum image
    r2 = np.sum(d**2, axis=-1)
    np.fill_diagonal(r2, np.inf)                       # no force of an atom on itself
    inside = r2 < RC**2
    inv6 = np.where(inside, r2**-3, 0.0)
    U = 0.5 * np.sum(np.where(inside, 4 * (inv6**2 - inv6) - V_RC, 0.0))   # each pair counted twice
    F = -np.sum((24 * (2 * inv6**2 - inv6) / r2)[:, :, None] * d, axis=1)
    return F, U, r2

F, U, r2 = forces(r_lattice)
print(f"V(rc) = {V_RC:.4f} ε")
print(f"lattice: U/N = {U / N:.3f} ε, largest total force component {np.abs(F.sum(axis=0)).max():.1e}")

V(rc) = -0.0163 ε
lattice: U/N = -6.296 ε, largest total force component 1.0e-14

The largest total force component is \(1.0 \times 10^{-14}\), where one nearest neighbor alone pulls with about \(2\,\varepsilon/\sigma\), so Newton's third law survives the minimum image. The function also returns the squared distances, which the loop of Step 3 reuses for the histogram.

Step 3: Melt the lattice and set the temperature by rescaling velocities

The loop is velocity Verlet in place, as in its cheat sheet. Positions are never wrapped back into the box: the minimum image looks only at differences, and the unwrapped positions give each atom's true displacement. Every ten steps the velocities are multiplied by \(\sqrt{T_\text{target}/T}\), which sets the temperature back to 0.75. Every tenth step of a run may also add its pair distances to a histogram, whose bins of 0.02 σ are the shells of \(g(r)\):

iu = np.triu_indices(N, k=1)                           # each pair once
n_bins = int(L / 2 / 0.02)
EDGES = np.linspace(0, 0.02 * n_bins, n_bins + 1)     # bins of 0.02 σ up to L/2

def run(r, v, n_steps, rescale_every=0, sample_every=0):
    r, v = r.copy(), v.copy()
    F, U, r2 = forces(r)
    kin, pot = np.empty(n_steps), np.empty(n_steps)   # per atom, per step
    hist, frames = np.zeros(n_bins), 0
    for k in range(1, n_steps + 1):
        v += 0.5 * DT * F
        r += DT * v
        F, U, r2 = forces(r)
        v += 0.5 * DT * F
        kin[k - 1], pot[k - 1] = 0.5 * np.sum(v**2) / N, U / N
        if rescale_every and k % rescale_every == 0:
            v *= np.sqrt(T_TARGET / temperature(v))
        if sample_every and k % sample_every == 0:
            hist += np.histogram(np.sqrt(r2[iu]), bins=EDGES)[0]
            frames += 1
    return r, v, kin, pot, hist, frames

r_eq, v_eq, kin, pot, _, _ = run(r_lattice, v, 1000, rescale_every=10)
T_eq = kin * 2 * N / (3 * (N - 1))
print(f"lowest T {T_eq.min():.2f} at step {T_eq.argmin() + 1}, mean T over the last 100 steps {T_eq[-100:].mean():.3f}")
print(f"U/N = {pot[-100:].mean():.3f} ε (lattice {U / N:.3f} ε)")
print(f"mean squared displacement from the lattice sites {np.mean(np.sum((r_eq - r_lattice)**2, axis=1)):.2f} σ²")

lowest T 0.39 at step 20, mean T over the last 100 steps 0.756
U/N = -5.158 ε (lattice -6.296 ε)
mean squared displacement from the lattice sites 0.93 σ²

The lattice starts at the bottom of its potential with all the energy kinetic. Each atom swings out of its site like an oscillator, and within a quarter swing about half the kinetic energy has turned into potential: the temperature falls to 0.39 at step 20, 0.1 τ, before the rescalings pull it back. The melting comes later and slower. After 1,000 steps, 5 τ, the mean squared displacement \(\langle |\mathbf r(t) - \mathbf r(0)|^2 \rangle\) is 0.93 σ². A crystal holds its atoms to a few hundredths of σ²; this one has melted.

Rescaling is how you reach a temperature, not how you measure at it. Each rescaling adds or removes energy by hand, so the run is no longer the motion Newton's equations give and its energy record checks nothing. It also pins the temperature where 256 free atoms would let it fluctuate. Switch it off before measuring.