Velocity Scaling: Isokinetics and the Berendsen Thermostat

First, velocity scaling schemes do not strictly follow the canonical ensemble, though in practice, the amount they deviate from canonical is quite small. (This can be measured by comparing the velocity distribution function with a Gaussian.) It is relatively easy to implement the second class of schemes, because they can be ``dropped'' in to existing codes using almost any integrator. However, they suffer the drawback that they are not time-reversible or deterministic, properties that become important in some advanced MD techniques. The third class are slightly more difficult to implement, but do not suffer from such drawbacks as time-irreversibility.

We have in effect already encountered simple velocity scaling in mdlj.c, in the initialization function. Here, particle velocities are chosen randomly from $[-0.5,0.5]$ and the rescaled to result in a desired temperature given by the relation:

$$\frac{3}{2}Nk_BT = \frac{1}{2}\sum_im_iv_i^2$$ (176)

We could, if we wanted to, turn this into a dynamic scheme for continually keeping the velocities scaled such that the total kinetic energy is constant. We can measure the instaneous temperature immediately after a velocity update, and call it $T_i$. Eq. 176 indicates that if we scale velocities by a constant $\lambda$, where

$$\lambda = \sqrt{\left(T/T_i\right)}$$ (177)

we will be left with a system at temperature $T$. Velocity scaling to maintain constant $\mathscr{K}$ is called an isokinetic thermostat. Such a thermostat cannot be used to conduct a simulation in the canonical ensemble, but is perfectly fine to use in a warmup or initialization phase. We could perform velocity rescaling at every step, or only every few steps. As a suggested exercise, modify mdlj.c to perform velocity scaling to a user-specified setpoint temperature every $m$ time steps, where $m$ is a user-defined interval between velocity scaling events. Begin with system at $T$ = 1.0, and command it to jump to $T$ = 3.0 after 1,000 steps. How does the system behave, and is it sensitive to your choice of $m$?

Another popular velocity scaling thermostat is that of Berendsen [10]. Here, the scale factor is given by

$$\lambda = \left[1+\frac{\Delta t}{\tau_T}\left(\frac{T}{T_0} - 1\right)\right]^{\frac{1}{2}}$$ (178)

Here, $T_0$ is the setpoint temperature, $\Delta t$ is the integration time step, and $\tau _T$ is a constant called the ``rise time'' of the thermostat. It describes the strength of the coupling of the system to a hypothetical heat bath. The larger $\tau _T$, the weaker the coupling; in other words, the larger $\tau _T$, the longer it takes to achieve a given $T_0$ after an instantaneous change from some previous $T_0$. The code mdlj_ber.c implements the Berendsen thermostat. As a brief exercise, you can experiment with this code to get a feeling for how various values of the rise time affect the response of the system when the setpoint temperature is changed instantaneously from 1.0 to 2.0. Below is a lin-log plot of just such an experiment with $N$ = 256 particles at a density of 0.5. Each curve corresponds to a different value of $\tau _T$, and they increase by factors of 10. The corresponding time at which the setpoint $T$ is reached is also seen to increase by the same factor.

Instantaneous temperature, $T$, vs. time in an MD simulation of 256 particles at a density of 0.5, with temperature controlled by the Berendsen thermostat [10], for various values of the thermostat ``rise time,'' $\tau _T$.

Though relatively simple, velocity scaling thermostats are not recommended for use in production MD runs because they do not strictly conform to the canonical ensemble.