MD Thermostats

When chemists and materials scientists want to understand how molecules move, fold or react, they run a molecular dynamics (MD) simulation: Newton's equations of motion, integrated step by step for millions of atoms. The positions and velocities evolve, and statistics collected along the way tell you thermodynamic properties — how stable a protein fold is, how quickly an ion diffuses, what a drug binding event looks like at the atomic scale.

But there is a catch. The real world does not simulate in a vacuum. Experiments are usually done at a fixed temperature — a water bath, a furnace, room temperature in a lab. In a bare simulation, energy conservation means temperature can wander: numerical rounding, initialization errors and small imbalances let the kinetic energy creep up or collapse. The simulation drifts away from the physical ensemble you actually care about.

The solution is a thermostat — an algorithmic device that nudges the particle velocities so the system's average kinetic energy stays pinned to a target temperature $T_0$. The key relationship is the equipartition theorem:

$$\langle E_k \rangle = \frac{3}{2} N k_B T$$

where $N$ is the number of particles and $k_B$ is Boltzmann's constant. Two thermostats dominate in practice: the Berendsen thermostat (1984, simple and stable) and the Nosé–Hoover thermostat (1984–1985, physically rigorous). They are not equally good — but understanding the difference requires knowing what "correct" even means for a simulation ensemble.

The canvas below simulates a small box of particles. Without a thermostat the temperature drifts freely. Choose Berendsen or Nosé–Hoover, set your target temperature, and watch the measured temperature lock onto $T_0$.

<!-- {{c_html_intro}} -->
<div class="controls">
  <label>{{lbl_mode}}
    <select id="mode">
      <option value="none">{{opt_none}}</option>
      <option value="berendsen">Berendsen</option>
      <option value="nh">{{opt_nh}}</option>
    </select>
  </label>
  <label>{{lbl_target}} <b id="t0disp">1.0</b>
    <input type="range" id="t0slider" min="0.3" max="2.5" step="0.05" value="1.0">
  </label>
  <button id="resetBtn" type="button">{{btn_reset}}</button>
</div>
<canvas id="cv" width="480" height="220" aria-label="{{cv_aria}}"></canvas>
<div class="chart-wrap">
  <canvas id="chart" width="480" height="140" aria-label="{{chart_aria}}"></canvas>
</div>
<div id="status" class="status"></div>

/* {{c_css_intro}} */
* { box-sizing: border-box; margin: 0; padding: 0; }
body { font-family: system-ui, sans-serif; font-size: 14px; color: #222; background: #fff; }
.controls { display: flex; flex-wrap: wrap; gap: .6rem .9rem; padding: .5rem 0 .6rem; align-items: center; }
label { display: flex; align-items: center; gap: .3rem; font-size: .85rem; }
select, input[type=range] { cursor: pointer; }
select { padding: .2rem .4rem; border: 1px solid #adb1b8; border-radius: 6px; background: #f5f6f8; }
input[type=range] { width: 90px; }
button { font: 600 13px system-ui; padding: .3rem .7rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 6px; cursor: pointer; }
canvas { display: block; border: 1px solid #dde2e8; border-radius: 6px; background: #f8f9fb; }
.chart-wrap { margin-top: .45rem; }
.status { margin-top: .4rem; font-size: .9rem; font-weight: 600; min-height: 1.4em; }
.status.ok { color: #0a7d33; }
.status.drift { color: #c92f3c; }
.status.settling { color: #e07e00; }

// Code not found

Notice how Berendsen converges smoothly but too tightly — the temperature barely fluctuates. Nosé–Hoover takes a little longer to settle and then oscillates around $T_0$. Those fluctuations are not noise: they are thermodynamically correct. Real systems at fixed temperature do fluctuate; only the average is pinned. The Berendsen thermostat suppresses those fluctuations and therefore does not sample the true NVT (constant $N$, $V$, $T$) ensemble — a subtle but important flaw for computing free energies or heat capacities.

Both thermostats share the goal of keeping the instantaneous temperature

$$T(t) = \frac{2 E_k(t)}{3 N k_B}$$

near $T_0$, but they take very different routes.

Berendsen (1984). At each time step, multiply every velocity by a rescaling factor:

$$\lambda = \sqrt{1 + \frac{\Delta t}{\tau} \left(\frac{T_0}{T(t)} - 1\right)}$$

where $\tau$ is a coupling time constant. When $T > T_0$ the factor shrinks velocities; when $T < T_0$ it inflates them. The result is exponential relaxation toward $T_0$ with time constant $\tau$. It is robust, easy to implement, and almost never crashes — which is why it is used routinely for equilibration. The problem is that it clamps fluctuations: the system is forced toward exactly $T_0$ rather than being allowed to fluctuate around it. This means it does not sample the canonical (NVT) ensemble correctly, making it unreliable for equilibrium thermodynamic averages.

Nosé–Hoover (1984–1985). Shuichi Nosé introduced a fictional extra degree of freedom $s$ — a "heat bath variable" — with its own mass $Q$ and conjugate momentum. The extended Lagrangian generates equations of motion that include a friction-like term $\xi$ for the velocities:

$$\dot{\xi} = \frac{1}{Q}\left(\sum_i m_i v_i^2 - 3 N k_B T_0\right)$$

When the kinetic energy is too high, $\xi$ grows and slows the particles; when too low, $\xi$ shrinks and speeds them up. Crucially, the time-averaged distribution of positions and momenta is exactly the canonical ensemble. The "mass" $Q$ controls how quickly the thermostat responds: too small and it oscillates wildly; too large and it barely couples. Getting $Q$ right — usually set to $3 N k_B T_0 \tau^2$ — is the main tuning challenge.

The lesson from both: thermostats are not a free lunch. Controlling temperature requires injecting or removing energy from the system, which changes its dynamics. The art is doing so in a way that preserves the statistical mechanics you care about.

Thermostat choice is not an academic detail — it shapes the physical validity of every result:

• Drug discovery and protein folding: free energy calculations (binding affinities, folding stabilities) require the canonical ensemble. Berendsen gives wrong free energies; Nosé–Hoover or similar methods are mandatory. See protein folding for the broader complexity of predicting structure.
• Materials design: simulating how a polymer crystallizes, how a metal deforms, or how a battery electrode swells depends on sampling at the right temperature. Wrong ensemble → wrong phase diagram.
• Chemical reaction rates: transition-state theory relies on the Boltzmann distribution over energies. A Berendsen-equilibrated system has the right mean energy but wrong fluctuations, leading to rate errors.
• Climate and atmospheric chemistry: coarse-grained MD-style models of atmospheric aerosols or ocean water use similar thermostat ideas to represent heat exchange with the environment.
• Teaching statistical mechanics: the contrast between Berendsen and Nosé–Hoover is one of the clearest demonstrations that getting the right average is not enough — the fluctuations carry real physics too.

The topic also connects directly to quantum simulation: at low temperatures, quantum nuclear effects matter and classical thermostats must be extended to path-integral MD with quantum heat baths.

Keeping simulated atoms at a target temperature sounds like a plumbing problem — a dial you set and forget. But the two main solutions, proposed within months of each other in 1984, turned out to encode a deep lesson about what it means to simulate thermodynamics correctly.

Berendsen's rescaling is fast and forgiving: perfect for equilibrating a system before the real measurement begins, disastrous if you use it during the measurement. Nosé–Hoover is slower to converge but generates the true canonical ensemble — the one where fluctuations carry real physical meaning.

Forty years later, every major MD package (GROMACS, NAMD, LAMMPS, OpenMM) ships both. The right choice still depends on what you are measuring. That is not a failure of the field — it is the field doing its job: knowing which problem you are actually solving, and picking the tool that matches.