How Thermostats Influence Dynamics Across Time Scales: A Systematic Study from Fast Motions to Slow Transitions

TL;DR

Thermostatting can bias dynamical properties in molecular dynamics simulations. The authors systematically benchmark NVE, deterministic thermostats, velocity-rescale, and stochastic Langevin thermostats across liquids, examining time-correlation functions from vibrational to slow conformational scales, including diffusion coefficient $D$, shear viscosity $\eta$, vibrational density of states (vDOS), and Markov state model (MSM) kinetics. They find that deterministic thermostats and the velocity-rescale scheme reproduce NVE reference data across observables, while strongly coupled stochastic thermostats distort diffusion, viscosity, and MSM timescales; moderate stochastic coupling with coupling time around $\tau_T \sim 1~\mathrm{ps}$ restores near-NVE behavior while still sampling the canonical ensemble. The results yield practical guidelines: for accurate dynamical properties or reliable MSMs, choose thermostat schemes and coupling strengths that balance energy fluctuations with minimal dynamical distortion, typically 1–10 ps for stochastic thermostats.

Abstract

Reliable dynamical properties from molecular dynamics simulations require careful control of thermostatting artifacts. We systematically assess how NVE, deterministic thermostats, velocity-rescale dynamics, and stochastic Langevin-type thermostats affect time-correlation functions across liquids of varying complexity. The analysis spans vibrational spectra, velocity and pressure autocorrelations, diffusion coefficients, shear viscosities, and Markov state models. Deterministic thermostats and velocity-rescale dynamics closely reproduce NVE reference data over all observables. In contrast, strongly coupled stochastic thermostats (tau less 1 ps) systematically distort dynamical properties. By constrast, moderate stochastic coupling (tau eq. 1 ps) restores near-NVE behavior while maintaining canonical sampling. Our results provide practical guidelines for selecting thermostat schemes when accurate dynamical properties or Markov models are required.

Figures (19)

• Figure 1: Energy distribution of TIP3P water in NVE and NVT ensemble with different thermostats with a simulation time of 60 ns. The coupling time $\tau_T$ is 1 ps for Berendsen and velocity rescale and 2 ps for GROMACS stochastic dynamics GSD. The chain length for Nose-Hoover is 1 with a coupling time of 4 ps. NVT is calculated using a normal distribution and eq. \ref{['eq:sigma']} with a Cv of 4.14 kJ/kgKMao2012
• Figure 2: Self diffusion of water. a) VACF of TIP3P water for different thermostat coupling times $\tau_T$ of VR thermostat and SD thermostat. b) Apparent self-diffusion coefficient of TIP3P and xTB water for various thermostats and thermostat coupling times $\tau_T$. Error bars show standard deviations obtained from block averaging.
• Figure 3: Shear viscosity. PACF of TIP3P water (a) and anilin (b) for different thermostat coupling times $\tau_T$ of VR thermostat and SD thermostat. c)Long tail of the PACF of anilin. d) Shear viscosities of pentane, TIP3P water, anilin and glycerol for different thermostats and coupling times $\tau_T$. The error is the standard deviation obtained from block averaging. Note the logarithmic scale of the viscosity-axis.
• Figure 4: a) Atomic vibrational autocorrelation function (VACF) of liquid water obtained from AIMD simulations using various coupling times $\tau_T$ for both velocity rescaling (VR) and stochastic dynamics (SD, here: OBABO). A logarithmic scale is used for lag times greater than 0.1 ps. The shaded areas represent the envelopes of ideal damped oscillators with decay times of $\tau_T=0.01$ and $0.1$ ps, respectively. b) Vibrational density of states (vDOS), computed as the Fourier transform of the VACF.
• Figure 5: Markov state model for a) pentane. b) Stationary distribution at lag time $\tau=20.1$ ps. c) Slowest dynamical processes ($\mathbf{l}_1$, $\mathbf{l}_2$) for GSD with different coupling times $\tau_T=0.01$ ps, $\tau_T=1$ ps and $\tau_T=10$ ps, with shared color bar. d) MSM correlation matrix elements $C_{5,5}$ as function of the lag time ($0.1-59.1$ ps). e) MSM implied timescales as mean with standard deviation over converged lag time ($20.1-59.1$ ps) for $\lambda_1$ and $\lambda_2$.