Lab: analysis by hand

The setup

Two systems. Part A (steps 1–4) uses the 9-atom 2D Lennard-Jones box of Labs 2 and 3 (reduced units, T = 0.72, a gentle Bussi thermostat so the dynamics are realistic), plus a few made-up numbers small enough to do by hand. Part B (steps 5–7) uses 2 cations and 2 anions moving in 1D (nm, ps), and a Brownian model of 100 ions.

All positions used for displacements must be unwrapped: the atom is followed across the periodic walls instead of being put back in the box (step 2 shows why).

Step 1

Radial distribution function: counting pairs

Why: g(r) tells how the density around an atom differs from the average: g = 1 means "as in an ideal gas", a peak means a shell of neighbours. It is the bridge to X-ray and neutron scattering, and its first peak gives the coordination number. The recipe is counting, plus one division.

What the numbers say: The factor 2 counts each pair for both atoms; the ideal-gas count 2πr·dr·ρ grows with r, which is why g → 1 at large distance even though there are more and more neighbours. A single frame is useless (36 pairs over 36 bins: mostly 0 or 3), so g(r) is always averaged over many frames, and it converges quickly because every frame contributes N(N−1)/2 pairs. Limitations: only r < L/2 is meaningful (here 1.8), and the bin width trades resolution against noise. In 3D the shell volume is 4πr²dr.

Step 2

Mean squared displacement: all time origins

Why: The diffusion coefficient D follows from how far atoms wander: in d dimensions MSD = 2dDt at long times (Einstein). A trajectory of n frames gives many displacements for each time lag, if we use every frame as a possible starting point (time origin).

What the numbers say: All origins turn one trajectory into many samples, so lags 1–3 already give a steady slope, while the single-origin MSD jumps around (0.09, 0.01, 0.36). The price is the long-lag end: lag 5 has a single sample, so the last part of any MSD curve is noise and must not be fitted. Also fit only the straight (diffusive) part, after the short-time ballistic region. With wrapped coordinates the MSD falls from 1.44 to 0.04 the moment the atom crosses a wall: always unwrap (LAMMPS dumps xu yu zu, GROMACS -pbc nojump).

Step 3

Velocity autocorrelation: the same D a second way

Why: Displacement is the time integral of velocity, so D can also come from how long an atom remembers its velocity: D = (1/d) ∫ ⟨v(0)·v(t)⟩ dt (Green–Kubo). Getting the same D from MSD and VACF is a strong check of the analysis.

What the numbers say: For the liquid, the two routes agree within 1–2 % (0.059–0.061 vs 0.0598): the analysis is consistent. The VACF also shows physics the MSD hides: it goes negative after 0.25 time units, the atom bouncing back from its cage of neighbours. Limitations: the integral must be stopped somewhere, and the noisy tail makes the answer drift (0.0587 → 0.0605); and velocities must be saved very often (here every 10 steps). The toy series shows the other danger: its VACF has not decayed by lag 4, so its "D" is not yet meaningful.

Step 4

Van Hove function: how displacements are distributed

Why: The MSD is just one number per time lag, the average of r². The self part of the van Hove function, G_s(r, t), is the whole histogram of displacements. If every atom diffuses the same way it is a Gaussian; if some atoms jump while others are stuck, it grows tails. The non-Gaussian parameter measures that with two averages.

What the numbers say: Two "hoppers" out of ten carry 80 % of ⟨x²⟩ and push α₂ from about 0 to 0.20. The MSD alone cannot tell this situation (a few fast ions, many stuck) from uniform diffusion; G_s and α₂ can. Our simple LJ fluid stays nearly Gaussian (α₂ ≈ 0.03–0.04). In ionic liquids and glassy electrolytes, α₂ peaks at intermediate times, which is the signature of cages and hopping. The distinct part G_d(r, t), where a different atom is found at time t, is the time-dependent g(r) and shows how neighbours replace each other. Limitation: the tails are rare events, so G_s and α₂ need many more samples than the MSD.

Step 5

Conductivity, Onsager coefficients and the Haven ratio: 2 + 2 ions

Why: Conductivity is a collective property: it depends on how the total charge moves, ΣqΔx, not on how far each ion goes. If a cation and an anion travel together, they carry no current, even though both diffuse. The Nernst–Einstein estimate σ_NE uses only the self terms; the Haven ratio σ_NE/σ measures how much correlations change the result.

What the numbers say: The pair c1–a1 adds 0.115 to the self sum and exactly 0 to ΔM: a neutral pair diffuses but carries no current, so σ is lower than σ_NE and the Haven ratio is above 1. The Onsager decomposition shows where it comes from: L₊₋ > 0 means cations and anions tend to move the same way. It also changes the transference number, from 0.436 (self diffusion only) to 0.345 (with correlations), which matters for battery electrolytes. These are lag-1 numbers for illustration; in practice σ and the L's come from the long-time slopes, exactly as D comes from the MSD slope. The identity L₊₊ + L₋₋ − 2L₊₋ = ⟨ΔM²⟩ is a good check of any implementation.

Step 6

Units: from nm²/ps to S/m

Why: Simulation output comes in nm, ps and e; experiments report S/m or mS/cm. A unit slip here gives errors of factors of 10⁶, so it is worth one careful calculation.

What the numbers say: Every factor is visible: the charge squared, the thermal energy, the ion density and the sum of diffusion coefficients. Typical ionic-liquid numbers (D ~ 10⁻¹¹ m²/s, a few ion pairs per nm³) give a few mS/cm, as measured. Remember the conversion 1 m²/s = 10⁶ nm²/ps, and that σ (from ⟨|ΣqΔr|²⟩) has a 6 (3D) in the denominator just like D from the MSD. Here n counts ion pairs, so each species has density n; for a salt like Li₂X the charges and densities must be counted species by species.

Step 7

Statistics: why collective quantities converge so slowly

Why: D is an average over every ion; σ comes from one quantity per time origin, the total charge displacement of the whole box. With N ions, the self average has about N times more independent samples. How much noisier is σ in practice?

What the numbers say: The same runs give D to 2.4 % but σ only to 26 %, about √N = 10 times worse, and a single run can report a Haven ratio anywhere from 0.91 to 3.05 when the true value is 1.52. To make σ as precise as D you need about N times more sampling (longer runs or several independent runs; a bigger box does not help σ, because ΔM is still one number per time origin). That is why conductivities and Onsager coefficients from short simulations often scatter, and why papers should report error bars from independent runs or block averaging. The self quantities (D, t₊ from self diffusion) are the cheap ones; the collective ones carry the correlations but cost far more.

Summary