Lesson 6.1 · 14 min

Radial distribution function

The radial distribution function g(r) tells you how likely you are to find a neighbour at distance r, compared with a random (ideal-gas) arrangement. It is the first thing to compute from any trajectory: it shows the structure, gives coordination numbers and can be compared directly with X-ray and neutron scattering.

What g(r) means and how to compute it from a trajectory (histogram + ideal-gas normalization)
Reading g(r): excluded core, neighbour shells, gas vs liquid vs solid
Coordination numbers, the potential of mean force, and the link to scattering

The idea

Sit on one atom and ask: how many neighbours are at distance r? Count the atoms in a thin spherical shell between r and r + dr, and divide by the number you would expect if the atoms were scattered at random with the same density ρ (an ideal gas). That ratio is the radial distribution function:

g(r) = (neighbours in the shell) / (ρ · 4πr² dr)

Average over every atom as the centre and over many frames, and g(r) becomes smooth.

Reading it (code tab, Part 2; Lennard-Jones liquid near its triple point):

g = 0 for r < 0.9σ: atoms cannot overlap.
A tall first peak (g ≈ 3 at r ≈ 1.09σ): the first shell of neighbours, three times more likely than at random.
Smaller peaks at about 2σ and 3σ, then g → 1: beyond a few molecular diameters a liquid has no memory of the central atom.

Gas, liquid, solid. A dilute gas has one broad peak, close to e^−u(r)/kT. A liquid has a few decaying peaks. A crystal has sharp peaks at the lattice distances that never die out (Part 3).

Coordination number. Integrating ρ g(r) 4πr² up to the first minimum counts the neighbours in the first shell: about 12 for the LJ liquid and exactly 12 for fcc. For battery electrolytes this is how you count how many anion oxygens or solvent molecules surround a Li⁺ ion, using the Li–O g(r).

Potential of mean force. w(r) = −kT ln g(r) is the effective free energy between two atoms, averaged over everything else (Part 4). It differs from the bare pair potential because of the surrounding atoms.

Experiments. X-ray and neutron scattering measure the structure factor S(k), which is a Fourier transform of g(r) − 1. That makes g(r) one of the most direct checks of a force field.

Key idea. g(r) = local density at distance r / average density. Zero inside the core, peaks for the neighbour shells, 1 far away. Integrate the first peak to get coordination numbers.

The algorithm

choose dr (about 0.01-0.02 σ, or 0.002 nm); bins up to r_max ≤ L/2
h[:] = 0
for each saved frame:
    for each pair i < j:
        r = |minimum image of (x_i - x_j)|
        if r < r_max: h[floor(r / dr)] += 2        # the pair counts for both atoms
for each bin k (centre r_k):
    shell = 4/3 π ((r_k + dr/2)³ - (r_k - dr/2)³)
    g[k] = h[k] / (n_frames · N · ρ · shell)
n_1 = Σ_{r_k < r_min} 4π ρ r_k² g[k] dr               # coordination number

Never go beyond half the box length: past L/2 the minimum image no longer sees every neighbour.
For two species A and B, count only A–B pairs and divide by N_A ρ_B; the coordination number of B around A uses ρ_B.
Equilibrate first, and use uncorrelated frames (every 10–100 steps). The first peak converges fast; the tail and dilute species need many frames.
Tools: GROMACS gmx rdf, LAMMPS compute rdf + fix ave/time, VMD, MDAnalysis InterRDF, and TRAVIS.