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

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?

• maybe someone can add to this in a formal answer. One way you can check if the timestep is okay is to check if there is any drift in a constant energy simulation. If there is that can mean that the integrator is not behaving time-reversibly. There is also the rule of thumb, the timestep should be faster than the period of the fastest vibration by at least 2. – Cody Aldaz May 10 at 19:30
• That "rule of thumb" in Cody's comment, is sometimes called Nyquist's theorem. – Nike Dattani May 10 at 20:42
• @CodyAldaz I don't have much experience with MD, do you want to answer it? – Nike Dattani May 11 at 17:25

## 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 $$x(t)}$$ contains no frequencies higher than B hertz, it is completely determined by giving its ordinates at a series of points spaced $$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
• The system is a cyclohexane molecule – Cody Aldaz May 11 at 18:46
• A reasonable rule of thumb is that the 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 relevant "conserved quantity" depends on the ensemble used, and will only be the total energy for NVE. Light nuclei typically need shorter time steps, perhaps as low as 0.25 fs for accurate hydrogen dynamics, whereas for more massive nuclei you can often use timesteps around 2 fs or greater. – Phil Hasnip May 14 at 2:03
• @PhilHasnip This would be a good additional answer – Cody Aldaz May 14 at 2:14
• I've expanded it into an additional answer. – Phil Hasnip May 15 at 2:14

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.