A common, avoidable source of error in molecular dynamics integrators

In constrained molecular dynamics simulations using some of the most popular molecular dynamics codes, calculation of the velocities of constrained particles is based solely on the differences in particle positions during two successive time steps. This creates a numerical instability that the authors’ show to be signicant in a typical single-precision floating-point simulation. They describe a simple modification that eliminates this source of instability and demonstrate that this change substantially reduces the energy drift of a sample single-precision NVE simulation.