Lab: packing and equilibrating a liquid

The setup

The recipe used by every real setup (Packmol, gmx insert-molecules, LAMMPS create_atoms random):

System: N = 256 LJ atoms, cutoff 2.5σ, k_BT = 0.722 and ρ* = 0.8442 (the liquid near the triple point, the classic test state). Bussi thermostat (τ = 0.5), C-rescale barostat, Δt = 0.005. Argon: σ = 3.405 Å, ε = 0.996 kJ/mol.

Step 1

Box size from the target density

Why: Experiments give density in g/cm³; a simulation needs a number of molecules and a box length. The conversion is one line, but getting it wrong by a few percent shifts every property you compute (pressure most of all).

What the numbers say: Pure liquids with a known density can be built at the right size directly. For a solution, a mixture or a new force field, the density is part of the answer: build at a guess, then let NPT (Step 4) find it. Our test state ρ* = 0.8442 is the same liquid as argon at 1.42 g/cm³.

Step 2

Random positions: why they fail, and why a tolerance only partly helps

Why: The obvious way to fill a box is to put each atom at a random point. In a dense liquid, that makes many atoms land almost on top of each other, and the LJ repulsion grows as r⁻¹².

Insertion with a tolerance: a new atom is refused if it is closer than d to any placed atom (give up after 20 000 refusals in a row).

What the numbers say: Fully random positions give forces 10¹⁵ times too large and the first step sends atoms flying. A tolerance removes the worst clashes, but a tolerance equal to the atom size cannot reach a liquid density: random filling jams long before. That is why Packmol uses a tolerance smaller than the molecules (2 Å) and why a minimization always follows. A lattice start (Lab 2.7) avoids the problem entirely, at the price of having to melt it.

Step 3

Steepest descent by hand

Why: Minimization moves every atom downhill along its force until no force is large. GROMACS uses a simple, robust rule: move the atom with the largest force by a step h, others in proportion; grow h when the energy falls, cut it hard when it rises.

The random 256-atom box, stopped at the GROMACS default emtol = 10 kJ/mol/nm = 10/(0.996/0.3405) = 3.42 ε/σ:

What the numbers say: Steepest descent removed every clash in 160 iterations; the start that exploded now conserves energy. But a minimized structure is not an equilibrium liquid: it sits at T = 0, at an energy (−6.00) far below the liquid at T = 0.722 (−5.20), and velocities given at T lose about half their energy to the potential. That is why minimization is always followed by a thermostatted run, not by a measurement.

Step 4

NVT from two starts, then NPT to find the density

Why: The thermostat brings the temperature to the target within a few τ. Whether the structure follows depends on the start: a lattice has to melt, a minimized random box is already disordered. The FCC order parameter (1 on the lattice, about 1/N in a liquid) shows how far along the melting is.

What the numbers say: The temperature is right after t ≈ 1 in both runs, but the lattice start keeps crystal order until t ≈ 10: temperature alone says nothing about whether the structure has relaxed. NPT started 17% too dilute and found the density to 0.25% (1 error bar) within t ≈ 10. This is how the density of the NaCl solution in Step 1 would be found.

Step 5

When is it equilibrated? Three tests

Why: Everything before equilibrium must be thrown away, and everything after kept. Throwing away too little biases the average; too much wastes data. The tests below replace “it looks flat”.

(i) Choose t₀ to maximize the number of independent samples after it: N_eff = (n − n₀)/g, with g the statistical inefficiency of the data after t₀ (Lab 6.8).

(ii) Half-split: compare the means of the first and second halves of what is kept.

(iii) Independent replicas (different starts and seeds), t > 12 kept:

What the numbers say: Different quantities equilibrate at different times: the energy of the lattice start at t₀ ≈ 8, its structure at t₀ ≈ 10, the random start already at 2. Discard by the slowest quantity you care about. The half-split test can be fooled by the very trend it is meant to detect, so look at the trace as well. Replicas from different starts are the strongest test: if they agree within their error bars, no memory of the start is left.

Summary