Diffusion & Ionic Transport¶

This page is a self-contained, textbook-style introduction to how MolPy turns an equilibrium molecular-dynamics (MD) trajectory into transport coefficients — diffusion coefficients, Onsager phenomenological coefficients, and ionic conductivity. It starts from the random walk and builds up to the collective correlation functions used in modern electrolyte analysis. No background beyond undergraduate statistical mechanics and a little linear algebra is assumed.

It is the conceptual companion to the Dielectric Spectroscopy guide: that page derives the frequency-dependent response and the two conductivity routes in depth; this page focuses on diffusion and the displacement picture, and points back to the dielectric page for the spectral machinery.

As with all MolPy analyses, the heavy numerics run in Rust (molrs); the MolPy layer extracts coordinates, unwraps periodic images, builds the collective quantities, and returns typed result objects.

1. The random walk and the Einstein relation¶

A particle in a liquid is kicked around by its neighbours. Over a short time it moves ballistically (it remembers its velocity), but after many uncorrelated kicks its motion becomes a random walk: the direction is forgotten while the spread keeps growing. The natural measure of that spread is the mean-squared displacement (MSD):

For a diffusive process the MSD grows linearly in time, and the slope defines the self-diffusion coefficient \(D\) (Einstein, 1905):

The factor \(1/6 = 1/(2d)\) with \(d=3\) counts the three independent directions the walker can spread into.

1.1 The three regimes¶

A real MSD curve has three parts, and only the middle one is physical diffusion:

• Ballistic (short \(\tau\)): \(\mathrm{MSD}\propto\tau^2\) — the particle still moves at roughly constant velocity.
• Diffusive (intermediate \(\tau\)): \(\mathrm{MSD}\propto\tau\) — the linear regime; fit the slope here.
• Noisy (long \(\tau\)): few time origins remain, so the estimate is dominated by statistical scatter.

Choosing the linear window is the single most important judgement call in a diffusion calculation. The transport computes return the full correlation curve so you can inspect it before fitting; IonicConductivity exposes a fraction-based fit window (fit_start_frac/fit_end_frac) for exactly this reason (see §6).

1.2 Averaging over time origins¶

At equilibrium the dynamics are stationary, so every frame can serve as a time origin \(t\). Averaging over all origins (a "windowed" MSD) uses the data far more efficiently than measuring displacement from frame 0 only:

The number of usable origins shrinks as \(\tau\) grows — which is exactly why the long-\(\tau\) tail is noisy.

1.3 Minimum-image unwrapping (the usual trap)¶

MD runs in a periodic box: when an atom leaves one face it re-enters the opposite one, so its stored coordinate jumps by a box length \(L\). A naive MSD built from such coordinates registers a spurious \(L\)-sized displacement at every crossing. The fix is to accumulate minimum-image steps,

valid as long as a particle moves less than \(L/2\) per frame. MolPy does this with Box.delta(p1, p2, minimum_image=True); it is shared by every displacement-based transport compute. (See Dielectric §2.1 for the same argument in the dipole context.)

2. Self vs distinct diffusion: the MDC¶

A single diffusion coefficient hides a lot of physics. In a multi-component system the displacements of different particles are correlated — ions drag their counter-ions, solvent flows around them. The mean displacement correlation (MDC) generalises the MSD to expose this.¹

Self (tag "3") — the ordinary MSD of one species, giving the self-diffusion coefficient \(D^\mathrm{s}_\alpha\):

Distinct (tag "3,4") — the cross-correlation between the displacements of species \(\alpha\) and \(\beta\), giving the distinct-diffusion coefficient \(D^\mathrm{d}_{\alpha\beta}\):

The distinct term is a collective quantity: it is dominated by how the species move together, not by any single particle.

from molpy.compute import MCDCompute

mdc = MCDCompute(tags=["3", "4", "3,4"], max_dt=20.0, dt=0.01)
result = mdc(trajectory)
result.correlations["3"]    # self MSD of species 3, vs lag time
result.correlations["3,4"]  # distinct cross-correlation of 3 and 4

3. Onsager phenomenological coefficients¶

The cleanest way to describe coupled transport in an electrolyte is Onsager's phenomenological coefficients \(\Omega_{\alpha\beta}\) (also written \(L_{\alpha\beta}\)). They are the proper, normalized version of the collective displacement correlation: \(\Omega_{\alpha\beta}\) relates the flux of species \(\alpha\) to the thermodynamic driving force on species \(\beta\).

Define the collective coordinate of a species — the summed (unwrapped) position of all its atoms,

The collective displacement correlation and the Onsager coefficient are then

• The diagonal \(\Omega_{\alpha\alpha}\) is the collective MSD of species \(\alpha\) — it contains the self term plus the like-ion cross terms.
• The off-diagonal \(\Omega_{\alpha\beta}\) (\(\alpha\ne\beta\)) captures cation–anion coupling. A negative value (anticorrelated drift) is the signature of ion pairing.

Onsager returns the correlation curves \(\mathrm{corr}_{\alpha\beta}(\tau)\); you take the long-time slope and apply the \(1/(6 k_B T V N_A)\) prefactor for the coefficient itself.

from molpy.compute import Onsager

ons = Onsager(tags=["1,1", "1,2", "2,2"], max_dt=20.0, dt=0.01)
result = ons(trajectory)
result.correlations["1,2"]  # cation-anion collective correlation L_12(tau)

3.1 From Onsager coefficients to conductivity¶

The ionic conductivity is a weighted sum of the Onsager coefficients over the ion charges \(z_\alpha\):

If the off-diagonal (distinct) terms vanish — ions moving independently — this collapses to the Nernst–Einstein estimate \(\sigma_\text{NE}\) built from the self-diffusion coefficients alone. The ratio \(\sigma/\sigma_\text{NE}\) (the ionicity or Haven ratio) measures how strongly ion correlations suppress (or enhance) conduction — a number you cannot get from single-particle diffusion alone, which is the whole reason the Onsager picture exists.

4. Ionic conductivity: the two equivalent routes¶

Conductivity can be obtained from the collective charge transport directly, without going through individual \(\Omega_{\alpha\beta}\). There are two equivalent routes (a general consequence of the fluctuation–dissipation theorem; see Dielectric §1.3).

4.1 Einstein route — polarization MSD (PMSD)¶

Build the collective charge displacement (a.k.a. the translational dipole) of the ions,

and measure its MSD. Its long-time slope gives the conductivity by an Einstein relation:

from molpy.compute import PMSDCompute, IonicConductivity

# The PMSD curve itself:
pmsd = PMSDCompute(cation_type=1, anion_type=2, max_dt=30.0, dt=0.01)(trajectory)

# The conductivity (Einstein-Helfand), fitted and converted to S/m:
sigma = IonicConductivity(dt=0.01, temperature=298.15, max_correlation_time=1000)(ion_trajectory)
sigma.sigma  # S/m

PMSDCompute returns the curve; IonicConductivity does the fit and unit conversion. See Dielectric §7 for the full derivation, the diffusive-window caveat, and the SI prefactor.

4.2 Green–Kubo route — current autocorrelation (JACF)¶

Equivalently, integrate the autocorrelation of the charge current \(\mathbf{J}(t)=\sum_a q_a\mathbf{v}_a(t)\):

The integrand \(C(\tau)=\langle\mathbf{J}(0)\cdot\mathbf{J}(\tau)\rangle\) is the current autocorrelation function (JACF); the factor \(1/3\) is the Green–Kubo analogue of the Einstein \(1/6\) (one comes from integrating the ACF, the other from differentiating the MSD). The single-particle analogue — the Green–Kubo route to \(D\) via the velocity autocorrelation — is derived in Velocity Autocorrelation & VDOS.

from molpy.compute import JACF

jacf = JACF(cation_type=1, anion_type=2, max_dt=30.0, dt=0.01, temperature=298.15)
result = jacf(trajectory)        # requires per-atom velocities vx, vy, vz
result.jacf                      # <J(0).J(t)>
result.sigma                     # DC conductivity, S/m (the GK integral)
result.sigma_running             # running integral sigma(tau), to check convergence

The Einstein (PMSD) and Green–Kubo (JACF) routes are mathematically identical; in practice the Einstein route is more robust at coarse sampling, while the JACF exposes the current memory directly and lets you watch the integral converge. Dielectric §8 derives the frequency-dependent generalisation \(\sigma(\omega)\).

5. Reading the results¶

Cross-checks. The Einstein (IonicConductivity/PMSDCompute) and Green–Kubo (JACF) conductivities must agree within statistics. The Nernst–Einstein estimate (from MCDCompute self terms) should exceed the correlated conductivity from Onsager/JACF when ion pairing is significant — their ratio is the ionicity.

6. Parameters and hyperparameters¶

6.1 Parameters and their meaning¶

MCDCompute, Onsager and PMSDCompute return raw correlation curves — the \(1/(2d\,\tau)\) slope, the \(1/(6 k_B T V N_A)\) Onsager prefactor, and any unit conversion are applied by you, which is what keeps the fit window under your eyes. Only IonicConductivity and JACF fit and convert internally.

6.2 Hyperparameter effects¶

• Fit window placement. Fitting too early includes the ballistic/caged head where the local exponent \(d\ln\mathrm{MSD}/d\ln\tau \ne 1\) — the extracted slope is systematically biased (check that the log–log slope is ≈ 1 across your window before trusting \(D\)). Fitting too late trades bias for variance: the number of independent time origins falls as \(N_\text{origins} = N_\text{frames} - \tau/\Delta t\), so the statistical error grows roughly like \(\sqrt{\tau/T_\text{traj}}\). Move fit_start_frac/fit_end_frac (or your manual window) and report the spread — for few-carrier collective quantities it is rarely below several percent.
• max_dt vs trajectory length. Lags approaching the trajectory length average over almost no origins; keep max_dt a small fraction of the run (the ≤ n_frames / 5 rule from IonicConductivity is a good default for every curve on this page). Doubling the trajectory beats doubling max_dt.
• Frame spacing dt. The Einstein/MSD routes are robust at coarse dt — displacement is cumulative. The Green–Kubo JACF is not: the current ACF decays in ~0.1–1 ps, so coarse frames under-resolve the integrand and \(\sigma\) drifts. The same sampling-rate trade-offs as for the VACF apply — see vacf.md §6.
• Dimensionality. The prefactors on this page assume \(d = 3\) (\(1/(2d) = 1/6\) Einstein, \(1/3\) Green–Kubo). For quasi-2-D systems (confined films, membranes) apply the \(d = 2\) factors to the raw curves yourself — using the 3-D factor on 2-D motion underestimates \(D\) by \(2/3\).
• Centre-of-mass drift (center_of_mass). Collective quantities (Onsager, PMSD, JACF) are a single realization per species — a net COM drift adds a coherent \(\propto \tau^2\) contamination that never averages out. Remove it (pass the mass map, or pre-clean the trajectory) before quoting off-diagonal coefficients.
• temperature / volume. \(\sigma \propto 1/(V\,T)\): a 3 % error in either rescales the conductivity linearly. Use the production-run averages, not the target thermostat values, for NPT data.
• Green–Kubo integration limit (JACF.max_dt). Quote \(\sigma\) at the plateau of result.sigma_running, never at the last lag — integrating ACF tail noise makes the estimate drift linearly (same plateau rule as the running-integral \(D(\tau)\) in vacf.md §2).

The scattered warnings in §7 are the failure modes of exactly these knobs.

7. Pitfalls checklist¶

1. No unwrapping → boundary crossings inject \(L\)-sized jumps; every MSD is garbage. (MolPy unwraps automatically via Box.delta.)
2. Fitting outside the diffusive window → the ballistic head or the noisy tail biases \(D\) and \(\sigma\). Always inspect the curve first.
3. Too few carriers / too short a trajectory → collective quantities (PMSD, Onsager off-diagonal, JACF) are intrinsically noisy because there is only one collective coordinate per species. Report a range, not a digit.
4. Wrong velocity units in JACF → the current must be in \(e\cdot\)Å·ps⁻¹ (velocities in Å/ps); \(\sigma\) scales linearly, so a unit slip rescales the answer.
5. Ignoring distinct diffusion → quoting only Nernst–Einstein conductivity ignores ion correlations and typically overestimates \(\sigma\).

8. References¶

• A. Einstein, Ann. Phys. 322, 549 (1905) — the diffusion/MSD relation.
• M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, 2^nd ed. (2017) — MSD, time-origin averaging, transport coefficients.
• J.-P. Hansen, I. R. McDonald, Theory of Simple Liquids, 4^th ed. — Green–Kubo relations and the current correlation function.
• D. Frenkel, B. Smit, Understanding Molecular Simulation, 2^nd ed. (2002), §4.4 — Einstein relations for transport coefficients.
• L. Onsager, Phys. Rev. 37, 405 (1931); 38, 2265 (1931) — reciprocal relations and phenomenological coefficients.

See also¶

• Dielectric Spectroscopy — the spectral machinery (\(\varepsilon^*(\omega)\), autocorrelation, FFT) and the full conductivity derivations.
• Pair Persistence — residence times and the survival functions that resolve the pairing contribution to diffusion.
• Compute overview — the Compute → Result pattern.
• API reference: Compute.