Lesson 3.4 · 15 min
Nosé–Hoover thermostats
Nosé–Hoover adds one extra dynamical variable, a friction ξ, that brakes the atoms when they are too hot and pushes them when too cold. It is deterministic, samples the canonical ensemble, and conserves an extended energy. Chains of thermostats fix its one weakness.
- The Nosé–Hoover equations: a friction ξ driven by the temperature error
- The thermostat mass Q (or period τ), the conserved extended energy, and ringing
- Why one Nosé–Hoover thermostat can fail (ergodicity) and how chains fix it
The idea
Berendsen (Lesson 3.3) reaches the target temperature but gives the wrong fluctuations. Nosé (1984) and Hoover (1985) found a deterministic way to get them right. The heat bath is represented by one extra dynamical variable, a friction coefficient ξ, with its own equation of motion:
- The atoms feel an extra force −ξ·m·v. If ξ > 0 it brakes them; if ξ < 0 it pushes them.
- ξ itself grows when the system is hotter than T₀ and shrinks when it is colder. It is like a "thermal piston" with inertia Q.
Because ξ has inertia, it overshoots: T and ξ oscillate around the target with a period set by Q. Writing Q = NfkBT₀τ² makes τ that period. Typical choices are τ ≈ 0.1–1 ps (about 100 time steps or more). Too small a τ gives fast, violent oscillations; too large makes the coupling slow.
Why it works. One can prove that, if the dynamics is ergodic, the atoms sample exactly the canonical distribution. There is also a conserved quantity, an extended energy that includes the thermostat terms. Its drift is your energy check in NVT, just like E in NVE (Lesson 1.6).
Two weaknesses. (1) A big disturbance makes the thermostat ring for a long time; equilibrate first (Berendsen or V-rescale), then switch. (2) For small or very regular systems (a harmonic oscillator, a stiff crystal), the dynamics is not ergodic: trajectories stay on a narrow band and never fill the canonical distribution. The fix is a Nosé–Hoover chain (Martyna, Klein and Tuckerman, 1992): ξ₁ is itself thermostatted by ξ₂, and so on. A chain of 3–5 is standard.
The algorithm
A symmetric step (the thermostat acts for half a step before and after velocity Verlet):
Q = Nf * kB * T0 * tau**2
def thermostat_half(v, h): # advance xi and scale v over time h
xi += (h/2) * (sum(m v**2) - Nf kB T0) / Q
v *= exp(-xi * h)
eta += xi * h
xi += (h/2) * (sum(m v**2) - Nf kB T0) / Q
every step:
thermostat_half(v, dt/2)
velocity Verlet step
thermostat_half(v, dt/2)
H_ext = K + U + Q xi**2 / 2 + Nf kB T0 eta # should stay constant
- For a chain, each ξj is driven by the kinetic energy of the link below and damped by the link above. The code tab shows the full routine for a chain of any length.
- Do not start Nosé–Hoover far from the target temperature: it rings. Equilibrate with Berendsen or V-rescale first.
- The friction is a smooth, global force, so dynamics (diffusion, spectra) is barely disturbed when τ is large enough.
- LAMMPS:
fix nvt temp T T tau(a chain of 3 by default;tau≈ 100 × timestep). GROMACS:tcoupl = nose-hoover,nsttcouple,nh-chain-length(default 10; the leap-frogmdintegrator supports only 1).