Lesson 6.3 · 14 min

Velocity autocorrelation

The velocity autocorrelation function (VACF) asks how long an atom remembers its velocity. Its shape tells you whether atoms fly freely, rattle in a cage or vibrate in a lattice; its integral gives the diffusion coefficient (Green–Kubo), and its Fourier transform gives the vibrational spectrum.

The idea

Take an atom's velocity at some moment, v(0), and again a time t later, v(t). Their dot product is positive if the atom is still moving the same way, negative if it has turned around, and zero if there is no relation. Averaged over all atoms and all starting times this is the velocity autocorrelation function:

C(t) = ⟨ v(t) · v(0) ⟩,   C(0) = ⟨v²⟩ = 3kT/m

Three shapes (code tab, Parts 1 and 2):

Green–Kubo. The diffusion coefficient is the area under the VACF: D = ⅓∫C(t)dt. It is exactly equivalent to the MSD slope (they are two views of the same motion). For the LJ liquid both give D* ≈ 0.031 from the same run. The running integral overshoots, because the positive early part comes first, then settles once the negative dip has been added.

Spectrum. The Fourier transform of C(t) is the vibrational density of states (VDOS): which frequencies the atoms move at (Part 3). A solid has peaks and S(0) = 0; a liquid has a broad band and a nonzero value at ω = 0, proportional to D. For molecules the VDOS shows the bond-stretch and bend frequencies and can be compared with IR, Raman and inelastic neutron spectra (with the right weighting).

Long-time tails: in liquids C(t) decays as t−3/2 at long times (hydrodynamics), so the Green–Kubo integral converges slowly; noise in the tail is the main source of error.

Key idea. C(t) = ⟨v(t)·v(0)⟩ measures velocity memory. Its area gives D (Green–Kubo, same as the MSD), the negative dip shows caging, and its Fourier transform gives the vibrational spectrum.

The algorithm

save velocities every n steps (short interval: a few fs for molecules, ~0.01 τ here)
for each lag k = 0 ... K:
    C[k] = mean over atoms i and origins t0 of v_i(t0 + k) · v_i(t0)
D(t)  = (1/3) * cumulative trapezoid integral of C            # running Green-Kubo integral
D     = plateau value of D(t)
S(ω)  = (2 / C[0]) Σ_k C[k] w[k] cos(ω k Δt) Δt                 # VDOS, w = smooth window
  1. Run in NVE (or with a very weak thermostat) for dynamical properties; strong thermostats change the VACF.
  2. The FFT gives C(t) for all lags at once (Wiener–Khinchin): pad to twice the length, |FFT|², inverse FFT.
  3. For molecules, use atom (or centre-of-mass) velocities depending on what you want: COM for diffusion, atoms for spectra; mass-weight the VDOS for the total.
  4. Tools: MDAnalysis, TRAVIS (power spectra, IR/Raman with dipoles), GROMACS gmx velacc, LAMMPS compute vacf (single origin).