The time step

The idea

Every MD step assumes the force stays constant during Δt. That is only true if the atoms barely move in one step. So Δt must be short compared with the fastest motion in the system.

The fastest motions are bond vibrations involving hydrogen, the lightest atom. A C–H bond vibrates back and forth about once every 11 fs, and an O–H bond about every 9 fs. To follow a vibration well you need about 10–20 steps per period, which gives the familiar Δt ≈ 1 fs for all-atom MD.

• Too big: the steps overshoot. The energy drifts, and past a limit (about T/π) the simulation explodes.
• Too small: everything is accurate, but you need more steps. 1 ns at 1 fs is already a million steps.

Common choices in practice:

• 1 fs: all-atom, flexible bonds.
• 2 fs: bonds to hydrogen held fixed with constraints (SHAKE/LINCS, Lesson 4.4).
• 4 fs: hydrogen mass repartitioning, which makes H heavier.
• 10–40 fs: coarse-grained models with no light atoms at all.

The algorithm

How to choose and check Δt for a new system:

1. Find the fastest motion. Usually the stiffest bond to a light atom: O–H, N–H, C–H.
2. Estimate its period from its vibrational wavenumber: T = 1/(c ν̃). For 3000 cm⁻¹, T ≈ 11 fs.
3. Pick Δt ≈ T/20 to T/10. Never go near T/π, where the update becomes unstable.
4. Run a short NVE test (no thermostat) and plot the total energy. It should fluctuate but not drift.
5. If it drifts, halve Δt or constrain the bonds to hydrogen, then test again.

for dt in [0.5, 1.0, 2.0]:          # fs
    run_nve(dt, n_steps=10_000)
    plot(total_energy)            # flat = good, drifting = too big

The code tab does exactly this test for a single C–H bond.