An Explanation of the Flying Ice Cube Effect: How Simple Velocity Rescaling and the Berendsen Thermostat Lead to the Violation of the Equipartition Theorem

By integrating the classical Newtonian equations of motion, molecular dynamics simulations naturally sample the microcanonical ensemble due to conservation laws. For comparison with experiment, it is often desirable to sample... [ view full abstract ]

By integrating the classical Newtonian equations of motion, molecular dynamics simulations naturally sample the microcanonical ensemble due to conservation laws. For comparison with experiment, it is often desirable to sample constant-temperature ensembles for which thermostats are applied to perturb the Hamiltonian dynamics. In the 1990s, it was found that the simple velocity rescaling and Berendsen thermostat algorithms introduce an artifact termed “the flying ice cube effect,” whereby the equipartition theorem is violated by the draining of kinetic energy from high-frequency modes such as bond stretching into low-frequency modes such as center-of-mass translation.^1,2 The cause for this artifact is not well understood, but violation of the equipartition theorem indicates that the intended thermodynamic ensemble is not being ergodically sampled. It has been questioned whether this artifact appears in other thermostats, particularly those that rely on velocity rescaling like the Bussi thermostat. In addition, since the Gaussian isokinetic thermostat is equivalent to simple velocity rescaling within an order(timestep) error,³ it may be suspected to exhibit the artifact as well. Given the widespread use of these thermostats in MD simulations for equilibration, static and dynamic property studies,⁴ and multiple time step integrators,⁵ more understanding is warranted.

We have conducted simulations with moving harmonic oscillators (diatomic molecules) to elucidate the nature of the flying ice cube effect. We explicitly calculate the partitioning of kinetic energies between translational, rotational, and vibrational degrees of freedom, allowing us to determine which thermostats and conditions lead to the violation of equipartition, as well as the manner and degree to which they do so. We have found that while simple velocity rescaling and the Berendsen thermostat lead to the violation of the equipartition theorem, the Gaussian isokinetic thermostat does not; thus, it appears that the flying ice cube effect is due to the order(timestep) error in the velocity rescaling algorithm leading to a systematically non-ergodic sampling. Meanwhile, the Bussi thermostat does not lead to the violation of the equipartition theorem, but when we force the Bussi thermostat to rescale to kinetic energy distributions with deviations from the canonical distribution, we find that larger deviations lead to larger errors in energy equipartitioning; thus, it appears that the error introduced by velocity rescaling is rectified when velocities are rescaled to the correct distribution.

We also take the opportunity afforded by this discussion of the flying ice cube effect to point out the continued dangers of molecular simulation studies using simple velocity rescaling and the Berendsen thermostat, with the latter continuing to be “the most commonly used algorithm for constant temperature MD of biomolecules”⁶. As a case study, we examine a highly-cited publication which used the Berendsen thermostat and reported a result that can only be replicated when using a flawed thermostat.

¹Lemak and Balabaev, Mol. Simul., 1994, doi.org/10.1080/08927029408021981

²Harvey, Tan and Cheatham, J. Comput. Chem., 1998, doi.org/10.1002/(SICI)1096-987X(199805)19:7<726::AID-JCC4>3.0.CO;2-S

³Nosé, Prog. Theor. Phys. Suppl., 1991, doi.org/10.1143/PTPS.103.1

⁴Minary, Martyna and Tuckerman, J. Chem. Phys., 2003, doi.org/10.1063/1.1534582

⁵Minary, Tuckerman and Martyna, Phys. Rev. Lett., 2004, doi.org/10.1103/PhysRevLett.93.150201

⁶Cooke and Schmidler, J. Chem. Phys., 2008, doi.org/10.1063/1.2989802