Mean squared displacement and diffusion

The idea

Follow each atom from time 0 to time t and measure how far it has moved, r(t) − r(0). The average of the square of this displacement over all atoms is the mean squared displacement:

MSD(t) = ⟨ |r(t) − r(0)|² ⟩

Because every moment in an equilibrium run is an equally good starting point, we also average over all time origins. This makes the MSD much smoother than using t = 0 only.

Two regimes (code tab, Part 2):

• Short times (ballistic): before an atom hits its neighbours it flies freely, r = vt, so MSD = ⟨v²⟩t² = 3(kT/m)t². On a log–log plot the slope is 2.
• Long times (diffusive): after many collisions the motion is a random walk and the MSD grows linearly, MSD = 6Dt in 3D. On a log–log plot the slope is 1.

In between, a dense liquid shows a "cage" region where atoms rattle among their neighbours. In a solid the MSD stops growing (D = 0); in glasses and polymer electrolytes the cage can last nanoseconds.

Getting D. Fit a straight line to the MSD in the region where the log–log slope is 1, and divide the slope by 6 (by 4 in 2D, by 2 in 1D). For the LJ liquid at ρ* = 0.8442, T* = 0.722 this gives D* ≈ 0.03 (for argon, about 1.6 × 10⁻⁵ cm²/s).

Three classic mistakes:

• Wrapped coordinates. An atom that leaves the box re-enters on the other side, jumping by L. Always use unwrapped coordinates (Part 3).
• Fitting too early or too late. Too early includes the ballistic or cage region; too late uses lag times with few origins and large noise.
• Ignoring the box size. Hydrodynamic interactions with periodic images slow diffusion by about k_BT ξ/(6πηL). That is about 15 % for this small box and still a few percent for typical boxes (Part 4).

For ionic liquids and battery electrolytes, the MSD can stay sub-diffusive for many nanoseconds, so check the log–log slope before trusting a D. Self-diffusion is what pulsed-field-gradient NMR measures.

The algorithm

save unwrapped positions every n steps (or unwrap: add L each time an atom crosses the box)
for each lag k = 1 ... K:
    MSD[k] = mean over atoms i and origins t0 of |r_i(t0 + k) - r_i(t0)|²
check the log-log slope; pick t1..t2 where it is ≈ 1
D = slope of a linear fit of MSD(t) on [t1, t2] / 6
error: split the run into blocks (or use x, y, z separately) -> D ± std / sqrt(n_blocks)
finite size: D_inf = D(L) + 2.837297 k_B T / (6 π η L)   (cubic box; η = shear viscosity)

1. Remove centre-of-mass drift (or check it is zero) before computing the MSD.
2. The lag time should not exceed about a third of the run; long lags have few independent origins.
3. For mixtures compute one MSD per species (Li⁺, anion, solvent). Cross terms and collective motion belong to conductivity (Lesson 6.5).
4. Fast all-origin MSD uses FFTs (cost N_frames log N_frames), as in MDAnalysis EinsteinMSD(fft=True). Other tools: GROMACS gmx msd, LAMMPS compute msd (single origin!), TRAVIS.