← Lab Notes

Knowing when a molecular dynamics run is equilibrated

~9 min read · MD, equilibration, sampling, LAMMPS

Every graduate student who runs their first molecular dynamics simulation asks the same question, and almost no textbook answers it directly: when is the run equilibrated? "It finished without crashing" is not the answer. Neither is "I ran it for a nice round number of nanoseconds." Averaging over data the system collected while it was still relaxing toward equilibrium quietly biases every number you report. This is the checklist I actually use before I trust an average.

1. Watch the right observables, not just the energy

Equilibration is ensemble-dependent, so monitor the quantities that are free to move in your ensemble:

• Potential energy — the first thing to plateau in almost every run. Necessary, but never sufficient: a system can have a settled energy while its density or structure is still drifting.
• Temperature — should fluctuate around your setpoint. A slow drift means your thermostat coupling is too weak or your timestep is too large.
• Pressure — noisy by nature (its fluctuations scale poorly with system size), so judge its running average, not the instantaneous trace.
• Density / box volume — the decisive one for an NPT run. Cell volume relaxing to a plateau is usually the slowest and most honest equilibration signal you have. For adsorption work, the guest loading itself is the observable that must plateau.

A run is not equilibrated until all of the observables that can move have stopped drifting. The slowest one sets the clock.

2. The discard rule: don't eyeball it, cut it

Once a property has plateaued, you still have to decide how much of the front of the trajectory to throw away. Eyeballing is fine for a sanity check but not for a reported number. Two defensible methods:

Reverse cumulative averaging. Compute the average of the last 10%, 20%, 30%… of the trajectory. Where those averages stop changing within their error bars, you've found a safe cut point — everything before it is equilibration.

Block averaging also gives you the uncertainty you need to report at the same time:

# discard the first N steps, then average the rest in blocks
fix eqcheck all ave/time 100 50 20000 v_myTemp v_myPE v_myDensity &
    file thermo_blocks.dat ave running

Split the production segment into 5–10 contiguous blocks, average each block, and take the standard error across block means. If the block means themselves trend up or down across the trajectory, you discarded too little — the plateau you thought you saw is still sloping.

3. Energy drift is a correctness check, not an equilibration check

In an NVE run the total energy should be conserved. A steady drift in total energy over the trajectory is not a sign of ongoing equilibration — it's a sign your integration is broken: timestep too large, cutoffs too aggressive, or a thermostat left on where it shouldn't be. Quantify it:

# fractional energy drift over the run — want |drift| ≪ 1e-3 per ns
variable e0 equal $(etotal)
variable drift equal (etotal-${e0})/${e0}
thermo_style custom step temp etotal v_drift

If the drift is large, no amount of extra runtime fixes it — you have to shorten the timestep (or tighten the neighbor list) and rerun. Equilibrate only after the integrator conserves energy in a short NVE test.

4. Is the system actually sampling? Check the MSD

A settled energy can hide a system that is stuck. The mean-squared displacement tells you whether molecules are exploring configuration space or frozen in place:

compute msd all msd
fix msdout all ave/time 1000 1 1000 c_msd[4] file msd.dat

• A linear MSD in time means the system is diffusing — it is genuinely visiting new configurations, and time-averages will converge.
• A plateauing MSD means the molecules are caged. That is correct physics for a solid or a tightly adsorbed film, but for a fluid it means you are trapped in one basin and your averages describe that basin, not the ensemble.

You can watch exactly this distinction between a diffusing and a caged fluid, live, in the simulation playground — raise the temperature and the radial distribution function melts from sharp coordination shells toward the featureless g(r) → 1 of a gas.

5. Structural convergence: when the RDF stops moving

The final check is structural. Compute the radial distribution function over the first and second halves of your production window separately. If g(r) is superimposable between the two halves — same peak heights, same positions — the structure has converged. If the second-half peaks are still sharper (or the coordination number is still climbing), the system is still ordering and you need more time.

For adsorption specifically, the analogous check is the isotherm point: the guest uptake at fixed conditions must reach a stationary mean, and its block-to-block variation must be small compared to the differences you care about between state points.

6. The checklist

Before I average anything, all of these must be true:

• Every ensemble-free observable (energy, temperature, pressure, density/loading) has reached a stationary plateau.
• A short NVE test conserves total energy — negligible drift per nanosecond.
• The discard point is chosen by reverse-cumulative or block averaging, not by eye.
• Halving the production window leaves every reported average inside its error bar.
• The MSD is linear (fluid) or a physically-expected plateau (solid/adsorbed) — not accidentally caged.
• g(r) over the first and second halves of production is superimposable.
• Every reported number carries a block-averaged uncertainty.