Energy conservation: checking your simulation

The idea

With no thermostat (the NVE ensemble), Newton's laws conserve the total energy E = K + U. Kinetic and potential energy trade back and forth, but their sum must stay constant.

A numerical simulation never conserves it perfectly. A healthy run shows only small fluctuations around a flat line. What you must never see is drift (a steady climb or fall) or jumps. These are the usual culprits:

1. Time step too big. The integrator cannot follow the fastest motion (Lesson 1.4).
2. Forces that do not match the energy. The force must be exactly minus the slope of the potential. A typo, a missing factor or a wrong unit breaks this.
3. Jumps at the cutoff. If V(r) is simply cut off at r_c, a pair crossing the cutoff suddenly changes the energy. Shifting the potential so that V(r_c) = 0 removes the jump.

Other common causes are thermostats left on, broken periodic boundaries, or single-precision arithmetic in very long runs.

The algorithm

Test 1: energy conservation (NVE).

1. Turn the thermostat and barostat off.
2. Run a few thousand steps and record K, U and E = K + U.
3. Fit a straight line to E/N versus time. Its slope is the drift.
4. Compute the fluctuation ratio σ(E)/σ(K). A healthy run gives roughly 0.01 or less.

Test 2: forces versus energy (finite difference).

h = 1e-5
x[i] += h;      U_plus  = energy(x)
x[i] -= 2*h;    U_minus = energy(x)
x[i] += h
F_numeric = -(U_plus - U_minus) / (2*h)
assert abs(F[i] - F_numeric) < 1e-6 * abs(F_numeric)

The force test catches coding errors in the forces. It does not catch cutoff jumps or a time step that is too big. Those only show up in the energy test, so you need both.