4.2. Verlet Integration#
4.2.1. 1. Liouville’s Theorem#
How does the probability distribution evolve?
Define ( \Gamma = (q, p) ) to be the phase space coordinate.
Probability is a conserved quantity, so:
This is analogous to:
Apply the product rule and Hamilton’s equations:
Therefore:
This means phase space is incompressible.
Liouville’s Theorem: The phase-space distribution is constant along trajectories of the system.
Alternative form using Poisson bracket:
The Liouville operator is:
This operator is Hermitian (important property).
4.2.2. 2. Propagators: Symplectic, Time Reversible#
The phase-space trajectories follow from Hamilton’s equations and have a formal solution:
The exponential of the operator is a Taylor series and acts as a propagator.
This propagator is:
Unitary: ( e^{-t \mathcal{L}} = (e^{t \mathcal{L}})^{-1} )
Volume-preserving: Jacobian of transformation is identity
Symplectic: Canonical transformation
Time reversible
4.2.2.1. Not all integrators are symplectic or reversible#
4.2.2.1.1. Example: Explicit Euler#
( x(t + \Delta t) = x(t) + \Delta t \cdot p(t) )
( p(t + \Delta t) = p(t) + \Delta t \cdot f(x(t)) )
Not symplectic or reversible.
4.2.3. 3. Trotter Decomposition#
Split the Liouville operator:
Where:
( \mathcal{L}_1 = \sum_i \frac{\partial H}{\partial p_i} \frac{\partial}{\partial q_i} )
( \mathcal{L}_2 = -\sum_i \frac{\partial H}{\partial q_i} \frac{\partial}{\partial p_i} )
In general, ( \mathcal{L}_1 ) and ( \mathcal{L}_2 ) do not commute.
Using Trotter decomposition:
This gives a single-step propagator:
Symplectic: Each sub-propagator is unitary
Reversible: Symmetric operator sequence Higher-order integrators exist but are more complex.
4.2.4. 4. Velocity Verlet Algorithm#
Apply Trotter decomposition to derive Velocity Verlet:
4.2.4.1. Steps#
Half-step velocity update:
\[ v(t + \frac{\Delta t}{2}) = v(t) + \frac{\Delta t}{2m} F(t) \]Full-step position update:
\[ x(t + \Delta t) = x(t) + \Delta t \cdot v(t + \frac{\Delta t}{2}) \]Half-step velocity update with new force:
\[ v(t + \Delta t) = v(t + \frac{\Delta t}{2}) + \frac{\Delta t}{2m} F(t + \Delta t) \]
4.2.4.2. Equivalent form#
Symplectic
Time reversible
Widely used in MD packages
Other variants: basic Verlet, leapfrog — less convenient than Velocity Verlet.
4.2.4.3. #
4.2.4.4. Derivation of the Verlet algorithm#
Let us make a Taylor solution of trajectory about the current position:
Substitute in for -t: $\( r(t-h) = r(t) - v(t) h + 1/2 a(t) h2 - b(t) h3 + O(h4) \)$
Add the two expansions:
Velocities are eliminated. One can estimate them as:
Time reversal invariance and momentum conservation are built in. This implies that the energy has no drift. Probably the dynamics generated by the Verlet algorithm has a constant of motion which is nearly but not exactly equal to the total energy.
Prove that the Verlet algorithm is time-reversal invariant. That is, show that if we use \(r(t+h)\) and \(r(t)\) as inputs we get \(r(t-h)\) (up to round-off errors). Why does this form have a possible problem with round-off errors?
Velocity Verlet is algebraically equivalent but avoids taking differences between large numbers.
Prove that this is equivalent.
Higher order algorithms are available (Gear etc.) These are useful only when high accuracy is needed and the potential and several derivatives are smooth.
Berendsen and van Gunsteren (in MD simulation of Statistical-Mechanical Systems , 1986) is a study of energy conservation versus time step for the Verlet algorithm and higher order Gear algorithms on a protein. It is shown that for large time steps the Verlet algorithm is considerably more stable. Higher order algorithms have only a the marginal increase in efficiency and actually get worse beyond 7th order.
The following quote is from their article:
How accurate a simulation has to be in order to provide reliable statistical results is still a matter of debate. The purpose is never to produce reliable trajectories. Whatever accuracy the forces and the algorithms have, significant deviations from the exact solution will inevitably occur. Also in the real physical world individual trajectories have no meaning: in a quantum mechanical sense they do not exist and even classically they become unpredictable in the long run due to the nonisolated character of any real system. So what counts are statistical averages. It seems that very little systematic evaluation of algorithms has been done with this in mind. We have the impression that a noise as high as 10% of the kinetic energy fluctuation is still acceptable, although the accuracy of fluctuations may not be sufficient to obtain thermodynamic date from them. With such a level of inaccuracy the Verlet of leap frog algorithm is always to be preferred.
4.2.4.5. Criteria for selecting an integrator#
simplicity (How long does it take to write and debug?)
efficiency (How fast to advance a given system?)
stability (what is the long-term energy conservation?)
reliability (Can it handle a variety of temperatures, densities, potentials?)
vectorization/parallelization (What is the potential efficiency?)
Some additional questions:
Which algorithm will enable me to answer some physical question as quickly as possible? (i.e. minimize the scientist’s time assuming adequate computational resources).
Which algorithm is most efficient? That is which algorithm will calculate some physical quantity to a given accuracy in the smallest amount of CPU time or the least memory (on a given machine)? This is the situation for a mature problem and for writing a proposal to get computer time. You must justify that you are using the machine as efficiently as possible.
Have we balanced statistical errors (going as h-1/2 ) with time step errors which go as hp where p=2,3…? It doesn’t help to have an extremely accurate trajectory if the algorithm is so slow that it cannot sample enough of phase space.
Characteristics of simulations which affect choice of estimator.
Potentials are highly non-linear with discontinuous higher derivatives either at the origin or at the box edge. Harmonic oscillator is not a good test problem!
Small changes in accuracy lead to totally different trajectories. (the mixing or ergodic property)
We need low accuracy because the potentials are not very realistic. Universality saves us: a badly integrated system is probably similar to our original system. This is not true in the field of non-linear dynamics or, for example, in studying the solar system (we know gravity is correct.)
CPU time is totally dominated by the calculation of forces.
Memory considerations are not too important these days (but beware cache misses.)
Energy conservation is important. Energy stability is roughly equivalent to time-reversal invariance. A good test for your code is whether the algorithm can retrace its steps. The energy surface becomes a little fuzzy if energy is not conserved. Rule of thumb: allow 0.01kT fluctuation in the total energy.
Other symmetries may be important, momentum etc.
Note
The nearly universal choice for an integrator is the velocity-Verlet or Verlet algorithm (L. Verlet, Phys. Rev. 159, 98 (1967)).
4.2.5. References#
Reversible multiple time scale molecular dynamics, J. Chem. Phys. 97, 1990 (1992); https://doi.org/10.1063/1.463137, M. Tuckerman, B. J. Berne and G. J. Martyna