Canonical sampling through velocity rescaling

This paper addresses a central problem in molecular dynamics (MD): how to generate trajectories that sample the canonical (constant number of particles, volume, and temperature; NVT) distribution efficiently and accurately, while retaining good numerical stability and (ideally) preserving dynamical properties. The authors focus on thermostats—algorithms that control temperature by modifying particle velocities—and ask whether a simple, velocity-based stochastic scheme can be constructed that (i) samples the correct canonical ensemble, (ii) admits a conserved-like diagnostic quantity to assess sampling accuracy under finite time steps, and (iii) performs well across a range of thermostat parameters without excessive tuning.

The significance is twofold. First, many widely used thermostats (e.g., Berendsen, Nosé-Hoover) either do not correspond to a well-defined ensemble in general or can suffer from ergodicity issues, meaning that time averages may not match ensemble averages for certain systems. Second, even when a thermostat is formally correct, finite integration time steps introduce errors; without a robust diagnostic, it is difficult to know whether the simulation is truly sampling the intended distribution. The authors position their method as an extension of velocity-rescaling ideas: rather than deterministically forcing the kinetic energy to a target value (as in Berendsen-like schemes) or randomly reassigning individual particle velocities (as in Andersen), they rescale all particle velocities by a single random factor chosen so that the kinetic energy follows the correct canonical distribution.

Methodologically, the paper develops two related thermostat formulations. The simplest version (Section II A) rescales all particle velocities by a factor $\alpha =\sqrt{K_t/K}$, where $K$ is the instantaneous kinetic energy and $K_t$ is drawn from the canonical equilibrium distribution for kinetic energy. For a system with $N_f$ degrees of freedom and inverse temperature $\beta$, the target distribution is proportional to $\bar{P}(K_t)d{K}_t\propto K_t^{\frac{N_f}{2}-1}{e}^{-\beta K_t}d{K}_t$. Between rescaling events, the system evolves under Hamilton’s equations using an ergodic microcanonical dynamics assumption. Because Hamiltonian evolution preserves the microcanonical measure and the rescaling step is constructed to preserve the canonical distribution, the combined process samples the canonical ensemble (with a caveat that the center-of-mass momentum and angular momentum slices are preserved; the authors argue this effect is negligible for their purposes).

The authors then introduce a smoother, more general scheme (Section II B) to reduce fast fluctuations induced by frequent full rescaling. The key idea is to evolve the kinetic energy $K$ through an auxiliary stochastic differential equation over a time interval corresponding to one MD time step, and then apply a rescaling consistent with the updated kinetic energy. The auxiliary dynamics is defined so that the kinetic energy marginal remains canonical. The most general one-dimensional Fokker–Planck form with zero-current stationary solution is given, leading to a stochastic differential equation for $K$. They choose a specific diffusion function $D(K)=\frac{2K\bar{K}}{N_f\tau }$, introducing a thermostat relaxation time $\tau$. This yields the explicit auxiliary dynamics $$dK=\frac{\bar{K}-K}{\tau }dt+2\sqrt{\frac{K\bar{K}}{N_f}}dW\sqrt{\tau },$$ where $dW$ is a Wiener noise increment. In the limit $\tau \to 0$, the stochastic evolution thermalizes instantly and the method reduces to the simpler stochastic velocity-rescaling thermostat. In the limit $\tau \to \infty$, Hamiltonian dynamics is recovered. The algorithm is designed to preserve bond constraints automatically when a single global scaling factor is applied to all velocities, and it conserves total linear momentum (and angular momentum for non-periodic systems).

A major methodological contribution is the introduction of an “effective energy” $\tilde{H}$ used to diagnose sampling accuracy under finite time steps. For discrete-time dynamics targeting a distribution $\bar{P}(x)$, they define weights $w_i$ based on the ratio of transition probabilities for forward and momentum-reversed moves. For canonical targets $\bar{P}(x)\propto e^{-\beta H(x)}$, this leads to an effective energy ${\tilde{H}}_i=-\frac{1}{\beta }\ln w_i$. In their thermostat, the velocity-rescaling step is constructed to satisfy detailed balance exactly (when the auxiliary kinetic-energy equation is solved exactly), so changes in $\tilde{H}$ arise primarily from the Hamiltonian integration step. Therefore, if the integrator is accurate, $\tilde{H}$ should remain nearly constant, analogous to energy conservation in NVE simulations. This diagnostic also supports reweighting or hybrid Monte Carlo correction strategies.

For numerical evaluation, the authors test the method on Lennard-Jones and water models in both solid and liquid phases. Lennard-Jones: 256 atoms in a cubic box, using an argon parameterization; simulations in fcc solid at 20 K and liquid at 120 K. Water: TIP4P with rigid molecules, particle-mesh Ewald for electrostatics; simulations in cells of 360 molecules for ice Ih and liquid water. They compute energy fluctuations (kinetic and potential) over 1 ns trajectories for a range of $\tau$ values spanning three orders of magnitude, and they compare to Nosé–Hoover and Berendsen thermostats. They also assess dynamical properties: vibrational density of states for hydrogen atoms in ice Ih via Fourier transform of velocity autocorrelation, and self-diffusion coefficient for liquid water via Einstein’s relation from mean-square displacement.

Key findings are reported with specific quantitative comparisons. For the time-step diagnostic, they show that for Lennard-Jones at 120 K with $\tau =0.1$ ps, using $\Delta t=5$ fs yields good conservation of $\tilde{H}$ with fluctuations around $0.3$ kJ/mol, while the root-mean-square fluctuations in the conventional energy $H$ are on the order of $16$ kJ/mol. When the time step is increased to $\Delta t=40$ fs, $\tilde{H}$ exhibits a systematic drift, indicating inaccurate sampling even though the thermostat keeps the simulation stable (in contrast to NVE where the system can “explode”).

Regarding ensemble correctness via energy fluctuations, the authors plot square fluctuations of kinetic energy $\Delta K^2$ and potential energy $\Delta U^2$, normalized by the ideal-gas kinetic-energy fluctuation $N_f{k}_B^2{T}^2/2$. For Lennard-Jones, they find that their thermostat reproduces the correct canonical fluctuation behavior across the full $\tau$ range for both liquid and solid. In contrast, Berendsen’s thermostat shows strong dependence on $\tau$: as $\tau \to 0$, it approaches the isokinetic ensemble, causing kinetic-energy fluctuations to vanish while configurational fluctuations approach the canonical limit. Nosé–Hoover matches their thermostat for the liquid over a wide $\tau$ range but fails for the solid when $\tau >0.2$ ps due to ergodicity problems (leading to poor sampling). For water and ice, they report that Nosé–Hoover becomes inefficient for ice when $\tau >0.2$ ps and for water when $\tau >2$ ps, whereas their thermostat remains fairly independent of $\tau$.

For dynamical properties, they show that the vibrational density of states of ice Ih is preserved when using their thermostat. For $\tau =2$ ps, there is no shift of the main peak frequencies relative to NVE, and intensity changes are within numerical errors. With very small $\tau =0.002$ ps, the low-frequency translational broad peak is slightly affected, but no fictitious peaks appear.

They quantify diffusion in liquid water (TIP4P) at 300 K. Table I reports self-diffusion coefficients $D$ (in units of $10^{-5}{\mathrm{cm}}^2/\mathrm{s}$) for different $\tau$: $\tau =0.002$ ps gives $3.63\pm 0.01$, $\tau =0.02$ ps gives $3.44\pm 0.06$, $\tau =0.2$ ps gives $3.51\pm 0.05$, and $\tau =2$ ps gives $3.53\pm 0.01$. The NVE reference value is $3.47\pm 0.03$. The authors conclude that their thermostat yields diffusion consistent with NVE and literature, with only marginal deviation at very small $\tau$.

Limitations include the authors’ reliance on an ergodicity assumption for Hamiltonian dynamics in the microcanonical ensemble; if ergodicity fails, canonical sampling may not hold. They also acknowledge (and partially neglect) that rescaling does not change center-of-mass momentum and, for non-periodic systems, angular momentum, so the algorithm samples a slice of the canonical ensemble unless those degrees of freedom are already appropriately distributed. Additionally, the effective-energy diagnostic depends on using the exact propagator for the auxiliary kinetic-energy equation; if the auxiliary dynamics is not solved exactly, detailed balance and the diagnostic’s interpretation could degrade. Finally, as with any MD thermostat, extreme parameter choices (very small $\tau$) can perturb dynamics.

Practically, the results matter for anyone running MD simulations where correct canonical sampling and reliable numerical integration are required—especially for small systems, systems with slow modes, or solids where ergodicity problems are common. The method is presented as easy to implement (global velocity rescaling with a stochastic factor), stable, and largely insensitive to the thermostat parameter $\tau$. The effective energy $\tilde{H}$ provides a concrete way to assess whether finite time-step errors are compromising ensemble fidelity, and it can guide time-step selection or enable reweighting/hybrid Monte Carlo corrections.

Overall, the paper’s core contribution is a thermostat that combines the simplicity of velocity rescaling with a theoretically justified stochastic kinetic-energy evolution that samples the canonical ensemble, plus a conserved-like diagnostic quantity to quantify integration-induced ensemble errors. The authors show that their velocity-rescaling thermostat samples the canonical NVT ensemble while remaining robust across thermostat parameters and offering an effective-energy measure of sampling accuracy.