Lesson 6.5 · 16 min

Ionic conductivity

Ionic conductivity measures how much current ions carry in an electric field. From an equilibrium MD run it follows from the fluctuations of the total charge displacement (Einstein–Helfand) or of the electric current (Green–Kubo). Comparing it with the Nernst–Einstein estimate from self-diffusion reveals ion pairing and correlated motion.

The Einstein–Helfand and Green–Kubo formulas for conductivity from an equilibrium run
The Nernst–Einstein approximation, ionicity and the Haven ratio, and why ion pairs lower σ
Why conductivity is so noisy, and how to convert to S/m and mS/cm

The idea

Apply a small electric field E and ions drift: cations one way, anions the other. The current density is J = σE, where σ is the ionic conductivity. You could simulate exactly that (non-equilibrium MD), but linear-response theory gives σ from an ordinary equilibrium run, by watching how the ions' charge spreads out spontaneously.

Einstein–Helfand. Add up every ion's displacement weighted by its charge, M(t) = Σ q_iΔr_i(t). Its mean square grows linearly in time, and the slope gives σ, exactly as the MSD gives D:

σ = (1 / 6VkT) · d/dt ⟨|Σ q_i Δr_i(t)|²⟩

The equivalent Green–Kubo form integrates the autocorrelation of the current J = Σq_iv_i.

Nernst–Einstein. Expanding the square gives self terms (each ion with itself, i.e. its own MSD) and cross terms (pairs of different ions). Keeping only the self terms gives the Nernst–Einstein estimate from the self-diffusion coefficients (Lesson 6.2):

σ_NE = (e² / VkT) (N₊q₊²D₊ + N₋q₋²D₋)

The cross terms describe correlated motion. If a cation and an anion move together as a neutral pair, both still diffuse (they count in D), but together they carry no current. Then σ < σ_NE (code tab, Part 1).

Ionicity (or degree of dissociation) = σ/σ_NE; its inverse is the Haven ratio H_R = σ_NE/σ.
For many ionic liquids and concentrated electrolytes, measured ionicities are roughly 0.5–0.8. Cation–cation and anion–anion correlations matter too, not only ion pairs: the full picture needs the Onsager coefficients (Lesson 6.6).

Statistics. D averages over N ions, but M(t) is one collective number per frame, so σ converges about N times more slowly (Part 2). Use long runs, many independent replicas, and block averaging, and always quote an error bar. In molecular systems remember that σ is also affected by the force field (for example, polarizability or charge scaling, Lesson 5.4).

Key idea. σ = slope of ⟨|Σq_iΔr_i|²⟩ / 6VkT. Keeping only self terms gives Nernst–Einstein. Their ratio (ionicity) measures correlated ion motion. Conductivity needs much more sampling than diffusion.

The algorithm

save UNWRAPPED positions every n steps
M(t) = Σ_i q_i u_i(t)                                   # one 3-vector per frame
for each lag k:  MSD_M[k] = < |M(t0 + k) - M(t0)|² >_t0       # all time origins
σ    = slope of MSD_M (linear region) / (6 V kB T)
σ_NE = (N+ q+² D+ + N- q-² D-) / (V kB T)                    # D from the per-species MSD
ionicity = σ / σ_NE;   Haven ratio = σ_NE / σ
errors: blocks or independent replicas (≥ 5);  check that the slope is linear
SI units: σ [S/m] = σ_reduced · e² / (kB T · length · time) with your MD units

Use the same lag-time window for σ and for the D values in σ_NE.
Check that the total charge is exactly zero. Then M(t) does not change when all ions shift together, so any centre-of-mass drift cancels.
Green–Kubo alternative: save velocities often, compute ⟨J(t)·J(0)⟩, integrate, and divide by 3VkT.
Non-equilibrium check: apply a weak field E, measure the steady drift current J, and use σ = J/E (field small enough to be linear).
Tools: TRAVIS, MDAnalysis plus your own script (as here), GROMACS gmx current, LAMMPS via compute of the charge-weighted displacement.