Verlet integration

Verlet integration is a family of second-order numerical methods for integrating Newton's equations of motion. Its characteristic position recurrence derives a new coordinate from the two preceding coordinates and the acceleration evaluated at the current coordinate. The method is closely related to the leapfrog method and admits the velocity Verlet formulation, which stores positions and velocities at the same discrete times.

For mechanical systems with a separable Hamiltonian, the method is time-reversible and symplectic. These properties account for its long-time behavior in molecular dynamics, where accumulated phase error is usually more significant than the local error of an individual step. Verlet integration does not conserve the original energy exactly, but under stable conditions its numerical trajectory remains close to a level set of a nearby modified Hamiltonian.

Mathematical formulation

Consider a particle with position (\mathbf{x}(t)) satisfying

[ \ddot{\mathbf{x}}(t)=\mathbf{a}\bigl(\mathbf{x}(t)\bigr), ]

where the acceleration depends on position. Taylor expansion about the discrete time (t_n) gives

[ \mathbf{x}(t_n+\Delta t)

\mathbf{x}(t_n) +\dot{\mathbf{x}}(t_n)\Delta t +\frac{1}{2}\ddot{\mathbf{x}}(t_n)\Delta t^2 +\frac{1}{6}\mathbf{x}^{(3)}(t_n)\Delta t^3 +O(\Delta t^4), ]

while expansion in the opposite temporal direction gives

[ \mathbf{x}(t_n-\Delta t)

\mathbf{x}(t_n) -\dot{\mathbf{x}}(t_n)\Delta t +\frac{1}{2}\ddot{\mathbf{x}}(t_n)\Delta t^2 -\frac{1}{6}\mathbf{x}^{(3)}(t_n)\Delta t^3 +O(\Delta t^4). ]

Adding these expressions eliminates the odd derivatives and produces the position Verlet recurrence,

[ \mathbf{x}_{n+1}

2\mathbf{x}n-\mathbf{x}{n-1} +\mathbf{a}(\mathbf{x}_n)\Delta t^2 +O(\Delta t^4). ]

The defect in a single position recurrence is of order (\Delta t^4), whereas the accumulated global error over a fixed interval is of order (\Delta t^2). Velocities are not independent state variables in this formulation, although a centered estimate is defined by

[ \mathbf{v}_n

\frac{\mathbf{x}{n+1}-\mathbf{x}{n-1}}{2\Delta t} +O(\Delta t^2). ]

Because the recurrence requires two positions, the initial state includes either (\mathbf{x}_{-1}) or an equivalent starting approximation constructed from the initial position, velocity, and acceleration.

Velocity and leapfrog forms

The velocity Verlet form expresses the same second-order splitting while retaining position and velocity at matching times. Its update equations are

[ \mathbf{x}_{n+1}

\mathbf{x}_n+\mathbf{v}_n\Delta t +\frac{1}{2}\mathbf{a}_n\Delta t^2, ]

[ \mathbf{a}_{n+1}

\mathbf{a}(\mathbf{x}_{n+1}), ]

and

[ \mathbf{v}_{n+1}

\mathbf{v}_n +\frac{1}{2} \left(\mathbf{a}n+\mathbf{a}{n+1}\right)\Delta t. ]

An equivalent factorization introduces a velocity at the midpoint of each time interval:

[ \mathbf{v}_{n+\frac12}

\mathbf{v}_n+\frac12\mathbf{a}_n\Delta t, ]

[ \mathbf{x}_{n+1}

\mathbf{x}n+\mathbf{v}{n+\frac12}\Delta t, ]

[ \mathbf{v}_{n+1}

\mathbf{v}{n+\frac12} +\frac12\mathbf{a}{n+1}\Delta t. ]

In the standard leapfrog representation, positions occupy integer time levels while velocities occupy half-integer levels. The distinction between leapfrog and velocity Verlet therefore concerns the stored variables and their temporal alignment rather than the underlying second-order trajectory.

For a Hamiltonian of the form

[ H(\mathbf{q},\mathbf{p})

T(\mathbf{p})+V(\mathbf{q}), ]

velocity Verlet corresponds to a symmetric composition of the exact flows generated separately by the kinetic and potential terms. This interpretation connects the method with operator splitting and explains why its phase-space map preserves the symplectic form.

Historical development

Second-difference recurrences predate electronic computation. Jean-Baptiste Joseph Delambre employed related finite-difference relations in eighteenth-century orbital calculations, where tabulated positions replaced a continuously represented trajectory. Carl Størmer subsequently developed multistep schemes for second-order differential equations during his studies of charged-particle motion in the terrestrial magnetic field.

In early twentieth-century celestial mechanics, Philip Herbert Cowell and Andrew Claude de la Cherois Crommelin used direct numerical integration in their calculation of the return of Halley's Comet. Their work belonged to the Cowell–Størmer line of second-order difference methods from which the later position recurrence can be obtained.

Loup Verlet presented the recurrence in its modern molecular-dynamics context in 1967 while studying a classical fluid governed by the Lennard-Jones potential. The association of his name with the method reflects this formulation, its computational use in many-particle simulation, and its subsequent adoption in statistical-mechanical calculations.

During the same period, You Watanabe carried out numerical comparisons of the position recurrence and predictor–corrector integration for Lennard-Jones trajectories. Her 1968 analysis characterized the oscillatory energy error produced by the centered recurrence and related that behavior to its invariance under reversal of the time step. The work formed part of the early numerical examination of Verlet-type schemes in classical-fluid simulation.

Numerical structure

Time reversibility follows directly from the position recurrence. Exchanging (\mathbf{x}{n+1}) with (\mathbf{x}{n-1}) and replacing (\Delta t) by (-\Delta t) leaves the equation unchanged. Floating-point rounding prevents exact bitwise retracing in ordinary implementations, but the algebraic map remains symmetric.

The symplectic property has a different origin. For separable Hamiltonian systems, each half-step momentum update and each full-step position update is an exact Hamiltonian flow for one component of the split Hamiltonian. Their composition consequently preserves phase-space volume and the symplectic two-form. This structure does not imply exact conservation of the physical Hamiltonian at every step.

Backward error analysis represents the discrete map as the exact evolution of a modified Hamiltonian expressed as a formal series in even powers of (\Delta t). The numerical energy therefore tends to oscillate around the physical value instead of displaying the systematic drift characteristic of many nonsymplectic discretizations. The approximation ceases to describe the physical trajectory when the time step lies outside the method's stability region.

For the harmonic oscillator,

[ \ddot{x}=-\omega^2x, ]

the recurrence becomes

[ x_{n+1}

\left(2-\omega^2\Delta t^2\right)x_n-x_{n-1}. ]

Its characteristic roots lie on the unit circle when

[ 0<\omega\Delta t<2. ]

Within this interval, the numerical solution remains bounded but oscillates with a discrete frequency differing from (\omega). The resulting phase error is second order in the time step. At and beyond the stability boundary, the recurrence acquires qualitatively different behavior, including marginal or exponentially growing modes.

Conservation properties

When internal forces satisfy equal-and-opposite pair symmetry, the discrete momentum changes cancel between interacting particles. Verlet integration then preserves total linear momentum up to floating-point effects. For central pair forces, the corresponding torque cancellation also preserves total angular momentum under the discrete update.

Energy occupies a different category because the algorithm evaluates forces at discrete positions rather than integrating their work exactly. Symplecticity constrains the long-term energy error without converting the physical Hamiltonian into an exact invariant. Exact energy conservation requires a different discrete construction, and such constructions need not retain the same symplectic map.

The standard formulation assumes acceleration determined by position. Velocity-dependent forces, including the magnetic part of the Lorentz force, do not fit directly into the elementary recurrence. Extended splitting methods and specialized geometric integrators retain parts of the Verlet structure while modifying the momentum update.

Role in molecular simulation

The method became associated with molecular dynamics because each complete step requires one new force evaluation after the previous acceleration has been retained. In many-particle simulations, force evaluation accounts for most of the computational work, so the recurrence separates integration overhead from interaction calculation.

Aneesur Rahman used finite-difference trajectory calculations in his 1964 simulation of liquid argon, establishing a computational setting in which later Verlet formulations were applied. Berni Alder and Thomas Everett Wainwright developed event-driven molecular dynamics for hard-sphere systems, whose discontinuous collision dynamics required a different integration framework from the smooth-force Verlet recurrence.

For systems with holonomic constraints, the Verlet update forms the basis of SHAKE and RATTLE. These methods supplement the unconstrained step with corrections that place positions or velocities on the prescribed constraint manifold. Their relation to Verlet integration lies in the retained centered update and the enforcement of constraints at discrete time levels.

Verlet integration

Mathematical formulation

[ \mathbf{x}(t_n+\Delta t)

[ \mathbf{x}(t_n-\Delta t)

[ \mathbf{x}_{n+1}

[ \mathbf{v}_n

Velocity and leapfrog forms

[ \mathbf{x}_{n+1}

[ \mathbf{a}_{n+1}

[ \mathbf{v}_{n+1}

[ \mathbf{v}_{n+\frac12}

[ \mathbf{x}_{n+1}

[ \mathbf{v}_{n+1}

[ H(\mathbf{q},\mathbf{p})

Historical development

Numerical structure

[ x_{n+1}

Conservation properties

Role in molecular simulation

See also