Lab: barostats by hand

The setup

Same 9 atoms in a 2D periodic box. The temperature is held at T₀ = 0.72 by a Langevin thermostat on the atoms (Lab 3.7 showed it gives the right fluctuations for any system). Now the box side L is allowed to change, so that the pressure is P₀ = 3.0.

Why P₀ = 3.0 and not 0? In two dimensions, Lennard-Jones atoms at T = 0.72 are above the critical temperature: at zero pressure they would fly apart like a gas and the box would grow without limit. A pressure of 3 keeps them a dense fluid.

Step 1

Pressure of one configuration, by hand

Why: A barostat can only react to a pressure it can measure. In a simulation there are no walls to push on, so the pressure is computed from two things we already have: how fast the atoms move (kinetic part) and how hard they push or pull on each other (virial part).

Now let the atoms move in the same box (NVT, 290 time units):

What the numbers say: The frozen lattice says −1.343 (the box wants to shrink); the moving atoms say +1.074 (it wants to grow). A single configuration is not "the pressure". Worse, the instantaneous pressure jumps by more than its own average every step. A barostat that reacted fully to each step would shake the box violently, so every barostat responds slowly, through a coupling time τ_P or a piston mass W. It also means that a pressure from a short run is unreliable: always quote averages over long runs.

Step 2

How much should the box breathe?

Why: At constant pressure the box area is not fixed; it fluctuates. As for temperature in Lab 3.7, we need the correct amount before we can judge any barostat. It follows from the equation of state P(A): how the average pressure depends on the box area.

Run NVT at a series of areas (each 300 time units):

What the numbers say: Where the box pressure is above P₀ the box tends to grow, where it is below it tends to shrink, and the exponential tells us by how much: this gives the full distribution of areas, the amber curve in the widget. The box should fluctuate by σ_A = 0.564 (5 % of its area) around ⟨A⟩ = 11.006. Note that the curve is lopsided: P(A) is steep at small A and flat at large A, so the average area is larger than the area where P = P₀ (about 10.83). This is the test for every barostat below, and the size of the fluctuation is how compressibilities are computed from simulations.

Step 3

Berendsen barostat

Why: The simplest idea, the same as the Berendsen thermostat: scale the box (and all coordinates) a little each step, in the direction that brings P towards P₀. The user supplies a coupling time τ_P and a guess of the compressibility β.

What the numbers say: Merit: smooth and safe, it never overshoots, so it is a good way to bring the density to the right value (step 6). Limitation 1: the box fluctuates by 0.247 (τ_P = 0.5) or 0.110 (τ_P = 5) instead of 0.564, so any compressibility computed from it would be 5 to 25 times too small (σ² enters). Limitation 2: the average area is wrong too, 10.83 instead of 11.01, a 1.6 % error in density. Berendsen holds the box near the point where P = P₀ and never explores the lopsided tail of step 2. Weaker coupling makes the fluctuation smaller, not better. Not for production.

Step 4

C-rescale: Berendsen plus the missing randomness

Why: Exactly the fix of the Bussi thermostat, applied to the box (Bernetti and Bussi, 2020). Keep Berendsen's smooth pull, add a random kick of the right size, and a small correction T₀/A. Working with ln A makes the update a simple sum.

What the numbers say: Merit: the noise is 4.5 times larger than the drift in a single step, yet on average the drift wins, so equilibration is as smooth as Berendsen's, and the fluctuation (0.587) and average (11.04) now agree with the correct values within the noise of a 700-time-unit run. It is first order (no piston to oscillate) and works for any τ_P. Limitation: like Bussi, it is stochastic, so trajectories are reproducible only with the same seed, and the compressibility β is still a user guess (it sets the speed, not the result). GROMACS: pcoupl = C-rescale (available since version 2021).

Step 5

The piston: Parrinello–Rahman (Andersen in 1 dimension)

Why: Instead of nudging the box, give it a mass W and let it move like a piston under Newton's law: the force on the piston is the pressure difference. This is deterministic, like Nosé–Hoover, and with a correctly thermostatted system it samples the right volumes. The new question is how heavy the piston should be.

What the numbers say: Merit: all three masses give the correct fluctuation (0.56–0.58) and average (11.05), and the dynamics are smooth and deterministic. The hand formula gets the oscillation period within 20 % (the rest is because P(A) is not a straight line). The mass only sets how fast the box breathes: W = 0.01 oscillates every 58 steps, which is too close to the atomic motions and disturbs them; W = 2 needs about 5 time units per breath, so a run must be many periods long to sample the volume. A common starting point: LAMMPS suggests a damping time of about 1000 time steps. The real limitation appears when you start far from equilibrium (next step). The full Parrinello–Rahman method lets the box shape change too (solids, phase transitions); MTK (LAMMPS fix npt) is the version with Nosé–Hoover chains and exact equations.

Step 6

Starting far away: equilibration

Why: Production barostats are judged by the fluctuations; equilibration barostats by how safely they reach the target. We start from a box that is far too big, A = 14 (P ≈ 1.1, well below 3.0), and watch the area.

What the numbers say: Berendsen glides smoothly into the right range within 1–2 time units and stays there (with the wrong, too-small fluctuations). C-rescale is just as fast and then fluctuates correctly. The light piston (W = 0.1) gathers speed and overshoots to A = 9.99 within half a time unit, where the pressure is more than twice the target; in a real system with a big initial mismatch this can crash the run or melt a solid. The heavy piston (W = 2) swings back to 12.76 at t = 5 and takes many periods to settle. Practical rule: equilibrate with Berendsen or C-rescale, then switch to a piston for production, or simply use C-rescale for both.

Step 7

What the cutoff leaves out of the pressure

Why: Our forces stop at r_c = 1.5, but the real attraction continues beyond it. Every pair beyond the cutoff would pull the box inward a little. Because a barostat sets the box size from the pressure, a missing piece of pressure becomes a wrong density.

What the numbers say: With our short cutoff the missing attraction is −1.157, a third of the target pressure, and the box comes out 4.2 % too large. Even the standard 3D choice r_c = 2.5 misses −0.68 in a dense liquid, which is comparable to the pressure itself. Merit of the tail correction: it is one formula, added once. Limitation: it assumes the atoms beyond r_c are spread evenly, which fails at interfaces (liquid–vapour, membranes); there use a longer cutoff or a long-range method (LAMMPS pair_modify tail yes, GROMACS DispCorr = EnerPres).

Summary