MD Thermostats: Interactive Particle Demo

Introduction

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.

Try It

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$ .

<!-- Controls: thermostat mode selector, target temperature slider, reset button -->
<div class="controls">
  <label>Thermostat:
    <select id="mode">
      <option value="none">None (free drift)</option>
      <option value="berendsen">Berendsen</option>
      <option value="nh">Nosé–Hoover</option>
    </select>
  </label>
  <label>Target T₀: <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">Reset</button>
</div>
<canvas id="cv" width="480" height="220" aria-label="Particle simulation canvas"></canvas>
<div class="chart-wrap">
  <canvas id="chart" width="480" height="140" aria-label="Temperature history chart"></canvas>
</div>
<div id="status" class="status"></div>

/* Layout and styling for the MD thermostat demo */
* { 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; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
/* MD thermostat demo — Berendsen velocity rescaling vs Nosé–Hoover heat bath */

/* Simulation parameters */
const N_PART = 40;
const BOX_W = 480, BOX_H = 220;
const RADIUS = 5;
const DT = 0.018;
const TAU_BEREN = 0.5;
const TAU_NH = 0.4;
const HISTORY_LEN = 200;

/* Simulation state variables */
let pos, vel, xi, T0, mode, tHistory, frameId;

function randGauss() {
  /* Box–Muller transform for Gaussian random numbers */
  let u = 0, v = 0;
  while (u === 0) u = Math.random();
  while (v === 0) v = Math.random();
  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
}

function initSystem() {
  /* Initialize particle positions and velocities */
  pos = [];
  vel = [];
  xi = 0;
  T0 = parseFloat(document.getElementById('t0slider').value);
  mode = document.getElementById('mode').value;
  tHistory = [];

  /* Arrange particles on a regular grid with small random offsets */
  const cols = Math.ceil(Math.sqrt(N_PART));
  const rows = Math.ceil(N_PART / cols);
  const dx = (BOX_W - 2 * RADIUS * 2) / cols;
  const dy = (BOX_H - 2 * RADIUS * 2) / rows;
  let idx = 0;
  for (let r = 0; r < rows && idx < N_PART; r++) {
    for (let c = 0; c < cols && idx < N_PART; c++) {
      pos.push([
        RADIUS * 2 + (c + 0.5) * dx + (Math.random() - 0.5) * dx * 0.3,
        RADIUS * 2 + (r + 0.5) * dy + (Math.random() - 0.5) * dy * 0.3
      ]);
      /* Maxwell–Boltzmann initial velocities */
      vel.push([randGauss() * Math.sqrt(T0 * 0.35), randGauss() * Math.sqrt(T0 * 0.35)]);
      idx++;
    }
  }
  /* Remove center-of-mass drift */
  let vxMean = 0, vyMean = 0;
  vel.forEach(v => { vxMean += v[0]; vyMean += v[1]; });
  vxMean /= N_PART; vyMean /= N_PART;
  vel.forEach(v => { v[0] -= vxMean; v[1] -= vyMean; });
}

function kineticTemp() {
  /* Instantaneous kinetic temperature: T = <v^2> / 2 (2D, kB=m=1) */
  let ek = 0;
  vel.forEach(v => { ek += v[0] * v[0] + v[1] * v[1]; });
  return ek / (2 * N_PART);
}

function step() {
  /* Euler integration: move particles */
  for (let i = 0; i < N_PART; i++) {
    pos[i][0] += vel[i][0] * DT;
    pos[i][1] += vel[i][1] * DT;
    /* Reflect off walls (elastic bounce) */
    if (pos[i][0] < RADIUS) { pos[i][0] = RADIUS; vel[i][0] = Math.abs(vel[i][0]); }
    if (pos[i][0] > BOX_W - RADIUS) { pos[i][0] = BOX_W - RADIUS; vel[i][0] = -Math.abs(vel[i][0]); }
    if (pos[i][1] < RADIUS) { pos[i][1] = RADIUS; vel[i][1] = Math.abs(vel[i][1]); }
    if (pos[i][1] > BOX_H - RADIUS) { pos[i][1] = BOX_H - RADIUS; vel[i][1] = -Math.abs(vel[i][1]); }
  }
  const T = kineticTemp();
  if (mode === 'berendsen') {
    /* Berendsen velocity rescaling: lambda couples T toward T0 */
    const lam = Math.sqrt(1 + (DT / TAU_BEREN) * (T0 / Math.max(T, 0.001) - 1));
    vel.forEach(v => { v[0] *= lam; v[1] *= lam; });
  } else if (mode === 'nh') {
    /* Nosé–Hoover friction: xi drives kinetic energy toward N*kB*T0 */
    const Q = 3 * N_PART * T0 * TAU_NH * TAU_NH;
    xi += DT * (2 * N_PART * T - 2 * N_PART * T0) / Q;
    vel.forEach(v => { v[0] -= xi * v[0] * DT; v[1] -= xi * v[1] * DT; });
  }
  /* Record temperature history for the chart */
  tHistory.push(kineticTemp());
  if (tHistory.length > HISTORY_LEN) tHistory.shift();
}

function drawParticles() {
  const cv = document.getElementById('cv');
  const ctx = cv.getContext('2d');
  ctx.clearRect(0, 0, BOX_W, BOX_H);
  const T = kineticTemp();
  /* Color particles by speed: blue = slow, orange = fast */
  pos.forEach((p, i) => {
    const speed = Math.sqrt(vel[i][0] * vel[i][0] + vel[i][1] * vel[i][1]);
    const hot = Math.min(speed / (3 * Math.sqrt(T0 + 0.1)), 1);
    ctx.beginPath();
    ctx.arc(p[0], p[1], RADIUS, 0, 2 * Math.PI);
    ctx.fillStyle = `hsl(${220 - hot * 200},70%,${55 - hot * 15}%)`;
    ctx.fill();
  });
}

function drawChart() {
  const cv = document.getElementById('chart');
  const ctx = cv.getContext('2d');
  const W = cv.width, H = cv.height;
  ctx.clearRect(0, 0, W, H);
  /* Draw temperature chart with axes and gridlines */
  const pad = { l: 38, r: 8, t: 8, b: 22 };
  const iW = W - pad.l - pad.r;
  const iH = H - pad.t - pad.b;
  const maxT = Math.max(T0 * 2.2, ...tHistory, 0.1);
  ctx.strokeStyle = '#ccd2da';
  ctx.lineWidth = 1;
  /* Horizontal gridlines at fractions of T0 */
  [0.25, 0.5, 0.75, 1.0].forEach(f => {
    const y = pad.t + iH * (1 - f * T0 * 2 / maxT);
    if (y < pad.t || y > pad.t + iH) return;
    ctx.beginPath(); ctx.moveTo(pad.l, y); ctx.lineTo(pad.l + iW, y); ctx.stroke();
    ctx.fillStyle = '#888'; ctx.font = '10px system-ui';
    ctx.fillText((f * T0 * 2).toFixed(2), 2, y + 3);
  });
  /* Red dashed line at target temperature T0 */
  ctx.strokeStyle = '#e63946';
  ctx.lineWidth = 1.2;
  ctx.setLineDash([4, 3]);
  const ty = pad.t + iH * (1 - T0 / maxT);
  ctx.beginPath(); ctx.moveTo(pad.l, ty); ctx.lineTo(pad.l + iW, ty); ctx.stroke();
  ctx.setLineDash([]);
  ctx.fillStyle = '#e63946'; ctx.font = '10px system-ui';
  ctx.fillText('T₀', pad.l + iW - 14, ty - 3);
  /* Blue line: measured temperature over time */
  if (tHistory.length > 1) {
    ctx.strokeStyle = '#1d6aa5';
    ctx.lineWidth = 1.6;
    ctx.beginPath();
    tHistory.forEach((t, i) => {
      const x = pad.l + (i / (HISTORY_LEN - 1)) * iW;
      const y = pad.t + iH * (1 - t / maxT);
      i === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
    });
    ctx.stroke();
  }
  /* Axis labels */
  ctx.fillStyle = '#444'; ctx.font = '10px system-ui';
  ctx.fillText('T', 4, pad.t + 4);
  ctx.fillText('time →', pad.l + iW / 2 - 12, H - 4);
}

function updateStatus() {
  const el = document.getElementById('status');
  const T = kineticTemp();
  const diff = Math.abs(T - T0) / T0;
  if (mode === 'none') {
    el.textContent = i18nFmt('No thermostat — T = %T% (drifting)', { T: T.toFixed(3) });
    el.className = 'status drift';
  } else if (diff < 0.05) {
    el.textContent = i18nFmt('Stable at T = %T% (target %T0%)', { T: T.toFixed(3), T0: T0.toFixed(2) });
    el.className = 'status ok';
  } else {
    el.textContent = i18nFmt('Settling… T = %T% → %T0%', { T: T.toFixed(3), T0: T0.toFixed(2) });
    el.className = 'status settling';
  }
}

function loop() {
  for (let i = 0; i < 3; i++) step();
  drawParticles();
  drawChart();
  updateStatus();
  frameId = requestAnimationFrame(loop);
}

function reset() {
  if (frameId) cancelAnimationFrame(frameId);
  initSystem();
  loop();
}

document.getElementById('resetBtn').addEventListener('click', reset);
document.getElementById('mode').addEventListener('change', () => {
  mode = document.getElementById('mode').value;
  xi = 0;
});
document.getElementById('t0slider').addEventListener('input', () => {
  T0 = parseFloat(document.getElementById('t0slider').value);
  document.getElementById('t0disp').textContent = T0.toFixed(2);
});

reset();

<!-- Controls: thermostat mode selector, target temperature slider, reset button -->
<div class="controls">
  <label>Thermostat:
    <select id="mode">
      <option value="none">None (free drift)</option>
      <option value="berendsen">Berendsen</option>
      <option value="nh">Nosé–Hoover</option>
    </select>
  </label>
  <label>Target T₀: <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">Reset</button>
</div>
<canvas id="cv" width="480" height="220" aria-label="Particle simulation canvas"></canvas>
<div class="chart-wrap">
  <canvas id="chart" width="480" height="140" aria-label="Temperature history chart"></canvas>
</div>
<div id="status" class="status"></div>

/* Layout and styling for the MD thermostat demo */
* { 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; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
/* MD thermostat demo — Berendsen velocity rescaling vs Nosé–Hoover heat bath */

/* Simulation parameters */
const N_PART = 40;
const BOX_W = 480, BOX_H = 220;
const RADIUS = 5;
const DT = 0.018;
const TAU_BEREN = 0.5;
const TAU_NH = 0.4;
const HISTORY_LEN = 200;

/* Simulation state variables */
let pos, vel, xi, T0, mode, tHistory, frameId;

function randGauss() {
  /* Box–Muller transform for Gaussian random numbers */
  let u = 0, v = 0;
  while (u === 0) u = Math.random();
  while (v === 0) v = Math.random();
  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
}

function initSystem() {
  /* Initialize particle positions and velocities */
  pos = [];
  vel = [];
  xi = 0;
  T0 = parseFloat(document.getElementById('t0slider').value);
  mode = document.getElementById('mode').value;
  tHistory = [];

  /* Arrange particles on a regular grid with small random offsets */
  const cols = Math.ceil(Math.sqrt(N_PART));
  const rows = Math.ceil(N_PART / cols);
  const dx = (BOX_W - 2 * RADIUS * 2) / cols;
  const dy = (BOX_H - 2 * RADIUS * 2) / rows;
  let idx = 0;
  for (let r = 0; r < rows && idx < N_PART; r++) {
    for (let c = 0; c < cols && idx < N_PART; c++) {
      pos.push([
        RADIUS * 2 + (c + 0.5) * dx + (Math.random() - 0.5) * dx * 0.3,
        RADIUS * 2 + (r + 0.5) * dy + (Math.random() - 0.5) * dy * 0.3
      ]);
      /* Maxwell–Boltzmann initial velocities */
      vel.push([randGauss() * Math.sqrt(T0 * 0.35), randGauss() * Math.sqrt(T0 * 0.35)]);
      idx++;
    }
  }
  /* Remove center-of-mass drift */
  let vxMean = 0, vyMean = 0;
  vel.forEach(v => { vxMean += v[0]; vyMean += v[1]; });
  vxMean /= N_PART; vyMean /= N_PART;
  vel.forEach(v => { v[0] -= vxMean; v[1] -= vyMean; });
}

function kineticTemp() {
  /* Instantaneous kinetic temperature: T = <v^2> / 2 (2D, kB=m=1) */
  let ek = 0;
  vel.forEach(v => { ek += v[0] * v[0] + v[1] * v[1]; });
  return ek / (2 * N_PART);
}

function step() {
  /* Euler integration: move particles */
  for (let i = 0; i < N_PART; i++) {
    pos[i][0] += vel[i][0] * DT;
    pos[i][1] += vel[i][1] * DT;
    /* Reflect off walls (elastic bounce) */
    if (pos[i][0] < RADIUS) { pos[i][0] = RADIUS; vel[i][0] = Math.abs(vel[i][0]); }
    if (pos[i][0] > BOX_W - RADIUS) { pos[i][0] = BOX_W - RADIUS; vel[i][0] = -Math.abs(vel[i][0]); }
    if (pos[i][1] < RADIUS) { pos[i][1] = RADIUS; vel[i][1] = Math.abs(vel[i][1]); }
    if (pos[i][1] > BOX_H - RADIUS) { pos[i][1] = BOX_H - RADIUS; vel[i][1] = -Math.abs(vel[i][1]); }
  }
  const T = kineticTemp();
  if (mode === 'berendsen') {
    /* Berendsen velocity rescaling: lambda couples T toward T0 */
    const lam = Math.sqrt(1 + (DT / TAU_BEREN) * (T0 / Math.max(T, 0.001) - 1));
    vel.forEach(v => { v[0] *= lam; v[1] *= lam; });
  } else if (mode === 'nh') {
    /* Nosé–Hoover friction: xi drives kinetic energy toward N*kB*T0 */
    const Q = 3 * N_PART * T0 * TAU_NH * TAU_NH;
    xi += DT * (2 * N_PART * T - 2 * N_PART * T0) / Q;
    vel.forEach(v => { v[0] -= xi * v[0] * DT; v[1] -= xi * v[1] * DT; });
  }
  /* Record temperature history for the chart */
  tHistory.push(kineticTemp());
  if (tHistory.length > HISTORY_LEN) tHistory.shift();
}

function drawParticles() {
  const cv = document.getElementById('cv');
  const ctx = cv.getContext('2d');
  ctx.clearRect(0, 0, BOX_W, BOX_H);
  const T = kineticTemp();
  /* Color particles by speed: blue = slow, orange = fast */
  pos.forEach((p, i) => {
    const speed = Math.sqrt(vel[i][0] * vel[i][0] + vel[i][1] * vel[i][1]);
    const hot = Math.min(speed / (3 * Math.sqrt(T0 + 0.1)), 1);
    ctx.beginPath();
    ctx.arc(p[0], p[1], RADIUS, 0, 2 * Math.PI);
    ctx.fillStyle = `hsl(${220 - hot * 200},70%,${55 - hot * 15}%)`;
    ctx.fill();
  });
}

function drawChart() {
  const cv = document.getElementById('chart');
  const ctx = cv.getContext('2d');
  const W = cv.width, H = cv.height;
  ctx.clearRect(0, 0, W, H);
  /* Draw temperature chart with axes and gridlines */
  const pad = { l: 38, r: 8, t: 8, b: 22 };
  const iW = W - pad.l - pad.r;
  const iH = H - pad.t - pad.b;
  const maxT = Math.max(T0 * 2.2, ...tHistory, 0.1);
  ctx.strokeStyle = '#ccd2da';
  ctx.lineWidth = 1;
  /* Horizontal gridlines at fractions of T0 */
  [0.25, 0.5, 0.75, 1.0].forEach(f => {
    const y = pad.t + iH * (1 - f * T0 * 2 / maxT);
    if (y < pad.t || y > pad.t + iH) return;
    ctx.beginPath(); ctx.moveTo(pad.l, y); ctx.lineTo(pad.l + iW, y); ctx.stroke();
    ctx.fillStyle = '#888'; ctx.font = '10px system-ui';
    ctx.fillText((f * T0 * 2).toFixed(2), 2, y + 3);
  });
  /* Red dashed line at target temperature T0 */
  ctx.strokeStyle = '#e63946';
  ctx.lineWidth = 1.2;
  ctx.setLineDash([4, 3]);
  const ty = pad.t + iH * (1 - T0 / maxT);
  ctx.beginPath(); ctx.moveTo(pad.l, ty); ctx.lineTo(pad.l + iW, ty); ctx.stroke();
  ctx.setLineDash([]);
  ctx.fillStyle = '#e63946'; ctx.font = '10px system-ui';
  ctx.fillText('T₀', pad.l + iW - 14, ty - 3);
  /* Blue line: measured temperature over time */
  if (tHistory.length > 1) {
    ctx.strokeStyle = '#1d6aa5';
    ctx.lineWidth = 1.6;
    ctx.beginPath();
    tHistory.forEach((t, i) => {
      const x = pad.l + (i / (HISTORY_LEN - 1)) * iW;
      const y = pad.t + iH * (1 - t / maxT);
      i === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
    });
    ctx.stroke();
  }
  /* Axis labels */
  ctx.fillStyle = '#444'; ctx.font = '10px system-ui';
  ctx.fillText('T', 4, pad.t + 4);
  ctx.fillText('time →', pad.l + iW / 2 - 12, H - 4);
}

function updateStatus() {
  const el = document.getElementById('status');
  const T = kineticTemp();
  const diff = Math.abs(T - T0) / T0;
  if (mode === 'none') {
    el.textContent = i18nFmt('No thermostat — T = %T% (drifting)', { T: T.toFixed(3) });
    el.className = 'status drift';
  } else if (diff < 0.05) {
    el.textContent = i18nFmt('Stable at T = %T% (target %T0%)', { T: T.toFixed(3), T0: T0.toFixed(2) });
    el.className = 'status ok';
  } else {
    el.textContent = i18nFmt('Settling… T = %T% → %T0%', { T: T.toFixed(3), T0: T0.toFixed(2) });
    el.className = 'status settling';
  }
}

function loop() {
  for (let i = 0; i < 3; i++) step();
  drawParticles();
  drawChart();
  updateStatus();
  frameId = requestAnimationFrame(loop);
}

function reset() {
  if (frameId) cancelAnimationFrame(frameId);
  initSystem();
  loop();
}

document.getElementById('resetBtn').addEventListener('click', reset);
document.getElementById('mode').addEventListener('change', () => {
  mode = document.getElementById('mode').value;
  xi = 0;
});
document.getElementById('t0slider').addEventListener('input', () => {
  T0 = parseFloat(document.getElementById('t0slider').value);
  document.getElementById('t0disp').textContent = T0.toFixed(2);
});

reset();

Result

Console

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.

The Real Complexity

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.

Where It Matters

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.

Conclusion

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.

Comments

Loading comments...

https://www.kipuhub.com/en/article/molecular-dynamics-thermostats/Content licensed under CC BY-NC 4.0.

Cookies & privacy

MD Thermostats

When atoms won't stay still at the temperature you want