Chapter 4
Temperature control
Context: Most molecular dynamics calculations are carried out in thermal equilibrium. Equilibrium is typically maintained by coupling the molecular calculation to a virtual heat bath, with which it exchanges energy but no particles. This chapter discusses properties of thermal equilibrium and introduces simple algorithms for heat-bath coupling.
Additional resources:
- Chapters 1-3 of Sethna, Entropy, Order Parameters, and Complexity
- Appendix ?? on statistical mechanics.
4.1 Introduction
In order to talk about temperature control, we need to discuss the properties of thermal equilibrium. This is the realm of statistical mechanics or statistical thermodynamics that is discussed in more detail in Chapter ?? and Appendix ??. A key outcome is that the velocity components are distributed according to a Boltzmann distribution. The velocity magnitude is then distributed according to a Maxwell-Boltzmann distribution.
A thermostat implicitly couples the atomistic system under investigation to a much larger heat bath. Because it is much larger, its temperature will not change when energy flows from and to the heat bath. The atomistic system becomes canonical and its statistics follows the canonical ensemble. An ensemble here describes which parameters are constrained, and the canonical ensemble is often also called the NVT-ensemble, because particle number \(N\), volume \(V\) and temperature \(T\) are constrained (fixed). An ideal thermostat guarantees relaxation of the distribution of atomic degrees of freedom to the canonical distribution function (see Chapter ??).
Here, we will start with a mechanistic treatment of thermostats and underpin it with more rigorous theory in Chapter ??. The present chapter teaches the basic concepts required for an implementation of simple thermostatting schemes. Thermostats can be roughly categorized into constraint methods (velocity rescaling and Berendsen), stochastic methods (Andersen, Langevin and dissipative particle dynamics) and extended system methods (Nosé-Hoover). Constraint and extended system methods are deterministic, i.e. they follow the same path when starting from the same initial state. In this chapter we will only discuss the simple constraint methods. We will come back to more advanced methods for temperature control later in these notes.
4.2 Simple themostatting schemes
4.2.1 Velocity rescaling
The crudest (and simplest) form of fixing the temperature in a molecular dynamics simulation to a value of \(T_0\) is by velocity rescaling. Since the instantaneous temperature is \begin {equation} \frac {3}{2} N k_B T = \sum _i \frac {1}{2} mv_i^2, \end {equation} we obtain a temperature of \(T_0\) if we rescale all velocities by \begin {equation} \label {eq:velocity-rescaling} \vec {v}_i \to \lambda \vec {v}_i\ \text { with } \ \lambda =\sqrt {\frac {T_0}{T}} \end {equation} after every time step. This is a very intrusive way of setting the temperature and should not be used in any practical situations, but it is a good illustration of how a simple constraint method works.
4.2.2 Berendsen thermostat
The Berendsen et al. (1984) thermostat uses a damping or acceleration term to control the temperature. The governing equations of motion of the Berendsen thermostat are \begin {equation} m \dot {\vec {v}}_i = \vec {f}_i + \frac {m}{2\tau } \left ( \frac {T_0}{T} -1 \right ) \vec {v}_i \label {eq:berendsen} \end {equation} where \(\tau \) is a relaxation time constant. The factor in front of the velocity is a damping coefficient. The coefficient vanishes for \(T = T_0\), Eq. \eqref{eq:berendsen} then reduces to Newton’s equation of motion. However, it has a positive sign (=speeds up particles) for \(T < T_0\) and has negative sign (=slows down particles) for \(T > T_0\). From Eq. \eqref{eq:berendsen} we can easily derive a differential equation for the evolution of the temperature: \begin {align} 3k_B \frac {dT}{dt} &= \sum _i m \vec {v}_i \cdot \dot {\vec {v}}_i \\ &= \sum _i \left [ \vec {v}_i \cdot \vec {f}_i + \frac {1}{2\tau } \left ( \frac {T_0}{T} - 1 \right ) m v_i^2 \right ] \\ &= - \frac {dE_\text {pot}}{dt} + \frac {3k_B(T_0 - T)}{\tau } \end {align}
This can be written as \begin {equation} \label {eq:berendsen-temperature-evolution} \frac {dT}{dt} = -\frac {T-T_0}{\tau } + S \end {equation} where \(S=-\frac {1}{3k_B} \frac {dE_\text {pot}}{dt}\) is the change of potential energy and constitutes an additional temperature (energy) source.
For \(S = 0\), this equation is solved by \begin {equation} \label {eq:berendsen-explicit-temperature-evolution} T(t) = T_0 + (T_1-T_0) e^{-t/\tau } \end {equation} The temperature relaxes exponentially from the intial value \(T_1\) towards \(T_0\). We directly see that \(\tau \) in Eq. \eqref{eq:berendsen} is indeed the relaxation time constant.
Note that Eq. \eqref{eq:berendsen-explicit-temperature-evolution} suggests an implementation of the Berendsen thermostat in terms of velocity rescaling. During at single time step \(\Delta t \ll \tau \), the temperature should from \(T\) to \(T_0 + (T-T_0)e^{-\Delta t / \tau }\). We can implement this as velocity rescaling, Eq. \eqref{eq:velocity-rescaling}, but with the target temperature \(T_0 + (T - T_0)e^{-\Delta t/\tau }\) instead of \(T_0\), i.e. with a scaling factor \begin {equation} \lambda = \sqrt {\frac {T_0}{T} + \left ( 1- \frac {T_0}{T} \right ) e^{-\frac {\Delta t}{\tau }}} \approx \sqrt {1+ \left ( \frac {T_0}{T} - 1 \right ) \frac {\Delta t}{\tau }} \end {equation} where \(T\) is the current (measure) temperature and \(T_0\) is the target temperature.
A Berendsen thermostat therefore constitutes a gentle way of rescaling velocities. The relaxation time \(\tau \) determines the strength of the coupling between thermal bath and atomistic system. The velocity rescaling limit \(\lambda \to \sqrt {T_0/T}\) is obtained as \(\tau \to \Delta t\). Thermostats should be tuned as weak as possible and as strong as necessary to disturb the system the least while still allowing it to reach the target temperature within the simulation time. There is the additional requirement \(\tau \gg \Delta t\) (where \(\Delta t\) is the time step), otherwise equation Eq. \eqref{eq:berendsen-temperature-evolution} will not be sampled properly numerically. The velocity rescaling thermostat discussed above is bad because it is very strong, but also because it violates \(\tau \gg \Delta t\).
4.3 Equilibrating a molecular simulation
A “happy” molecular dynamics simulation will nicely run at constant temperature. Simulations are only this happy once they are equilibrated and this equilibration implies that the positions \(\{\vec {r}_i\}\) are such that the system resides somewhere near a (potentially local) minimum in the potential energy landscape. When we set up a new simulation, we have to guess a set of \(\{\vec {r}_i\}\) that are often far away from this minimum. (For crystalline solids this guess is simple, since we typically know the crystal structure that we are interested in. For liquids, the guess is more difficult since the overall structure is disordered.) Since the forces \(\{\vec {f}_i\}\) point towards the minimum, the system will evolve in this direction and the potential energy \(E_\text {pot}\) will decrease over time, \(d E_\text {pot}/dt < 0\). Equation \eqref{eq:berendsen-temperature-evolution} tells us, that this leads to an increase in temperature since \(S>0\).
A common problem is that this temperature can be large enough to vaporize the system, i.e. the temperature increases above the vaporization point. The first step in any molecular dynamics simulation is hence to equilibrate the system while avoiding a temperature rise above the point of vaporization (or melting if you are setting up a solid). This can be achieved by running a calculation with a Berendsen thermostat and a strong coupling (i.e. a small \(\tau \)). Once the system has equilibrated, the value of \(\tau \) can be adjusted to a more reasonable relaxation time that does not disturb the calculation too much. Good values for \(\tau \) are between \(1\) ps and \(10\) ps.
Note that if we continuously pump energy into our system, for example because we deform it externally, then Eq. \eqref{eq:berendsen-temperature-evolution} acquires a non-zero source term, \(S > 0\). Assuming \(S\) is constant over time, the final temperature is shifted to \(T_0\,' = T_0 + S\tau \). This temperature offset gets smaller with increasing coupling strength \(1/\tau \).
Bibliography
H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak. Molecular dynamics with coupling to an external bath. J. Chem. Phys., 81(8):3684–3690, 1984. URL https://doi.org/10.1063/1.448118.