Introduction

Every atom in your body obeys Newton's second law: force equals mass times acceleration. If you know the forces on every atom at one moment, you can predict where they will be a tiny time-step later. Repeat that calculation billions of times and you have molecular dynamics — a computer movie of molecules folding, reacting, and colliding.

But Newton's equations are continuous, and computers are discrete. Every time-step introduces a small approximation error. The choice of integrator — the algorithm that advances positions and velocities — decides whether those errors stay bounded or spiral out of control.

The naive choice, Euler's method, is simple to write but fatally flawed: it systematically injects energy into the system. Particles speed up over time, the simulation heats up, and eventually everything flies apart. Loup Verlet solved this in 1967 with an integrator that is barely more complex but conserves a shadow energy near the true energy for as long as you care to run. Understanding why is one of the cleanest lessons in numerical analysis.

Try It: Orbit Under Two Integrators

Below, a particle orbits a fixed center under an inverse-square force (like gravity or Coulomb attraction). Both integrators start with identical conditions; only the update rule differs.

<div class="controls">
  <label>Time step (dt): <span id="dtVal">0.04</span>
    <input type="range" id="dtSlider" min="1" max="8" value="4" step="1">
  </label>
  <div class="btns">
    <button id="startBtn">Run</button>
    <button id="resetBtn" class="ghost">Reset</button>
  </div>
</div>
<div class="panels">
  <div class="orbit-panel">
    <div class="panel-label">Orbit path</div>
    <canvas id="orbitCanvas" width="260" height="260"></canvas>
  </div>
  <div class="energy-panel">
    <div class="panel-label">Energy over time</div>
    <canvas id="energyCanvas" width="260" height="130"></canvas>
    <div class="legend">
      <span class="dot euler"></span>Euler &nbsp;
      <span class="dot verlet"></span>Verlet
    </div>
  </div>
</div>
<div class="status" id="status">Press Run to start the simulation.</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; padding: 4px; }
.controls { display: flex; flex-wrap: wrap; align-items: center; gap: .6rem; margin-bottom: .6rem; }
.controls label { font-size: .85rem; display: flex; align-items: center; gap: .4rem; }
input[type=range] { width: 110px; accent-color: #1d3557; }
.btns { display: flex; gap: .4rem; }
button { font: 600 13px system-ui; padding: .35rem .8rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
.panels { display: flex; flex-wrap: wrap; gap: .8rem; }
.orbit-panel, .energy-panel { display: flex; flex-direction: column; gap: .2rem; }
.panel-label { font-size: .75rem; font-weight: 600; color: #555; text-transform: uppercase; letter-spacing: .04em; }
canvas { border: 1px solid #d0d7de; border-radius: 8px; background: #f8fafc; display: block; }
.legend { font-size: .8rem; margin-top: .3rem; display: flex; align-items: center; gap: .2rem; }
.dot { display: inline-block; width: 10px; height: 10px; border-radius: 50%; }
.dot.euler { background: #e63946; }
.dot.verlet { background: #1d7f47; }
.status { font-size: .85rem; font-weight: 600; margin-top: .4rem; min-height: 1.2em; color: #444; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Orbit demo: compare Euler and Verlet (position form) for a 2-D inverse-square
// central force.  F = -GM/r^2  directed toward origin;  GM = 1.
// Euler:  x += v*dt,  v += a*dt            (not symplectic → energy drifts up)
// Verlet: x_new = 2*x - x_prev + a*dt^2   (symplectic → shadow energy conserved)

const GM = 1.0;
const DT_TABLE = [0.005, 0.01, 0.02, 0.04, 0.06, 0.09, 0.13, 0.18];

// Shared initial conditions (circular orbit at radius r0).
const R0 = 0.5;
const V_CIRC = Math.sqrt(GM / R0); // speed for a perfect circular orbit

function accel(x, y) {
  const r2 = x * x + y * y;
  const r  = Math.sqrt(r2);
  const f  = -GM / r2;
  return [f * x / r, f * y / r];
}

function energy(x, y, vx, vy) {
  const r = Math.sqrt(x * x + y * y);
  return 0.5 * (vx * vx + vy * vy) - GM / r;  // KE + PE
}

// State
let running = false;
let animId = null;
let dt = 0.04;
let step = 0;

// Euler state
let ex, ey, evx, evy;
// Verlet state (needs prev + current position; velocity derived for energy)
let vx_pos, vy_pos, vx_prev, vy_prev, vvx, vvy;

const HIST_MAX = 400;  // energy history length
let eulerEnergy  = [];
let verletEnergy = [];
let eulerTrail  = [];
let verletTrail = [];
const TRAIL_MAX = 800;

function initState() {
  // Euler: start at (R0, 0), velocity = (0, V_CIRC)
  ex = R0; ey = 0; evx = 0; evy = V_CIRC;

  // Verlet: bootstrap previous position using one backward Euler step
  vx_pos = R0; vy_pos = 0;
  vvx = 0; vvy = V_CIRC;  // store velocity for energy calc
  const [ax0, ay0] = accel(vx_pos, vy_pos);
  // x_prev = x - v*dt + 0.5*a*dt^2  (from Taylor expansion)
  vx_prev = vx_pos - vvx * dt + 0.5 * ax0 * dt * dt;
  vy_prev = vy_pos - vvy * dt + 0.5 * ay0 * dt * dt;

  eulerEnergy  = [energy(ex, ey, evx, evy)];
  verletEnergy = [energy(vx_pos, vy_pos, vvx, vvy)];
  eulerTrail   = [[ex, ey]];
  verletTrail  = [[vx_pos, vy_pos]];
  step = 0;
}

function stepEuler() {
  const [ax, ay] = accel(ex, ey);
  ex  += evx * dt;
  ey  += evy * dt;
  evx += ax  * dt;
  evy += ay  * dt;
}

function stepVerlet() {
  const [ax, ay] = accel(vx_pos, vy_pos);
  const nx = 2 * vx_pos - vx_prev + ax * dt * dt;
  const ny = 2 * vy_pos - vy_prev + ay * dt * dt;
  // Approximate velocity for energy (central difference)
  vvx = (nx - vx_prev) / (2 * dt);
  vvy = (ny - vy_prev) / (2 * dt);
  vx_prev = vx_pos; vy_prev = vy_pos;
  vx_pos  = nx;     vy_pos  = ny;
}

// Canvas helpers
const orbitC  = document.getElementById('orbitCanvas');
const energyC = document.getElementById('energyCanvas');
const octx    = orbitC.getContext('2d');
const ectx    = energyC.getContext('2d');
const OW = orbitC.width, OH = orbitC.height;
const EW = energyC.width, EH = energyC.height;
const SCALE = 110; // world units → pixels; center = (OW/2, OH/2)

function worldToCanvas(x, y) {
  return [OW / 2 + x * SCALE, OH / 2 - y * SCALE];
}

function drawOrbit() {
  octx.clearRect(0, 0, OW, OH);
  // Center star
  octx.beginPath();
  octx.arc(OW / 2, OH / 2, 5, 0, Math.PI * 2);
  octx.fillStyle = '#f4a261';
  octx.fill();

  function drawTrail(trail, color) {
    if (trail.length < 2) return;
    octx.beginPath();
    const [x0, y0] = worldToCanvas(trail[0][0], trail[0][1]);
    octx.moveTo(x0, y0);
    for (let i = 1; i < trail.length; i++) {
      const [xi, yi] = worldToCanvas(trail[i][0], trail[i][1]);
      octx.lineTo(xi, yi);
    }
    octx.strokeStyle = color;
    octx.lineWidth = 1.5;
    octx.globalAlpha = 0.7;
    octx.stroke();
    octx.globalAlpha = 1;
  }

  drawTrail(eulerTrail, '#e63946');
  drawTrail(verletTrail, '#1d7f47');

  // Current positions
  function drawDot(trail, color) {
    if (!trail.length) return;
    const [x, y] = worldToCanvas(trail[trail.length - 1][0], trail[trail.length - 1][1]);
    octx.beginPath();
    octx.arc(x, y, 5, 0, Math.PI * 2);
    octx.fillStyle = color;
    octx.fill();
  }
  drawDot(eulerTrail, '#e63946');
  drawDot(verletTrail, '#1d7f47');
}

function drawEnergy() {
  ectx.clearRect(0, 0, EW, EH);
  const all = [...eulerEnergy, ...verletEnergy];
  const eMin = Math.min(...all) - 0.05;
  const eMax = Math.max(...all, eMin + 0.01) + 0.05;
  const range = eMax - eMin;

  function toECanvas(idx, e) {
    const n = Math.max(eulerEnergy.length, verletEnergy.length);
    const x = (idx / (HIST_MAX - 1)) * (EW - 20) + 10;
    const y = EH - 10 - ((e - eMin) / range) * (EH - 20);
    return [x, y];
  }

  function drawLine(arr, color) {
    if (arr.length < 2) return;
    ectx.beginPath();
    const [x0, y0] = toECanvas(0, arr[0]);
    ectx.moveTo(x0, y0);
    for (let i = 1; i < arr.length; i++) {
      const [xi, yi] = toECanvas(i, arr[i]);
      ectx.lineTo(xi, yi);
    }
    ectx.strokeStyle = color;
    ectx.lineWidth = 1.5;
    ectx.stroke();
  }

  // Zero-energy line (approximate equilibrium)
  const e0 = eulerEnergy[0];
  const [, yz] = toECanvas(0, e0);
  ectx.beginPath();
  ectx.moveTo(10, yz); ectx.lineTo(EW - 10, yz);
  ectx.strokeStyle = '#bbb';
  ectx.lineWidth = 0.8;
  ectx.setLineDash([4, 3]);
  ectx.stroke();
  ectx.setLineDash([]);

  drawLine(eulerEnergy, '#e63946');
  drawLine(verletEnergy, '#1d7f47');
}

function tick() {
  const STEPS_PER_FRAME = 6;
  for (let i = 0; i < STEPS_PER_FRAME; i++) {
    stepEuler();
    stepVerlet();
    step++;

    eulerTrail.push([ex, ey]);
    if (eulerTrail.length > TRAIL_MAX) eulerTrail.shift();
    verletTrail.push([vx_pos, vy_pos]);
    if (verletTrail.length > TRAIL_MAX) verletTrail.shift();

    if (step % 3 === 0) {
      eulerEnergy.push(energy(ex, ey, evx, evy));
      if (eulerEnergy.length > HIST_MAX) eulerEnergy.shift();
      verletEnergy.push(energy(vx_pos, vy_pos, vvx, vvy));
      if (verletEnergy.length > HIST_MAX) verletEnergy.shift();
    }
  }

  drawOrbit();
  drawEnergy();

  // Report last energies
  const ee = eulerEnergy[eulerEnergy.length - 1];
  const ve = verletEnergy[verletEnergy.length - 1];
  const e0 = eulerEnergy[0] || 1;
  const drift = ((ee - e0) / Math.abs(e0) * 100).toFixed(1);
  document.getElementById('status').textContent =
    fmt('Step %step% | Euler E: %ee% (drift %drift%%) | Verlet E: %ve%', { step: step, ee: ee.toFixed(4), drift: drift, ve: ve.toFixed(4) });

  if (running) animId = requestAnimationFrame(tick);
}

function reset() {
  running = false;
  if (animId) cancelAnimationFrame(animId);
  animId = null;
  dt = DT_TABLE[parseInt(document.getElementById('dtSlider').value) - 1];
  initState();
  drawOrbit();
  drawEnergy();
  document.getElementById('status').textContent = 'Press Run to start the simulation.';
  document.getElementById('startBtn').textContent = 'Run';
}

document.getElementById('dtSlider').addEventListener('input', function () {
  dt = DT_TABLE[parseInt(this.value) - 1];
  document.getElementById('dtVal').textContent = dt.toFixed(3);
  if (!running) {
    initState();
    drawOrbit();
    drawEnergy();
  }
});

document.getElementById('startBtn').addEventListener('click', function () {
  if (running) {
    running = false;
    if (animId) cancelAnimationFrame(animId);
    animId = null;
    this.textContent = 'Run';
  } else {
    running = true;
    this.textContent = 'Pause';
    tick();
  }
});

document.getElementById('resetBtn').addEventListener('click', reset);

// Init
reset();
Result
Console

Watch the energy panel. Euler's energy drifts upward monotonically — the orbit spirals outward and eventually the particle escapes. Verlet's energy oscillates around its starting value and stays bounded no matter how many steps you run. The step size slider lets you see that with a very small step Euler improves, but Verlet remains stable even with coarser steps.

The Real Complexity

Why does Verlet win? The answer lives in the geometry of physics, not just numerics.

  • Euler's method updates position with the current velocity, then updates velocity with the current acceleration: x←x+v⋅dtx ← x + v \cdot dt, v←v+a⋅dtv ← v + a \cdot dt. Each step draws a straight tangent line and overshoots the true curve, always in the direction that adds energy. Error accumulates one-way.
  • Verlet's update computes the new position from the last two positions: x(t+dt) = 2·x(t) − x(t−dt) + a(t)·$dt^{2}$. Velocity appears only implicitly. This symmetric stencil cancels the leading error term, giving second-order accuracy — and more importantly, the map from old state to new state has a special property: it is area-preserving in phase space (position × momentum). Mathematicians call this symplectic.
  • Symplectic integrators exactly preserve a shadow Hamiltonian — a slightly different energy function that is close to the true energy. Because a shadow energy is exactly conserved, the simulation never systematically drifts. Energy oscillates around the true value but does not grow without bound, even after a billion steps.
  • Velocity Verlet (the form used in modern MD codes) makes the velocity explicit: v(t+dt/2) = v(t) + a(t)·dt/2, then x(t+dt) = x(t) + v(t+dt/2)·dt, then v(t+dt) = v(t+dt/2) + a(t+dt)·dt/2. It is mathematically equivalent to Verlet but numerically more convenient — this is what GROMACS, LAMMPS, and NAMD all use.
  • Higher-order methods like Runge-Kutta 4 are not symplectic and show energy drift even though they are more accurate per step. For long simulations the lower-order but symplectic Verlet beats RK4 every time.

The lesson connects to non-convex optimization: the shape of the mathematical landscape — here, the symplectic geometry of phase space — governs which algorithms succeed over long runs, independently of per-step accuracy.

Where It Matters

The Verlet family is not an academic curiosity — it underlies every major molecular simulation today:

  • Drug discovery: simulations of protein-ligand binding identify drug candidates before synthesis. Billions of time-steps demand an integrator that won't drift. GROMACS, AMBER, and NAMD all default to Velocity Verlet or a closely related leapfrog variant.
  • Protein folding: D.E. Shaw Research's Anton supercomputer runs millisecond-scale MD simulations — impossible without a symplectic integrator that holds energy stable across ~101210^{12} steps.
  • Materials science: simulating crack propagation, radiation damage in reactor steel, or battery electrode degradation requires simulations running for nanoseconds, trillions of steps, at temperatures where energy drift would melt the lattice.
  • Quantum chemistry extensions: Car-Parrinello MD couples classical nuclei with quantum electrons; the integrator must remain symplectic to keep both subsystems stable simultaneously.
  • N-body gravity: the same Verlet logic appears in astrophysical simulations of star clusters and galaxy formation, where energy drift over millions of years would produce unphysical mergers.

See also how numerical simulation chooses integrators for differential equations in general, and how protein folding turns these simulations into predictions about life itself.

Conclusion

Loup Verlet's 1967 insight was not that he found a more accurate integrator — Runge-Kutta 4 has a smaller per-step error. His insight was that he found an integrator that respects the geometry of Hamiltonian mechanics. By preserving area in phase space, Verlet's method exactly conserves a shadow energy, and that is enough to keep a simulation physically meaningful for any number of steps.

The next time you see a molecule fold on screen, or a simulation of a drug binding to a protein receptor, remember: behind every frame is a stripped-down formula, x(t+dt) = 2·x(t) − x(t−dt) + a(t)·$dt^{2}$, marching billions of times without letting energy escape. Simple algebra, deep geometry, and the patience to step through time one tiny interval at a time.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

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