Lab: thermostats by hand

The setup

We reuse the 9 atoms of Lab 2 on the 3×3 lattice in the periodic box L = 3.6, with the same Gaussian velocities (drift removed). This time the velocities are scaled so that the system starts too hot, at T = 2.0, and the thermostat has to bring it down to T₀ = 0.72.

One idea runs through the whole lab. Temperature is an average property; the kinetic energy of a finite system at a fixed temperature is not constant, it fluctuates. A thermostat that removes these fluctuations gives the right average with the wrong physics.

Step 1

How much should the temperature fluctuate?

Why: Before we can judge a thermostat, we need the right answer. A system in contact with a heat bath (the canonical, NVT ensemble) exchanges energy with it all the time, so its kinetic energy goes up and down. How much it goes up and down is fixed by statistics, and we can compute it by hand.

What the numbers say: For our 9 atoms at T₀ = 0.72, the instantaneous temperature should wander by about ±0.25 (35 %) around 0.72. That is not noise to be removed: it is the correct physics. In a large system it is smaller (2.6 % for 1000 atoms) but still there, and quantities like the heat capacity are computed from it. This number, 0.354, is the test every thermostat below must pass.

Step 2

Simple velocity rescaling: right temperature, no fluctuations

Why: The most obvious idea: if T is 2.0 and we want 0.72, multiply every velocity by the same factor λ so that T becomes exactly 0.72. Since T is proportional to v², we need λ² = T₀/T.

What the numbers say: Merit: it reaches the target in one step, which is useful to remove a huge excess of heat at the very start. Limitation: applied every step, the temperature never fluctuates (0.000 instead of 0.354), so the system is not in the canonical ensemble. The sudden velocity jumps also distort the dynamics. Use it only for a first rough cooling, never for collecting data.

Step 3

Berendsen: a gentle version of rescaling

Why: Instead of jumping to T₀ in one step, Berendsen moves T only a small fraction dt/τ of the way each step. τ is a coupling time chosen by the user. This removes the violent jumps, but we must check whether the fluctuations come back.

(These hand steps show only the thermostat part of each step. In a real run the forces also move energy between kinetic and potential energy between thermostat steps; the long runs below include everything.) Now let it run, and watch two numbers in windows of 100 time units (20 000 steps): the kinetic-energy fluctuation, and the share of the kinetic energy in the sliding modes.

Each cell is sliding share / KE fluctuation.

What the numbers say: Merit: smooth, predictable cooling: the deviation from T₀ falls to 1/e in a time τ. That makes Berendsen excellent for equilibration. Limitation 1: even while things look healthy, the fluctuation is 0.15 (τ = 0.1) or 0.17 (τ = 1.0) instead of 0.354, so the ensemble is wrong, only less wrong for weak coupling. Limitation 2, much worse: with τ = 0.1, after about 250 time units the kinetic energy drains into the 4 sliding modes, until 99 % of it is there (fair share 25 %) and the fluctuation has collapsed to 0.002. The atoms inside the rows have stopped vibrating, they are frozen, while whole rows glide past each other: the thermometer still reads 0.72. Why: rescaling all velocities by a common factor every step breaks the balance of energy flow between motions of different speed (Harvey, Tan and Cheatham, 1998). The slow collective motions exchange energy with the rest very slowly, so each small bias adds up there, step after step. In big simulations this is the famous flying ice cube: the whole system moves as a block while it freezes inside. Weak coupling (τ = 1.0) did not show it within 1000 time units, but it still has the wrong fluctuation. This is why Berendsen must not be used for production runs.

Step 4

Nosé–Hoover: a friction that has its own equation

Why: We want the right average and the right fluctuations, and we do not want to force T at every step. Nosé–Hoover adds one extra variable, a friction ξ, which acts on all velocities. ξ grows when the system is too hot and shrinks (even becomes negative, i.e. heats) when it is too cold. Its "mass" Q sets how quickly it responds.

The friction builds up step by step, cooling gets faster, and later ξ overshoots: T drops below 0.72, ξ turns negative and heats the system back. T oscillates around T₀ instead of being forced onto it. The long runs:

What the numbers say: Merit: for the 9 atoms, Nosé–Hoover gets the fluctuation (0.350) and the energy sharing (0.249) right, and it is deterministic, so a run can be reproduced exactly. It is mathematically proven to give the canonical ensemble if the motion explores all states. Limitation: a very simple system, such as one harmonic oscillator (or a stiff, nearly harmonic solid), does not explore all states. Its fluctuation is 0.904 instead of 1.414 and fast velocities are too rare: in the widget you can see unvisited holes in the (x, v) plane. The fix is a chain: a second friction acting on the first one (NH chain, LAMMPS default tchain 3), which gives 1.428 and 0.132. Also, τ must not be too small or T oscillates strongly; 100 time steps is a common choice.

Step 5

Bussi (CSVR): Berendsen plus the missing randomness

Why: Berendsen fails because it pushes K towards its target smoothly, without any random kicks. Bussi, Donadio and Parrinello (2007) keep the same pull and add exactly the random part a heat bath would give, sized so that the fluctuations come out right.

What the numbers say: Merit: with R = 0 it is exactly Berendsen (smooth equilibration), and the random term restores the correct fluctuation (0.350) and energy sharing (0.250). It is robust even for the oscillator-like systems where plain NH fails, and since all velocities are scaled by one common factor, the diffusion (0.0597) is close to the Nosé–Hoover value (0.0575). Limitation: the run depends on the random numbers, so it is reproducible only with the same seed; total momentum is kept, but there is no conserved "energy" in the usual sense (codes report an effective one). In GROMACS it is v-rescale, in LAMMPS fix temp/csvr.

Step 6

Langevin: friction and random kicks on every atom

Why: A real heat bath is made of solvent molecules that slow an atom down (friction γ) and hit it at random. Langevin dynamics adds exactly these two forces to every atom separately. Because the kicks are on each atom, it does not matter whether the motion explores all states.

What the numbers say: Merit: the fluctuations and the energy sharing are right (0.330 and 0.330 vs 0.333), it works for any system, including the single oscillator (0.135 vs 0.134), and it is stable even far from equilibrium. Limitation: the friction slows the atoms down. Diffusion drops from 0.058 to 0.043 at γ = 1 and to 0.015 at γ = 10, four times too slow. Momentum is not conserved, so collective flow is damped too. Use Langevin for sampling structures and energies, or with a small γ; do not take diffusion coefficients or viscosities from a strongly coupled Langevin run.

Step 7

Pressure and barostats: the same story for the volume

Why: To simulate at constant pressure (NPT), we first need the pressure of a configuration, and then a rule that changes the box size. The pressure has a kinetic part (atoms hitting the walls) and a virial part (forces between atoms, Lesson 3.6).

What the numbers say: A negative pressure means the atoms pull on each other more than they push: at this density and temperature the box is too large, so the barostat shrinks it. Merit of Berendsen scaling: smooth and fast, ideal to bring the density to the right value. Limitation: exactly like the Berendsen thermostat, it suppresses the volume fluctuations, so the ensemble is not truly NPT and the compressibility comes out wrong. For production use a barostat with the correct fluctuations: C-rescale (Bernetti–Bussi, the barostat version of step 5), Parrinello–Rahman, or MTK (LAMMPS fix npt).

Summary