Molecular Dynamics Integrators

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.

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>{{lbl_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">{{btn_run}}</button>
    <button id="resetBtn" class="ghost">{{btn_reset}}</button>
  </div>
</div>
<div class="panels">
  <div class="orbit-panel">
    <div class="panel-label">{{lbl_orbit}}</div>
    <canvas id="orbitCanvas" width="260" height="260"></canvas>
  </div>
  <div class="energy-panel">
    <div class="panel-label">{{lbl_energy}}</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">{{status_initial}}</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; }

// Code not found

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.

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 \cdot dt$, $v ← 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.

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 ~$10^{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.

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.