Lesson 3.6 · 15 min

Pressure and barostats

Pressure in MD comes from the virial: atoms hitting the walls plus the forces between them. It is very noisy. A barostat lets the box breathe so the pressure is right on average. Berendsen gets the density right but the fluctuations wrong; C-rescale and Parrinello–Rahman/MTTK get both right.

How pressure is computed from the virial, and why the instantaneous value is so noisy
Berendsen vs stochastic cell rescaling (C-rescale) vs Parrinello–Rahman / MTTK
Choosing τ_P, the compressibility, and isotropic vs semi-isotropic vs anisotropic coupling

The idea

Where pressure comes from. Pressure has two parts. The kinetic part is the momentum carried by atoms moving through any surface (the ideal-gas term Nk_BT/V). The virial part comes from forces acting across that surface: repulsion pushes the box outwards, attraction pulls it in. MD computes both every step from the velocities and from Σ r·F over pairs, a quantity we already have from the force loop.

It is noisy. The instantaneous pressure fluctuates strongly, much more than the temperature. For a few thousand water molecules it swings by hundreds of bar from step to step, even when the average is exactly 1 bar. Only long averages mean anything. And like the energy (Lesson 2.5), a cutoff needs a tail correction for the pressure.

Barostats. To simulate at constant pressure (NPT, Lesson 3.2) the box volume becomes a variable. All positions are scaled with the box, so the atoms are not squeezed against one wall.

Berendsen: rescales the box so the pressure relaxes to P₀ with time constant τ_P. Robust and fast for equilibration, and the average density is fine, but the volume fluctuations are wrong (they depend on τ_P).
C-rescale (stochastic cell rescaling, Bernetti and Bussi, 2020): Berendsen plus a correctly chosen random term, exactly as V-rescale is for temperature. Correct NPT fluctuations, still first-order, so no oscillations. A very good default.
Parrinello–Rahman / MTTK: the box has a "mass" and its own equation of motion, like a piston. Correct ensemble and smooth dynamics, but it oscillates if started far from equilibrium; equilibrate with Berendsen or C-rescale first.

How the box may deform. Isotropic: all sides scale together (liquids). Semi-isotropic: the xy plane and z are coupled separately (membranes, interfaces). Anisotropic / fully flexible: each box vector can change, including shape (crystals, phase transitions).

Key idea. P = kinetic + virial, averaged over a long run. Equilibrate with Berendsen or C-rescale; produce with C-rescale or Parrinello–Rahman/MTTK. Choose the coupling (isotropic, semi-isotropic, anisotropic) to match the physics.

The algorithm

One NPT step (thermostat inside the MD step, barostat after it):

MD step with a thermostat
P = (2 K / 3 + W / 3) / V              # W = Σ r·F over pairs

# Berendsen
de = -(beta / tau_P) * (P0 - P) * dt

# C-rescale (adds the right noise; also scale v)
de = -(beta / tau_P) * (P0 - P - kB T / V) * dt
     + sqrt(2 kB T beta dt / (V tau_P)) * gaussian()

mu = exp(de / 3)                       # de = change of ln V
positions *= mu;  box *= mu
velocities /= mu                       # C-rescale only

β is the (estimated) isothermal compressibility; for water use 4.5 × 10⁻⁵ bar⁻¹. A wrong β only changes the effective coupling speed.
Typical τ_P: 1–2 ps for Berendsen and C-rescale, 2–5 ps for Parrinello–Rahman. τ_P should be several times τ_T.
Discard the part of the run where the volume is still drifting. Then average.
GROMACS: pcoupl = C-rescale or Parrinello-Rahman, pcoupltype = isotropic / semiisotropic / anisotropic, tau-p, compressibility, ref-p. LAMMPS: fix npt temp ... iso / aniso / tri ... (MTTK) or fix press/berendsen.