LESSON NOTES · 04

4. Derive velocity Verlet

Position in the course: Lesson 4 of 11. Complete the preceding derivation and use the explained exercises to check understanding.

Prerequisites: Hamiltonian trajectories; Taylor expansions

1. Expand the exact motion

Use \(h\) for an MD timestep, distinct from Planck’s constant. For smooth position-dependent acceleration, Taylor-expand position, velocity and acceleration. Average the acceleration at the two endpoints to reproduce the velocity expansion through second order.

\[ \begin{aligned}x(t+h)&=x+hv+\tfrac12h^2a+O(h^3),\\a(t+h)&=a+h\dot a+O(h^2),\\v(t+h)&=v+ha+\tfrac12h^2\dot a+O(h^3)=v+\tfrac h2[a(t)+a(t+h)]+O(h^3).\end{aligned} \]

2. Discard higher terms to define the algorithm

First update position, recompute forces there, then update velocity. Equivalently: half kick, full drift, recompute force, half kick. For a separable Hamiltonian the symmetric update is time reversible and symplectic, with global error \(O(h^2)\) over a fixed smooth time interval. Stable resolved conservative examples often have bounded oscillatory energy error, not exact conservation of the physical Hamiltonian.

\[ x_{n+1}=x_n+hv_n+\tfrac12h^2a_n,\quad a_{n+1}=F(x_{n+1})/m,\quad v_{n+1}=v_n+\tfrac h2(a_n+a_{n+1}). \]

3. Derive harmonic stability

For \(a=-\omega^2x\), eliminate velocity to get the central position recurrence. Insert an oscillatory solution \(x_n=e^{in\theta}\). A real phase requires \(0<h\omega<2\); the endpoint is marginal and unusable in practice. Stability alone does not guarantee accurate phase: the numerical frequency is \(\theta/h\).

\[ x_{n+1}-2x_n+x_{n-1}=-h^2\omega^2x_n,\quad2\cos\theta-2=-h^2\omega^2,\quad\cos\theta=1-\tfrac12(h\omega)^2. \]

Derive velocity Verlet — schematic under the stated model assumptions

Teaching schematic drawn from the equations or algorithm steps in this lesson; not measured data.

4. Worked checks: stability versus accuracy

At \(h\omega=1\), \(\theta=\pi/3\) rather than one radian per step: stable but phase inaccurate. If the asymptotic second-order trajectory error is \(0.04\) at \(h\), it should approach \(0.01\) at \(h/2\) for the same final time and error definition. Resolve the fastest physical mode, not just visible slow motion.

5. Assumptions, limits and checks

These equations omit constraints, velocity-dependent forces, thermostats and stochastic noise, which need suitable splitting or specialized updates. Sharp cutoffs, missed neighbors and energy–force inconsistency can spoil energy behavior even at small \(h\). Does bounded energy error prove trajectory accuracy indefinitely? No; phase and chaotic trajectory errors can grow despite good energy behavior.

6. Why reversibility does not imply exact energy

The half-kick/full-drift/half-kick map is a composition of exactly solvable Hamiltonian subflows. Reversing \(h\) reverses their order and undoes the map, which proves time reversibility under the stated smooth deterministic assumptions. Each subflow preserves phase-space volume and the symplectic form. This structure helps explain bounded energy errors, but it does not mean \(H(R,p)\) itself is constant at every step.

A nearby modified Hamiltonian can describe the discrete motion perturbatively over suitable time intervals. Its correction starts at order \(h^2\) for this symmetric second-order splitting. Smoothness, timestep stability and bounded motion are essential; neighbor-list mistakes, sharp force jumps or uncontrolled electronic force noise break the premise. Long chaotic trajectories should be compared statistically, since pointwise trajectory agreement is eventually lost even with good numerical methods.

Worked phase expansion: from \(\cos\theta=1-(h\omega)^2/2\), use \(\theta=2\arcsin(h\omega/2)\). The series gives \(\theta/h=\omega+h^2\omega^3/24+O(h^4)\). The numerical oscillator therefore runs slightly fast for small stable \(h\). After many periods this small frequency error accumulates into a substantial phase shift while the energy error can remain bounded.

Exercise: if the timestep halves, how does this leading frequency error change? It decreases by four. This is a stronger accuracy diagnostic than merely checking that the trajectory did not diverge. Compare at the same physical final time and keep force tolerances fixed so the timestep contribution can be isolated.

7. References and further reading

Derivations and numerical examples here are original teaching constructions, not copied passages or reported research data.

Quantum mechanics · Monte Carlo · Quantum Monte Carlo · CP2K · Quantum-Espresso

Previous lesson · Course overview · Next lesson