Lab 1.8 · 30 min

Lab: integrator extras (shadow energy, precision, RESPA)

Velocity Verlet keeps the energy steady for billions of steps, while “better” fourth-order methods slowly leak it away. The reason is a hidden conserved quantity, the shadow Hamiltonian. We compute it by hand, check it in an LJ cluster, see how single precision breaks energy conservation (and why it depends on where the atoms are), and speed up a stiff system with multiple time steps (RESPA), including the resonance that limits it.

Compute the shadow Hamiltonian that Verlet conserves, and see why Verlet does not drift but RK4 does
See how float32 rounding causes energy drift, and why it grows with the distance from the origin
Split fast and slow forces with RESPA, count the savings, and find the resonance that limits the outer step

The setup

Velocity Verlet (Lab 1.5), one step of size h, mass 1:

p ← p + (h/2) F(x) x ← x + h p p ← p + (h/2) F(x)

Systems: a harmonic oscillator (ω = 1), a 13-atom LJ cluster (the icosahedron, no cutoff, k_BT ≈ 0.25), a fast-plus-slow spring, and a liquid of 32 stiff dimers. All in LJ units.

Step 1

The shadow Hamiltonian, by hand

Why: Verlet does not conserve the true energy E exactly: it wobbles by about h². But for a harmonic oscillator it conserves a slightly different quantity exactly, forever. That hidden quantity is why the energy never drifts away.

E = p²/2 + x²/2 H̃ = p²/2 + (x²/2)(1 − h²/4) (conserved by Verlet) h = 0.5, start x = 1, p = 0 step x p E H̃ 0 +1.000000 +0.000000 0.500000 0.468750 1 +0.875000 −0.468750 0.492676 0.468750 2 +0.531250 −0.820312 0.477570 0.468750 3 +0.054688 −0.966797 0.468843 0.468750 4 −0.435547 −0.871582 0.474678 0.468750 step 1 in detail: p = 0 − 0.25 × 1 = −0.25 x = 1 + 0.5 × (−0.25) = 0.875 p = −0.25 − 0.25 × 0.875 = −0.46875 H̃ = 0.46875²/2 + (0.875²/2)(0.9375) = 0.468750

shadow frequency (2/h) arcsin(h/2) = 1.0107: 1.07% fast stable only while h ω < 2 h = 2.1: 1 − h²/4 = −0.10 < 0, no closed orbit after 20 steps |x| = 1.5×10⁵

What the numbers say: E moves between 0.469 and 0.5 but H̃ stays at 0.468750 to every digit. Verlet follows the exact orbit of a slightly different system (a spring 1 − h²/4 times as stiff). Its energy error is bounded, and its only lasting error is a small phase shift (1% at h = 0.5). Above h ω = 2 no shadow orbit exists and the motion blows up: the stability limit of Lesson 1.4.

Step 2

The shadow energy in an LJ cluster, and why RK4 drifts

Why: For any smooth force, Verlet conserves a shadow energy that differs from E by terms of order h². Adding the first correction should remove most of the wobble, and what is left should shrink as h⁴. A non-symplectic method has no shadow energy at all, however accurate it is per step.

H̃ = E + h² ( p·∇²U·p / 12 − |F|²/24 ) (reduces to the exact H̃ of Step 1 for a spring, to order h²) 13 LJ atoms, icosahedron (U = −44.3268), k_BT = 0.25, t = 50

h	std of E	std of H̃	ratio
0.008	1.16×10⁻³	7.5×10⁻⁶	154
0.004	2.99×10⁻⁴	5.0×10⁻⁷	594
0.002	7.70×10⁻⁵	3.2×10⁻⁸	2429

halving h: E wobble ÷ 3.88 (h²: ÷ 4) H̃ wobble ÷ 15–16 (h⁴: ÷ 16)

Verlet vs RK4 (fourth order) on the oscillator, h = 0.2, 100 000 steps:

method	E after 10 steps	E range over run	E at the end
Verlet	0.49675	0.49500 – 0.50000	0.49700
RK4	0.4999951	0.45768 – 0.50000	0.45768

RK4 energy factor per step: 1 − h⁶/72 = 0.999999111 after 100 000 steps: 0.999999111^100000 = 0.915

What the numbers say: One correction term makes the conserved quantity 150–2400 times steadier, with the h⁴ scaling expected. RK4 is a thousand times more accurate per step than Verlet, yet it loses 8.5% of the energy over the run, steadily and without limit; Verlet stays within 1%. Over long MD runs, having a conserved shadow energy matters more than the order of the method.

Step 3

Single vs double precision

Why: GPUs are much faster in single precision (float32, about 7 digits). Every position update is then rounded. If the rounding errors were truly random they would cancel, but a coordinate far from the origin has fewer digits left for the small step it takes each time.

float32 keeps 24 binary digits: spacing = 2⁻²³ × 2^(exponent) at 1: 1.19×10⁻⁷ at 100: 7.63×10⁻⁶ at 10: 9.54×10⁻⁷ at 1000: 6.10×10⁻⁵ float64 at 1000: 1.1×10⁻¹³ step size: v h ≈ 0.45 × 0.004 = 0.0018 per step relative rounding of that step: near 0: 1×10⁻⁷ / 0.0018 = 7×10⁻⁵ at 1000: 6×10⁻⁵ / 0.0018 = 3%

precision	cluster at	drift of E over t = 200	rms of E − E₀
float64	0	−1.9×10⁻⁶	1.21×10⁻³
float32	0	−3.7×10⁻⁶	1.20×10⁻³
float32	10	−2.8×10⁻³	0.97×10⁻³
float32	100	−1.4×10⁻²	6.1×10⁻³
float32	1000	+5.3×10⁻²	5.7×10⁻²

What the numbers say: Near the origin, float32 is as good as float64: the wobble is the h² shadow fluctuation, not rounding. At distance 1000 the drift is 27 000 times larger. That is why GROMACS mixed precision stores coordinates relative to the box, keeps the box small (nanometres), and accumulates some sums in double precision, and why long NVE runs for transport properties are often checked in double precision.

Step 4

RESPA by hand: fast and slow forces on different clocks

Why: The time step is set by the fastest motion (bond vibrations), but the most expensive forces (long-range electrostatics, non-bonded interactions) change slowly. RESPA (reversible reference system propagator) applies the slow force once per long step H and the fast force every short step h = H/n. It is still symplectic and time reversible.

outer step H: p ← p + (H/2) F_slow n inner Verlet steps of h = H/n with F_fast p ← p + (H/2) F_slow x'' = −(ω_f² + ω_s²) x, ω_f = 10, ω_s = 1, ω = 10.05 H = 0.1, n = 2 (h = 0.05), start x = 1, p = 0, E = 50.5 slow half kick: p = 0 − 0.05 × 1 × 1 = −0.05000 inner 1: kick p = −0.05 − 0.025 × 100 × 1 = −2.55000 drift x = 1 + 0.05 × (−2.55) = +0.87250 kick p = −2.55 − 0.025 × 100 × 0.8725 = −4.73125 inner 2: kick −6.91250, drift +0.52687, kick −8.22969 slow half kick: p = −8.25603, E = 48.10 (slow force: 1 evaluation; fast force: 2)

Resonance scan, inner step ≈ 0.005, t = 400:

H	H ω / π	max \|E/E₀ − 1\|
0.05	0.16	1.5×10⁻³
0.20	0.64	2.0×10⁻²
0.30	0.96	0.22
0.3126	1.00	10⁶ (blows up)
0.33	1.06	0.19
0.40	1.28	2.3×10⁻²
0.625	2.00	10⁶ (blows up)

What the numbers say: The inner steps are tiny, yet the run explodes when the outer step equals half a period of the fast motion (H = π/ω) or a multiple of it: the slow kicks then always hit the vibration at the same phase and pump energy into it. In practice the outer step must stay well below half the fastest period, about H < 0.3 × period.

Step 5

RESPA in a liquid of stiff molecules

Why: Now a real case: 32 dimers with stiff bonds (k = 2000) in a 3D LJ liquid. The bonds are cheap and fast; the LJ forces are expensive and slower. Counting LJ evaluations per unit time measures the cost.

bond: ω = √(k/μ) = √(2000/0.5) = 63.2, period 0.0993 Verlet limit h < 2/ω = 0.032; resonance H = π/ω = 0.050 64 atoms, ρ = 0.6, T = 1.5, NVE for t = 10

method	LJ evals / time	bond evals / time	std of E/N	drift E/N
Verlet h = 0.001	1000	1000	8.0×10⁻⁵	+2×10⁻⁶
Verlet h = 0.004	250	250	1.3×10⁻³	+3.8×10⁻⁴
Verlet h = 0.008	125	125	5.9×10⁻³	−1.1×10⁻³
RESPA H = 0.004, n = 4	250	1000	3.6×10⁻⁴	+1.2×10⁻⁴
RESPA H = 0.008, n = 8	125	1000	1.4×10⁻³	−6.8×10⁻⁴
RESPA H = 0.02, n = 20	50	1000	5.5×10⁻²	+0.18
RESPA H = 0.04, n = 40	25	1000	blew up

What the numbers say: RESPA with H = 0.008 matches Verlet with h = 0.004 at half the LJ cost, and at equal LJ cost it is 4 times steadier. The gain stops near H ≈ 0.01: beyond that the LJ forces themselves are too fast for the long step (H = 0.04 blows up, below the bond resonance at 0.05). In real codes the slow part is usually the Ewald reciprocal-space force (LAMMPS run_style respa, OpenMM MTSIntegrator); constraints (SHAKE, Lesson 4.4) are the other way to remove the fastest motions.

Summary

topic	key number	lesson
shadow Hamiltonian	H̃ = 0.468750 at every step	Verlet follows the exact orbit of a nearby system
stability	h ω < 2	no shadow orbit beyond it
shadow correction	E wobble ∝ h², H̃ wobble ∝ h⁴	the conserved quantity exists for real forces
RK4 vs Verlet	−8.5% vs bounded, 10⁵ steps	symplectic beats high order over long runs
float32	drift ×27 000 at distance 1000	keep coordinates small; check in double
RESPA resonance	blows up at H = π/ω	outer step well below half the fast period
RESPA savings	same accuracy, ½ the LJ cost	limited by the slow force’s own time scale