Lesson 3.3 · 13 min

The Berendsen thermostat

The simplest way to control temperature is to rescale all velocities. Berendsen does it gently, with a coupling time τ. It reaches the target smoothly, but it squeezes the natural fluctuations, so it is a tool for equilibration, not for production.

Velocity rescaling, and why rescaling exactly every step is too harsh
The Berendsen thermostat and its coupling time τ
Why Berendsen gives the wrong fluctuations, and how stochastic rescaling (V-rescale) fixes it

The idea

In NVT we want the temperature to be T₀ on average. The temperature is set by the kinetic energy (Lesson 3.1), and multiplying every velocity by λ multiplies T by λ². So the simplest thermostat is: measure T, then scale all velocities.

Exact rescaling picks λ = √(T₀/T) at every step. T is then exactly T₀ at all times. That sounds ideal but is wrong: a real system in contact with a heat bath has a fluctuating kinetic energy (spread √(2/N_f)·T). Exact rescaling removes the fluctuations entirely and disturbs the dynamics.

Berendsen (1984) scales more gently, so that T relaxes towards T₀ exponentially with a time constant τ:

dT/dt = (T₀ − T)/τ.

A small τ means strong coupling (fast, close to exact rescaling); a large τ means weak coupling (close to NVE). Typical values: τ = 0.1–1 ps. It is smooth, robust, and great for bringing a system to temperature.

The catch. Berendsen does not sample the canonical ensemble. The kinetic energy fluctuations come out too small, and averages that depend on them (heat capacity, some free energies) are wrong. In long runs it can also slowly pump energy from fast motions into slow ones, such as the overall translation or rotation of a molecule: the "flying ice cube".

The fix: stochastic velocity rescaling (Bussi, Donadio and Parrinello, 2007). Like Berendsen, but with a carefully chosen random term added to the kinetic energy, so it produces exactly the canonical distribution. It keeps Berendsen's simplicity. GROMACS calls it V-rescale; LAMMPS has fix temp/csvr.

Key idea. Berendsen relaxes T to the target with time constant τ. Use it to equilibrate. For production runs, use V-rescale (Bussi) or Nosé–Hoover (Lesson 3.4), which give correct fluctuations.

The algorithm

Every step, after the velocity Verlet update:

K  = 0.5 * sum(m * v**2)
K0 = 0.5 * Nf * kB * T0                 # target kinetic energy

# exact rescaling (too harsh)
lam = sqrt(K0 / K)

# Berendsen
lam = sqrt(1 + dt / tau * (K0 / K - 1))

# Bussi (V-rescale): Berendsen-like, plus a random term
c   = exp(-dt / tau)
R1  = gaussian();  S = chi_square(Nf - 1)
lam = sqrt(c + (1 - c) * K0 / (Nf * K) * (R1**2 + S)
           + 2 * R1 * sqrt(c * (1 - c) * K0 / (Nf * K)))

v *= lam

τ must be much larger than Δt (τ/Δt ≳ 50). Otherwise Berendsen is just exact rescaling.
The measured relaxation is a bit slower than τ, because some of the added energy flows into potential energy.
Thermostat groups separately when they exchange heat slowly (e.g. solute and solvent), and remove centre-of-mass motion to avoid a flying ice cube.
GROMACS: tcoupl = berendsen or v-rescale, with tau_t and ref_t. LAMMPS: fix temp/berendsen or fix temp/csvr, plus fix nve.