# Why aren't Runge-Kutta methods used for molecular dynamics simulations?

One of the most used schemes for solving ordinary differential equations numerically is the fourth-order Runge-Kutta method. Why isn't it used to integrate the equation of motion of particles in molecular dynamics?

• RK4 doesn't preserve conservation of energy, for one. IIRC it tends to lose energy over time. Symplectic integrators do better here.
– TLW
May 10, 2022 at 3:39
• What methods are used, then? May 10, 2022 at 14:01
• @TLW many (but not all) molecular dynamics applications do not need strict energy conservation, very often we use an explicit, significant damping term and use it to find the minimum energy configuration. In these cases nobody is interested in the accuracy of the path taken to the bottom, only that it successfully finds its way there. cf. FIRE: Fast Inertial Relaxation Engine for Optimization on All Scales If you can find an actual example of someone using a symplectic integrator in MD and of it "doing better" that would be interesting!
– uhoh
May 11, 2022 at 1:01
– uhoh
May 11, 2022 at 1:05

In molecular dynamics simulations, the overwhelming part of the computational time is spent evaluating forces. For this reason, since the very beginning of the method, the algorithms of choice were those requiring the minimum evaluation of forces per unit of time.

Higher-order algorithms are usually not necessary for two reasons:

• in molecular dynamics simulations, one is not interested in the maximum accuracy of the trajectories but in an efficient sampling of the phase space compatible with the physical constraints. For this reason, usually, lower-order but symplectic schemes are preferred.
• the increased accuracy of higher-order algorithms does not allow a proportional increase of the time step size. The reason is the stiff nature of the interatomic forces at short distances.

Based on these considerations, the fourth-order Runge-Kutta method is usually excluded. It is not time-reversible nor symplectic and requires four force evaluations per step. Even the existing symplectic Runge-Kutta methods are not appealing when compared with methods of the same order using fewer force evaluations.

• The Verlet method is 3rd order and symplectic without ANY additional force evaluations. RK4 is 4th order with 4 times the cost. Also, high accuracy is impossible when we are making crude approximations from quantum mechanics. May 11, 2022 at 6:57
• @KevinKostlan When we say that Verlet is 3rd order, it means that the local error is third order in the time step. The global error is 2nd order. RK4 has a fourth-order global and a fifth-order. local error. Therefore, the fair comparison should be between the global 2nd order Verlet and the global 2nd order RK that requires two evaluations of the force per step. Better than four, but worse than one :-) Finally, I would add that Verlet's method does not allow general velocity-dependent forces. Its generalization for this case (Groot-Warren) requires two force evaluations per step. May 11, 2022 at 8:34
• Good answer but might be useful to add a definition of "symplectic" for those unfamiliar with this term of art.
– WOPR
May 11, 2022 at 14:33

In molecular dynamics you have to evolve a large number of particles, so efficiency is important. The Verlet algorithm provides numerical stability, as well as other properties that are important in molecular dynamics: time reversibility and preservation of the symplectic form on phase space, at no significant additional computational cost over the most basic integrator, i.e., the Euler method. Therefore, if you want to use Runge Kutta (that is computationally more expensive than Euler), you first ask yourself if it would perform better than Verlet.

GiorgioP has the right answer, but I just want to emphasize a point. We use different time integration methods for different problems. For example, consider these two problems.

1. You want to know the location of satellite orbiting the Earth one year from now.
2. You want to know if one drug binds more strongly to its target protein than another drug.

In the first case, you use a high-order Runge-Kutta method, or something else that gives you very accurate trajectories. In the second case, you are more interested in getting accurate statistical properties. You don't need to know the exact position of a particular water molecule after 100 nanoseconds. Often for molecular simulation, we are only interested in equilibrium thermodynamic properties like free energies and radial distribution functions. To calculate such statistical properties, it is important that the integrator be reversible and symplectic so that equilibrium can be reached and so that we correctly sample a valid thermodynamic ensemble. Also, to get good statistical accuracy, we need a lot of sampling, so it's also important that the method be computationally efficient.

We usually don't want to exactly follow the trajectories of molecules at the nanoscale, since these usually can't be compared to experiment because of the sensitivity to initial conditions and perturbations from outside (thermal fluctuations, vibrations). In experiment, we can't control initial conditions at the molecular scale very well (certainly not to the sub-picometer precision that would be required for reproducible trajectories on the picosecond timescale) and can't isolate systems from thermal fluctuations of their surroundings. If we are comparing to experiments happening in beakers of water or inside living organisms, there is no way to control outside influences even if we wanted to.

Moreover, in molecular dynamics simulations, we often use thermostats and barostats that slightly modify the dynamics (and therefore trajectory) in artificial ways, so that we can sample constant temperature ($$NVT$$) or constant pressure ensembles ($$NpT$$). These ensembles are usually more useful for comparing to experiment than the ensemble given by applying Newton's laws unmodified (constant particle number, volume, and energy, $$NVE$$).

We may use enhanced sampling methods like replica exchange, which lead to even more dramatic departures from physically realistic trajectories. For example, in temperature replica exchange (also known as parallel tempering), you run 10 replicas of the system at different temperatures and periodically attempt to exchange coordinates between the different temperature replicas using the Metropolis criterion. When exchanges are accepted, the trajectory at a given temperature changes discontinuously (molecules may appear teleport around). If we are only interested average equilibrium quantities, we don't care about thes artificial jumps, only that we are accurately sampling an $$NVT$$ (or $$NpT$$) ensemble.