CSVR thermostat

One popular strategy to control temperature in NVT MD is to rescale atomic velocities ($\bold{v}_{i}$) at a certain predefined frequency by some factor $\alpha = \sqrt{\bar{K}/K}$ in such a way that the total kinetic energy of the system

$K= \frac{1}{2} \sum\limits_{i=1}^{N} m_i |\bold{v}_{i}|^2,$

is equal to the average kinetic energy corresponding to given temperature:

$\bar{K} = \frac{1}{2}N_f k_B T$

where $N_f$ is the number of degrees of freedom (e.g., $N_f = 3N -3$ in the case of 3D periodic systems) and $N$ is the number of atoms per the simulation cell. Such a method, however, suffers from several problems. First, the ensemble generated is not strictly canonical. Second, rescaling velocities creates discontinuities in trajectories. As a consequence, the method has no conserved quantity that could be used to guide the choice simulation parameters, such as the size of the integration step and also the evaluation of time correlations is problematic. Finally, the trajectories generated via a naive rescaling method often suffer from flying ice-cube problem, i.e., transfer of kinetic energy of a part of the vibrational degrees of freedom into translations and/or rotations, violating thus equipartition principle.

The canonical sampling through velocity rescaling (CSVR) proposed by Bussi et al. removes most of the difficulties of the naïve rescaling approach. Here, the term $\bar{K}$ is replaced by $K_{t}$ obtained for each time step by propagating in time via auxiliary dynamics

$dK = (\bar{K} - K) \frac{dt}{\tau} + 2\sqrt{\frac{K\bar{K}}{N_f}} \frac{dW}{\sqrt{\tau}}$

where $dW$ is a Wiener noise and $\tau$ determined the time scale of the thermostat. Importantly, the auxiliary dynamics generates cannonical distribution for kinetic energy:

$P(K_t) dK_t \propto K_t^{(N_f/2 - 1)} e^{-K_t/k_B T} dK_t$

In practice,