Molecular Dynamics × Active Learning
Learning a Lennard-Jones equation of state at the fewest simulations — by using calibrated uncertainty to decide what to run next, and how long to run it.
Abstract. A molecular-dynamics (MD) simulator is an expensive, noisy oracle whose noise can be paid down by running longer. This project treats the choice of which thermodynamic states to simulate — and for how long — as sequential experimental design, scored against an analytic equation of state so that error is measurable everywhere. On a supercritical Lennard-Jones fluid the naive variance-seeking rule fails exactly as predicted; a decision-theoretic rule is the most accurate and the most reliable; and a cost-aware policy spends long, precise runs only on the few states that need them.
| \(T^*,\ \rho^*\) | reduced temperature and number density, the two inputs |
| \(P^*,\ U^*\) | reduced pressure and potential energy per particle |
| \(\sigma,\varepsilon,m\) | LJ length, energy, mass scales (all set to 1) |
| \(r\) | pair separation |
| \(N,\ V\) | particle count (864) and box volume |
| \(\ell\) | production length of a run (steps) — a decision variable in Stage 2 |
| \(M,\ N_\text{eff}\) | recorded frames; effective (independent) samples |
| \(g\) | statistical inefficiency, \(N_\text{eff}=M/g\) |
| \(\sigma_n^2(\mathbf x)\) | observation-noise variance (aleatoric) |
| \(s_\text{epi},s_\text{tot}\) | epistemic and total predictive std of the surrogate |
| \(K(\mathbf x)\) | noise coefficient, \(\sigma_n^2=K/\ell\) |
| MD | molecular dynamics |
| LJ | Lennard-Jones (the pair potential) |
| NVT | fixed particles, volume, temperature ensemble |
| EOS | equation of state |
| KN | Kolafa–Nezbeda EOS (the analytic ground truth) |
| GP | Gaussian process (the surrogate) |
| LHS | Latin-hypercube (space-filling baseline) |
| ALC / IMSE | active-learning Cohn / integrated mean-square error |
| FCC | face-centred cubic (initial lattice) |
| PBC | periodic boundary conditions |
An MD simulator answers a physical question by integrating Newton's equations for many particles. Treated as a data source it has three awkward properties.
It is expensive. Each query costs seconds to hours, so the input space cannot be swept densely.
It is noisy, and unevenly so. Every observable is a finite-time average of a fluctuating, correlated signal. The size of that noise varies by more than an order of magnitude across the input domain.
The noise is a control variable. Run longer and it shrinks as \(1/\sqrt{N_\text{eff}}\). So how precisely to measure a point is a decision, weighed against compute — not a fixed property of the point.
The last two together are the whole problem: an oracle where you choose both where to query and how precisely. That is experimental design, not curve fitting.
A single-site Lennard-Jones fluid: structureless particles in a periodic box, interacting only through one pair potential,
\[ U(r) = 4\varepsilon\left[\left(\tfrac{\sigma}{r}\right)^{12} - \left(\tfrac{\sigma}{r}\right)^{6}\right]. \]The \(r^{-12}\) term is repulsion, the \(r^{-6}\) term attraction; the well sits at \(r=2^{1/6}\sigma\) with depth \(\varepsilon\). The potential is truncated at \(2.5\sigma\), with analytic long-range corrections applied to both energy and pressure.
Reduced units set \(\sigma=\varepsilon=m=1\), so the entire phase diagram collapses onto two dimensionless knobs: temperature \(T^*\) and density \(\rho^*=N\sigma^3/V\). We hold \(N=864\), \(V\), and \(T^*\) fixed — the NVT ensemble — with a Langevin thermostat coupling the system to the target temperature. The cubic box has side \(L=(N/\rho^*)^{1/3}\sigma\) and wraps under periodic boundaries; particles start on an FCC lattice and melt.
The pressure follows from the virial — the ideal-gas term plus the average of every pair force projected onto its separation,
\[ P = \rho\,k_B T + \frac{1}{3V}\Big\langle \sum_{i<j}\mathbf r_{ij}\cdot\mathbf f_{ij}\Big\rangle. \]Potential energy per particle \(U^*\) comes from the same trajectory; self-diffusion \(D^*\) is a later, harder estimator built on a mean-squared-displacement fit.
For each observable we want the response surface over the \((T^*,\rho^*)\) plane. The reason to start with a Lennard-Jones fluid is that this surface is known analytically — the Kolafa–Nezbeda equation of state, accurate to about a percent — so the surrogate's error can be evaluated everywhere on a dense grid rather than on a small held-out set.
One restriction: below the critical temperature \((T^*_c\approx1.31)\) an NVT box phase-separates and pressure stops being a smooth function of the inputs. We stay supercritical, \(T^*\ge1.35\), which keeps the hard near-critical curvature and drops the ill-posed region.
A simulation produces a trajectory; an observable is a time average that, by ergodicity, estimates the ensemble average,
\[ \langle A\rangle \approx \bar A = \frac1M\sum_{i=1}^{M} A(t_i). \]Two facts make the obvious error bar wrong. First, the opening transient is not drawn from equilibrium and must be discarded at a detected cutoff. Second, consecutive frames are correlated, so \(M\) frames are worth far fewer than \(M\) independent samples.
The variance of the mean is inflated by the statistical inefficiency \(g\):
\[ \operatorname{Var}(\bar A)=\frac{\sigma_A^2}{N_\text{eff}},\qquad N_\text{eff}=\frac{M}{g},\qquad g = 1 + 2\sum_{k\ge1}\rho(k), \]where \(\rho(k)\) is the autocorrelation at lag \(k\). Typically \(g\approx5\text{–}50\); reporting \(\sigma_A/\sqrt{M}\) instead makes the error bar 2–7\(\times\) too small. That number feeds the decision layer directly, so a wrong \(\sigma\) here corrupts everything downstream silently. The estimator was therefore built and validated first, against synthetic series of known \(g\).
The fluctuation size \(\sigma_A\) is itself a strong function of state — small in the dilute gas, large in the dense liquid (Fig 4). The surrogate must accept a per-point noise variance; a model with one global noise level would average away exactly the structure that makes the problem interesting.
The code is three layers with narrow interfaces. The oracle maps a run configuration to a trajectory; the estimator maps a trajectory to an observable with honest uncertainty; the decision layer maps all observations so far to the next query. The surrogate never touches trajectory data — only estimated values and their variances.
Each observable is modelled with a Gaussian process. Each datum is a noisy observation of a latent surface,
\[ y_i = f(\mathbf x_i) + \varepsilon_i,\qquad \varepsilon_i\sim\mathcal N\!\big(0,\ \sigma_n^2(\mathbf x_i)\big). \]The per-point noise enters as the GP's diagonal (not as a kernel term), which keeps it out of the latent posterior. The predictive variance then splits cleanly,
\[ \underbrace{s_\text{tot}^2(\mathbf x)}_{\text{total}} = \underbrace{s_\text{epi}^2(\mathbf x)}_{\text{reducible}} + \underbrace{\sigma_n^2(\mathbf x)}_{\text{irreducible}}. \]Epistemic variance shrinks with more data and vanishes at an observed point; aleatoric variance is the measurement noise at the chosen run length, and only a longer run shrinks it. Keeping the two apart is the whole game.
A space-filling design that ignores the surrogate — the honest, non-adaptive baseline.
\(a(\mathbf x)=s_\text{tot}^2\). Goes where predictive variance is largest — which, under heteroscedastic noise, is wherever the noise is largest. Expected to fail.
\(a(\mathbf x)=s_\text{epi}^2\). Goes where the function is least known, ignoring irreducible noise.
Scores a candidate by how much observing it would reduce integrated variance everywhere, not by how uncertain it is itself — a decision-theoretic criterion.
Domain \(T^*\in[1.35,3.0]\), \(\rho^*\in[0.05,0.9]\). Each strategy runs from a shared 8-point design; the score is the integrated absolute error against KN over a dense grid, \(\int|f_\text{GP}-P^*_\text{KN}|\), averaged over three seeds.
Worst at the end and stops improving past ~32 runs: it resamples the noisiest region, where total variance is dominated by irreducible noise that never shrinks. Reproduced across all three seeds.
Best at small budgets. On a smooth 2-D surface a space-filling design is close to optimal; adaptivity has little to exploit early.
Starts worst (it front-loads the low-noise, high-pressure corner), crosses over around 31 runs, and finishes lowest and with the tightest band by far — roughly 4\(\times\) more repeatable than LHS.
Tracks LHS throughout and ends just behind ALC. Ignoring the noise term keeps it out of the trap.
The single number hides the geography. Almost all error is in the high-density, high-temperature corner, where pressure is steepest. A uniform design under-resolves it; ALC concentrates there.
A floor everyone hits. The 864-particle, cutoff simulations differ from the full-potential EOS by a small fixed amount (~0.008). Every strategy piles up against it, so the contest is how fast you reach the floor, not how far below it you go. This — not a bug — is why the active-learning margin here is real but modest.
Stage 2 lets the policy choose each run's length \(\ell\) under a fixed compute budget. Longer runs cost more (\(c_0+c_1\ell\)) but return lower noise (\(\sigma_n^2=K(\mathbf x)/\ell\)). Maximizing variance-reduction-per-cost has a closed-form optimal length,
\[ \ell^\star(\mathbf x)=\sqrt{\frac{K(\mathbf x)\,c_0}{s_\text{epi}^2(\mathbf x)\,c_1}} \, . \]Run longer where the noise \(K\) is large; shorter where the function is already uncertain, since a cheap noisy sample there already tells you most of what a precise one would.
Against a fixed-length baseline at equal budget, the cost-aware policy ran 56 shorter simulations (mean ~3,300 steps, adaptively 1,500–7,900) versus 48 at 5,000, for slightly lower error at slightly lower cost. Spatially it banked ~2\(\times\) more production in the high-density corner than in the quiet majority. The answer to the open question — many cheap samples or few precise ones — is both, placed deliberately.
The headline that did not appear is itself the result: on a smooth 2-D surface, against a strong space-filling baseline and a finite-size floor, active learning's margin is modest. What holds cleanly is sharper: the naive rule fails in a specific, predicted way, and the principled rule is both the most accurate and — the property that matters on a one-shot campaign — the most reliable.
The regimes where adaptivity should dominate are the ones this benchmark lacks, and they set the direction: higher input dimension, where space-filling degrades quickly; a localized target such as tracing the coexistence-dome boundary rather than fitting a smooth field; and simply a larger budget with longer runs, to separate the strategies before the finite-size floor. Cost-aware allocation (§6) is the first step into that territory — and its length map is the first sign the policy reasons about precision the way a person spending compute would.