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 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:

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.

Key idea. Nosé–Hoover chains are the standard deterministic thermostat for production NVT. Pick τ ≈ 100 time steps or more, equilibrate before switching on, and monitor the extended energy.

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
  1. 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.
  2. Do not start Nosé–Hoover far from the target temperature: it rings. Equilibrate with Berendsen or V-rescale first.
  3. The friction is a smooth, global force, so dynamics (diffusion, spectra) is barely disturbed when τ is large enough.
  4. 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-frog md integrator supports only 1).