Verlet / leapfrog

Deep dive

How it works

In: positions x, velocities v, a force function a(x), fixed step h · Out: x and v at the next step · State: x, v and the last acceleration (kept so each step needs one force evaluation) · Cost: 1 force evaluation per step; O(h²) global error

Velocity Verlet (kick-drift-kick), with a₀ = a(x₀) computed once at the start:

v_half = v + h/2 · a(x)
x_new  = x + h · v_half
a_new  = a(x_new)              ← the only force evaluation
v_new  = v_half + h/2 · a_new

Leapfrog stores v at half steps (v_{n+½}) and is algebraically the same trajectory. Position (Störmer) Verlet eliminates v entirely: x_{n+1} = 2x_n − x_{n−1} + h²·a(x_n); velocity is implicit in the last two positions, which makes constraints trivial: move particles, then project positions (Jakobsen-style cloth).

Why energy stays bounded: the map is symplectic, so by backward error analysis it is the exact time-h flow of a modified Hamiltonian H̃ = H + O(h²). Energy H then oscillates within O(h²) of its start for exponentially long times, as long as h is constant.

Code

type V3 = [f64; 3];

pub struct NBody { pub pos: Vec<V3>, pub vel: Vec<V3>, pub mass: Vec<f64>, acc: Vec<V3>, pub g: f64, pub eps2: f64 }

impl NBody {
    pub fn new(pos: Vec<V3>, vel: Vec<V3>, mass: Vec<f64>, g: f64, softening: f64) -> Self {
        let n = pos.len();
        let mut s = NBody { pos, vel, mass, acc: vec![[0.0; 3]; n], g, eps2: softening * softening };
        s.compute_acc();
        s
    }

    /// Pairwise gravity, O(N²), each pair once (Newton's third law).
    fn compute_acc(&mut self) {
        let n = self.pos.len();
        self.acc.iter_mut().for_each(|a| *a = [0.0; 3]);
        for i in 0..n {
            for j in i + 1..n {
                let d: V3 = std::array::from_fn(|k| self.pos[j][k] - self.pos[i][k]);
                let r2 = d[0] * d[0] + d[1] * d[1] + d[2] * d[2] + self.eps2;
                let inv_r3 = self.g / (r2 * r2.sqrt());
                for k in 0..3 {
                    self.acc[i][k] += d[k] * inv_r3 * self.mass[j];
                    self.acc[j][k] -= d[k] * inv_r3 * self.mass[i];
                }
            }
        }
    }

    /// Kick-drift-kick; `self.acc` always holds a(pos) on entry and exit.
    pub fn step(&mut self, h: f64) {
        for (v, a) in self.vel.iter_mut().zip(&self.acc) { for k in 0..3 { v[k] += 0.5 * h * a[k]; } }
        for (x, v) in self.pos.iter_mut().zip(&self.vel) { for k in 0..3 { x[k] += h * v[k]; } }
        self.compute_acc();
        for (v, a) in self.vel.iter_mut().zip(&self.acc) { for k in 0..3 { v[k] += 0.5 * h * a[k]; } }
    }

    /// Yoshida (1990) 4th-order: three Verlet substeps with weights w1, w0, w1 (w0 < 0).
    pub fn step4(&mut self, h: f64) {
        let c = 2f64.powf(1.0 / 3.0);
        let w1 = 1.0 / (2.0 - c);
        let w0 = -c / (2.0 - c);
        self.step(w1 * h);
        self.step(w0 * h);
        self.step(w1 * h);
    }

    pub fn energy(&self) -> f64 {
        let n = self.pos.len();
        let mut e = 0.0;
        for i in 0..n {
            let v = self.vel[i];
            e += 0.5 * self.mass[i] * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
            for j in i + 1..n {
                let r2: f64 = (0..3).map(|k| (self.pos[j][k] - self.pos[i][k]).powi(2)).sum::<f64>() + self.eps2;
                e -= self.g * self.mass[i] * self.mass[j] / r2.sqrt();
            }
        }
        e
    }
}

/// Position (Stormer) Verlet for particles: velocity lives in (x - x_prev).
/// `damping` in [0, 1): 0 = none. Apply constraints to `x` after this call.
pub fn verlet_position(x: &mut [V3], x_prev: &mut [V3], acc: &[V3], h: f64, damping: f64) {
    for i in 0..x.len() {
        for k in 0..3 {
            let cur = x[i][k];
            x[i][k] += (cur - x_prev[i][k]) * (1.0 - damping) + acc[i][k] * h * h;
            x_prev[i][k] = cur;
        }
    }
}

Complexity and constants

• One force evaluation per step; for N-body / MD the force loop (O(N²) direct, O(N log N) tree or PME, O(N) with cutoffs and neighbor lists) is > 95 % of the cost.
• Two-body circular orbit, measured with the code above at 100 steps per orbit: relative energy error peaks at 3.9·10⁻⁶ in every 2,000-orbit window over 20,000 orbits (no drift). RK4 at the same cost (25 steps per orbit) drifts: 4 % energy error after 200 orbits, 44 % after 600, then the orbit breaks up.
• Stability: for the stiffest pair interaction of angular frequency ω, keep ωh well below 2; MD typically uses order of 10-20 steps per fastest period (e.g. 1-2 fs with bonds to hydrogen constrained).
• Yoshida 4: 3 force evaluations per step, error constant larger than the 2nd-order scheme for big h (it takes a negative substep); pays off for tight accuracy.

Efficiency

Cost: one force evaluation per step; in N-body and molecular dynamics the force loop is over 95 % of the cost (O(N²) direct, O(N log N) tree or PME, O(N) with cutoffs). Bound: the force kernel (compute for dense pair loops, memory for neighbor lists). Measured with this page's code on a two-body orbit at 100 steps per orbit: relative energy error 3.9·10⁻⁶ with no drift over 20,000 orbits.

Levers:

• Keep the end-of-step force for the next step — 1 force evaluation per step instead of 2 — always.
• Fixed h near 10-20 steps per fastest period (order of; MD 1-2 fs with hydrogen bonds constrained) — largest stable step — stiff pairs.
• Cutoffs, neighbor lists or tree/PME for the force loop — O(N²) to O(N) or O(N log N) — many particles.
• Structure-of-arrays layout — SIMD in the force loop — particle systems.
• Multiple time stepping (r-RESPA) — long step for slow forces — mixed fast/slow forces.

Naive → standard → state of the art:

• RK4 at 25 steps per orbit · 4 % energy error after 200 orbits, orbit breaks up by 600 (measured above)
• Velocity Verlet, fixed h · bounded energy error, same force cost
• Yoshida 4 (3 force evaluations per step) or r-RESPA · tighter accuracy or fewer expensive force calls

Implementation notes

• Keep h fixed. Changing h every step destroys the bounded-energy behaviour. Games: fixed simulation tick with an accumulator, render with interpolation.
• Store forces with the state so a restart/rollback does not recompute them and the step stays one evaluation.
• Initial half step for leapfrog form: v_{½} = v₀ + h/2·a₀; reporting velocities at integer steps needs the other half kick.
• Structure of arrays (separate x, y, z arrays) for SIMD in the force loop; softening ε avoids the 1/r² singularity in close encounters.
• Reversibility test: run N steps, negate velocities, run N steps; you should return to the start up to rounding. Great regression test for netcode determinism (lockstep, rollback).
• Thermostats/damping break symplecticity by design; apply them as separate operator-split substeps.
• Constraints (rigid bonds, cloth distance limits): SHAKE/RATTLE in MD, iterative position projection in games after verlet_position.

Pitfalls

• Computing forces twice per step (once at start, once at end) because the previous step's forces were not kept.
• Velocity-dependent forces evaluated with the half-step velocity without saying so: no longer the same method, properties lost.
• Damping in position Verlet that depends on h inconsistently: changing the frame rate changes the material.
• Reporting kinetic energy from half-step velocities in leapfrog form: an apparent O(h) energy oscillation that is not real.
• Variable dt from the frame time: energy drift and non-determinism; use a fixed step.

Variants and alternatives

In production

GROMACS (integrator = md is leap-frog, md-vv velocity Verlet), LAMMPS (velocity Verlet run_style verlet), NAMD, OpenMM (VerletIntegrator, LangevinMiddleIntegrator), REBOUND (leapfrog, WHFast), cloth/rope systems in games (Jakobsen's position-based approach in Hitman: Codename 47).

Implementations

Sources

• Loup Verlet, Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules, Physical Review 159(1), 1967. paper
• Ernst Hairer, Christian Lubich, Gerhard Wanner, Geometric Numerical Integration Illustrated by the Störmer-Verlet Method, Acta Numerica 12, 2003. paper
• Ernst Hairer, Christian Lubich, Gerhard Wanner, Geometric Numerical Integration, 2nd ed., Springer, 2006. book
• Haruo Yoshida, Construction of Higher Order Symplectic Integrators, Physics Letters A 150(5-7), 1990. paper
• Thomas Jakobsen, Advanced Character Physics, Game Developers Conference, 2001. paper