How should one choose the time step in a molecular dynamics integration?

Question

Choosing too small of a timestep leads to an unrealistic simulation time, whereas too big of a timestep leads to the system not being represented correctly (or, in the case of an algorithm like SHAKE, a SHAKE failure). Given a molecular system to be integrated in time, what are bases on which to decide what value of the time step is ideal? What can be used to assert that the choice is the correct one?

For example for hyperbolic partial differential equations, there is the Courant-Friedrichts-Lewy condition which aids in deciding the necessary timestep size in order to get convergence. Is there an equivalent for the equations of motion in molecular dynamics?

Answer 1

The Rule

The timestep should be less than the period of the fastest vibration by at least 2. In signal processing this is known as Nyquist's theorem.

If a function ${\displaystyle x(t)}$ contains no frequencies higher than B hertz, it is completely determined by giving its ordinates at a series of points spaced ${\displaystyle 1/(2B)}$ seconds apart.

The frequency of a C-H bond is around 3000 $cm^{-1}$. Converting to Hertz this is about 8.99e+13 $Hz$ or a period of 11 femtoseconds.

Therefore, we need a timestep of at least 5 fs but the integrator also introduces some error.

However, even when doing SHAKE (which removes most of the high-frequency vibrations) most MD stick with a 2 fs timestep. For example, see this CHARMM post.

---

So how do we check?

One way you can check if the timestep is okay is to check if there is any drift in a constant energy simulation (NVE). If there is that can mean that the integrator is not behaving time-reversibly. I ran the following with a time-step of 3 fs and no shake and the energy looks constant

I tried to sequentially increase the time-step to demonstrate the drift. However, the energy apparently deviated so fast from constant energy that the energy blew up and OpenMM complained ( this happened at a timestep of 4 fs)

Lastly, I wanted to update this post with this excellent open-access document:

In that document they give excellent advice on the choice of the time step:

• fluctuations of about 1 part in 5000 of the total system energy per twenty time steps are acceptable
• time step size is about 0.0333 to 0.01 of the smallest vibrational period in the simulation

Answer 2

The Nyquist sampling theorem states that the time step must be half or less of the period of the quickest dynamics. This is the absolute maximum time step that can capture the quickest dynamics at all, and it is usually recommended to choose a much smaller time step (often chosen to be between 0.1 and 0.2 of the shortest period).

Typical time steps range from 0.25 fs for systems with light nuclei (such as hydrogen), to 2 fs or greater for systems with more massive nuclei. (NB for quantum nuclear dynamics, e.g. using path-integrals, you may need even shorter time steps for light nuclei, especially at low temperatures.)

In order to work out how reasonable your time step is, it is common to monitor the appropriate "conserved quantity" over the course of a short simulation. The relevant "conserved quantity" depends on the ensemble used for the dynamics; for NVE it is simply the total energy, but for NVT, NPH etc. there are extra terms from the thermostat and barostat. Whatever the appropriate quantity is, it is conserved in the ideal dynamics, so monitoring the numerical drift in this quantity is a good indicator of the error in your dynamics. A reasonable rule of thumb is that the long-term drift in the "conserved quantity" should be less than 10 meV/atom/ps for qualitative results, and 1 meV/atom/ps for "publishable" results.

The time step alone is not sufficient in order to determine the accuracy of a molecular dynamics calculation. There are three other considerations which come to mind:

• The time-integration algorithm. It is very important that your integration algorithm is symplectic, which essentially means it is time-reversible (energy conservation is directly related to the time-reversal symmetry of the physical laws). If you use a non-symplectic integrator then you will need a much shorter time step, even if the integrator is supposedly more accurate (e.g. high-order predictor-corrector methods). Velocity-Verlet is probably the most common method used, and is symplectic.

• If you are doing ab initio dynamics, then the energy and forces are computed from an iterative solution to the quantum mechanical equations. Because the solutions are not converged perfectly there is numerical noise in the energies and forces, and this will cause additional drift in the conserved quantities.

• Again, for ab initio dynamics it is common to accelerate the iterative solution of the quantum mechanical equations by extrapolating the wavefunction, density etc. from the previous time step to the present time step. Many of these extrapolation schemes are not time-reversible, which breaks the overall symplecticity of the method, causing a greater drift in the conserved quantity and necessitating a smaller time step. This has been addressed by methods such as extended-Lagrangian dynamics, where the wavefunction and density degrees of freedom are also propagated with a symplectic algorithm.

Answer 3

Other answers explaining Nyquist's Sampling Theorem are completely correct. Read those before mine. I'd like to add some information about how details of how you run your simulation can affect the Nyquist frequency, and therefore, in practice, affect the fundamental question of how to choose your time step.

As one can easily imagine, tricks that let you simulate a longer time for the same computational effort have been a major area of research. Here are a few relevant points I know of (no effort to be comprehensive -- just the articles that spring to my mind!):

• Integration algorithms can affect reasonable step size. Some integrators (a.k.a. propagators) can exhibit fortuitous cancellation of errors that allows a larger step size. The main example of this that I know is the Langevin "BAOAB" splitting. For details, see this paper: https://doi.org/10.1093/amrx/abs010 (arXiv: https://arxiv.org/abs/1203.5428)
• You can add constraints to increase step size (and an integrator to exploit that). The purpose of constraints is usually to allow larger step sizes -- by removing faster vibrational motions, you change the Nyquist frequency. This can be further exploited by integration algorithms to allow even larger time step. For example, the article introducing the geodesic BAOAB Langevin integrator suggests that it should allow step sizes of up to 8 fs for biomolecular system! Here's the reference for g-BAOAB: https://royalsocietypublishing.org/doi/full/10.1098/rspa.2016.0138
• There is good evidence that you can redistribute mass to increase step size. In general, the highest frequency vibrations are those involving a hydrogen atom. For this reason, it is common to constrain bonds to hydrogen atoms. However, another approach is "hydrogen mass repartitioning" (HMR), where mass from the heavy atom is transferred to the hydrogen, keeping the same total mass without requiring constraints, and resulting in a slower Nyquist frequency. Here's a reference for that: https://pubs.acs.org/doi/10.1021/ct5010406

Answer 4

Having read the above answers I think it's important to point out why you generally use a time step below the upper limit provided by Nyquist's Theorem.

The theorem itself describes the maximum number of signals it's possible to send through a telegraph cable and still have them be discernable from one another. That makes the limiting case a sampling frequency that gives you two points along a wave period, one in the positive section and one in the negative section. This is fine for a signal being processed as Ones and Zeroes, but for atoms in molecular dynamics whose interactions are governed by a continuous pair potential this is extremely low resolution.

Therefore, a smaller time step is often chosen than the limiting frequency indicated by Nyquist's theorem, because such a low resolution can still produce unphysical results that cause a loss of energy conservation.