Extrapolating molecular dynamics simulations to zero time step and across thermodynamic space

The integration time step is a critical determinant of performance in molecular dynamics simulations, governing the trade-off between speed and fidelity. Although 2 fs remains the standard in atomistic biomolecular simulations, the push for performance has popularized a 4 fs time step with hydrogen mass repartitioning, often combined with multiple time stepping or mass rescaling. However, it is often unclear whether a chosen protocol is overly aggressive, as the apparent numerical stability of a trajectory can mask underlying thermodynamic inaccuracies. Increasing the time step will exacerbate systematic discretization errors, inherent to all numerical integration algorithms. In the widely used Verlet family of integrators, these errors manifest as O(Δt2) deviations in thermodynamic observables such as potential energy and volume, and for common Langevin splitting schemes, even temperature. We demonstrate that these deviations follow a simple, linear thermodynamic model, allowing for their rigorous removal by extrapolation to the zero time step limit. In turn, the time-step dependence provides us with estimates of the system heat capacity, compressibility, and thermal expansion coefficient. This framework allows us to construct consistent probability distributions of energy and volume across thermodynamic states, effectively recovering Boltzmann-consistent statistics at a target condition independent of time step. These considerations are particularly important for enhanced sampling methods, such as replica exchange and umbrella sampling, which rely on rigorous Boltzmann sampling and require accurate energies and temperatures for valid replica exchange probabilities and statistical reweighting.