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

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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.