9
$\begingroup$

With reference to the GROMACS manual 2023.2, the discussion of the Berendsen thermostat (p. 362) indicates that the time constant $\tau$ which governs the steepness of the relaxation toward a target temperature (eqn. 5.37: $\frac{\mathrm{d}T}{\mathrm{d}t} = \frac{(T_0 - T)}{\tau}$) is not equal to to the parameter $\tau_T$ found in the description of the re-scaling factor (eqn. 5.38): $$\lambda = \left[ 1 + \frac{n_{\mathrm{TC}}\Delta t}{\tau_T}\left(\frac{T_0}{T(t - \frac{1}{2} \Delta t)} - 1\right) \right] ^{\frac{1}{2}}.$$

They are related via eqn 5.39: $\tau = 2C_V\tau_t / N_{df}k$.

The reason cited is: "The reason that $\tau \neq \tau_T$ is that the kinetic energy change caused by scaling the velocities is partly redistributed between kinetic and potential energy and hence the change in temperature is less than the scaling energy."

I don't understand this explanation. Where does the redistribution of potential energies occur? Thanks.

$\endgroup$
2
  • 1
    $\begingroup$ Can you expand your question please? What is the re-scaling factor? $\endgroup$
    – MSwart
    Oct 19, 2023 at 10:56
  • $\begingroup$ The re-scaling factor is applies every $n_{\mathrm{TC}}$ steps to the velocities of all particles. This attempts to restrain the kinetic energies, and thus temperatures, to match the MB distribution of speeds. My confusion: scaling velocities seems to imply changes to the kinetic energy only. I don't believe the original paper discusses a second time-constant. (Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. (1984).) @MSwart $\endgroup$ Oct 19, 2023 at 20:00

2 Answers 2

4
$\begingroup$

The cited formula for $\lambda$ is a little bit non-standard compared to what other codes do because of the $n_\mathrm{TC}$ parameter. E.g., inside LAMMPS, $\lambda$ is obtained as $$ \lambda = \left[ 1 + \frac{\Delta t}{\tau} \left( \frac{T_0}{T} - 1\right) \right]^{1/2} $$ which more or less reduces to the GROMACS formula when setting $n_\mathrm{TC} = 1$, i.e., applied at every time step (although LAMMPS uses velocities at the end of the velocity verlet step, i.e., using $T(t + \Delta t)$, rather than the halfway point of the integration).

Now, regarding relaxation time. Suppose you've run NVE dynamics for some time, and reached equilibrium. Now, you're applying the thermostat exactly once, setting $\tau = \Delta t$. That's equivalent to scaling the kinetic energy by $T_0/T$, so you would assume you've achieved perfect thermalization within a single step.

However, the instantaneous coordinates of the system correspond to a state of temperature $T$. For example, if $T$ is high, bonds will on average be more stretched, and molecule geometries will in general be further from their minimum energy configurations. If you would now scale down velocities so that the kinetic energy corresponds to a lower $T_0$, you still have lots of excess energy stored in the potential energy component of the system. If you continue to run NVE MD, you'll notice that the average potential energy decreases, while the kinetic energy increases again. So in the end, the average temperature will have a value between the old $T$ and the target $T_0$.

On a broader note, this means that $\tau$ is really just a parameter, rather than truly expressing the relaxation time scale. A smaller $\tau$ will speed up equilibration, but the potential energy (and system configuration) will lag the manipulations of the kinetic energy. So setting $\tau = 10$ ps does not mean the system is fully equilibrated after 10 ps.

Finally, it's probably good to emphasize that nowadays there are not many reasons to still use a Berendsen thermostat for any production calculations except when trying to reproduce old simulations. This is because the Berendsen thermostat does not generate a correct canonical ensemble. All big simulations codes implement several algorithms that do, such as Nosé-Hoover chains, the Bussi-Donadio-Parrinello rescaling thermostat, and several types of (Generalized) Langevin thermostats (see also here). Even if correct sampling is not required, many of these newer methods are also more efficient. A wider discussion as to which approach is best-suited for a given problem is beyond the scope of the current question.

$\endgroup$
1
$\begingroup$

I can somewhat explain it but I agree that the literature is very frustrating to read. Basconi and Shirts have the modified time constant in their 2013 paper1 but only cite Berendsen's 1984 paper, which does not seem to make the same claim. I have not located the claim in any paper in between -- but it makes sense to me, after quite a lot of thought.


Imagine a canonical ensemble trajectory for some system at some temperature $T$, which has average total energy $\bar{E}$, average kinetic energy $\bar{K}$ and average potential energy $\bar{U}$. (In the canonical ensemble, heat is always exchanged with a bath, and all three values have varying instantaneous values.)

Now imagine a microcanonical ensemble trajectory for the same system, with a conserved energy equal to the canonical average: $E(t) = \bar{E} = $ some constant. To the extent that thermodynamic equivalence is valid, the microcanonical averages of the kinetic and potential energies $K$ and $U$ must be same as their canonical averages. But their instantaneous fluctuations and correlations will be different. In particular, the trajectory will always be in one of two conditions:

  1. Instantaneously $K > \bar{K}$ and $U < \bar{U}$ -- the system is kinetically "hot" and configurationally "cool", or

  2. Instantaneously $K < \bar{K}$ and $U > \bar{U}$ -- the system is kinetically "cool" and configurationally "hot".

Finally, consider a trajectory under Berendsen thermostatting, followed for a short enough time that it remains similar to a corresponding microcanonical trajectory. Suppose case 1 applies, and the system is kinetically "cool". Then the Berendsen thermostat will pump in kinetic energy -- at the same time that the system is configurationally "hot" and due to cool down! Conversely, when the system is kinetically "hot", the Berendsen thermostat pumps out energy, fighting a "cool" configuration.

And that's why the Berendsen time constant encompasses the entire systemic heat capacity, not just the kinetic.


As far as I can tell this is unique to the Berendsen thermostat because of its first-order coupling to the kinetic temperature -- that is, it adds to the kinetic energy's first derivative: $\dot{K} \mapsto \dot{K} - \alpha^2 (K - K_{\text{goal}})$. Other thermostats happen to avoid this in some way. The isokinetic thermostats (that simply force the kinetic energy to be exactly "right" every step) do not have a "timescale" to speak of. Stochastic thermostats have a random term that decorrelates the thermostat energy flow from the kinetic temperature.

And the Nose-Hoover thermostats add second-order coupling instead of first-order coupling ($\ddot{K} \mapsto \ddot{K} - \alpha^2 (K-K_{\text{goal}})$), so their timescale is really controlled by the "kinetic" heat capacity ($N_{\text{dof}}k_BT/2$) rather than the overall heat capacity. The Nose-Hoover thermostats effectively add "auxiliary variables" to the system, whose distribution functions really are appropriately harmonic2. So as far as I can tell, it's not just the thermostat's energy eventually "flowing" into the configurational temperature -- it's about how those fluctuations correlate.


1 Basconi and Shirts, JCTC 2013, 9, 7, 2887–2899

2 Bernhardt and Petersen, PCCP 2022, 24, 6383-6392

$\endgroup$
1
  • $\begingroup$ ah, that's indeed a very nice perspective; Berendsen only working on the kinetic energy, and ignoring the coupling with configurational degrees of freedom. Langevin-type methods (and also Bussi-Donadio-Parrinello) are stochastic, also breaking this first order coupling in another way $\endgroup$ Oct 24, 2023 at 12:37

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .