The Nosé–Hoover thermostat refers to a class of deterministic algorithms for molecular dynamics (MD) simulations at constant temperature. In such simulations, the system of interest interacts with a heat bath, and the dynamics of both the system and the reservoir are calculated over time. The thermostat is named after Shūichi Nosé, who in 1984 proposed the first method for simulations of the coupled system and heat bath,¹^,² and William G. Hoover—who in the following year improved Nosé’s method and presented the coupled dynamics in the form commonly used today.³ A remarkable feature of the Nosé–Hoover equations is that only a single fictitious particle is used to describe the heat bath. Consequently, the Nosé–Hoover thermostat is computationally inexpensive, and is widely used in constant temperature molecular dynamics simulations.

Why do we need a thermostat?

A “traditional” molecular dynamics simulation will sample the microcanonical ensemble, in which the total energy is fixed, rather than the canonical ensemble, in which the temperature is fixed. To see why, we consider a simulation involving \(N\) particles occupying a volume \(V\), where the position and momentum of the \(\alpha^{\textrm{th}}\) particle are denoted \(\boldsymbol{q}^\alpha\) and \(\boldsymbol{p}^\alpha\), respectively. The (micro)state of a classical system is captured by the set of all particle positions and momenta, for which we introduce the shorthand

\[\boldsymbol{\Gamma} \, := \, \Big(\mkern1mu \boldsymbol{q}^1, \, \boldsymbol{q}^2, \, \ldots \, \boldsymbol{q}^N, \ \boldsymbol{p}^1, \, \boldsymbol{p}^2, \, \ldots \, \boldsymbol{p}^N \mkern1mu\Big) ~.\]

If the system is isolated and governed by a Hamiltonian \(\mathcal{H} = \mathcal{H} (\boldsymbol{\Gamma})\), then the positions and momenta evolve in time according to Hamilton’s equations—written as

\[\boldsymbol{\dot{q}}^\alpha \, = \, \dfrac{\partial \mathcal{H}}{\partial \boldsymbol{p}^\alpha} \qquad\quad \text{and} \qquad\quad \boldsymbol{\dot{p}}^\alpha \, = \, - \, \dfrac{\partial \mathcal{H}}{\partial \boldsymbol{q}^\alpha} ~,\]

where a ‘dot’ accent denotes a time derivative. In this way, the system traces out a trajectory \(\boldsymbol{\Gamma} (t)\) in phase space. By taking the total time derivative of the Hamiltonian and then substituting the equations of motion, we find

\[\dot{\mathcal{H}} \ = \, \sum_{\alpha = 1}^N \bigg( \boldsymbol{\dot{p}}^\alpha \boldsymbol{\cdot} \mkern1mu \dfrac{\partial \mathcal{H}}{\partial \boldsymbol{p}^\alpha} \, + \, \boldsymbol{\dot{q}}^\alpha \boldsymbol{\cdot} \mkern1mu \dfrac{\partial \mathcal{H}}{\partial \boldsymbol{q}^\alpha} \bigg) \, = \, \sum_{\alpha = 1}^N \bigg( \boldsymbol{\dot{p}}^\alpha \boldsymbol{\cdot} \boldsymbol{\dot{q}}^\alpha \, - \, \boldsymbol{\dot{q}}^\alpha \boldsymbol{\cdot} \boldsymbol{\dot{p}}^\alpha \bigg) \, = \ 0 ~.\]

Accordingly, the value of the Hamiltonian—namely, the total (internal) energy \(E\)—is a constant of motion. We have thus shown that an MD simulation will sample the microcanonical ensemble, in which \(N\), \(V\), and \(E\) are constant, if particles evolve in time according to Hamilton’s equations.

The development of microcanonical MD simulations marks a significant achievement in statistical physics. However, the results of such simulations often cannot directly be compared to experimental measurements—in which the temperature \(T\), rather than the internal energy \(E\), is the relevant control parameter. There were consequently many efforts to formulate algorithms for constant temperature MD simulations, one of which is the Nosé–Hoover thermostat.

The development by Nosé

For an MD simulation to sample the canonical ensemble, the system of interest must be perturbed such that its total energy does not remain constant. Nosé’s great insight was to construct an extended, fictitious system involving one additional degree of freedom that represented a heat bath. If the fictitious system is constructed appropriately, then microcanonical simulations of the fictitious system lead to the real system sampling the canonical ensemble. In this fashion, one can simulate Hamilton’s equations in the extended system and guarantee the real system is at the desired temperature. In what follows, we present the main ideas of Nosé’s developments, albeit in a slightly different form from the original formulation.¹^,²

We begin by constructing a fictitious system, where all associated quantities are denoted with a ‘hat’ accent. At present, we do not know what these quantities represent, and we will determine their physical meaning subsequently. It is most natural to begin with the fictitious Lagrangian \(\hat{\mathcal{L}}\), which is a function of \(\{ \boldsymbol{\hat{q}}^\alpha \}\), \(\{ \mathrm{d} \boldsymbol{\hat{q}}^\alpha / \mathrm{d} \hat{t} \mkern1mu \}\), \(\hat{s}\), and \(\mathrm{d} \hat{s} / \mathrm{d} \hat{t}\), and is given by

\[\hat{\mathcal{L}} \, = \, \sum_{\alpha = 1}^N \dfrac{m \hat{s}^2}{2} \, \dfrac{\mathrm{d} \boldsymbol{\hat{q}}^\alpha}{\mathrm{d} \hat{t}} \boldsymbol{\cdot} \dfrac{\mathrm{d} \boldsymbol{\hat{q}}^\alpha}{\mathrm{d} \hat{t}} \, - \, U \big( \{\boldsymbol{\hat{q}}^\alpha\} \big) \, + \, \dfrac{Q}{2} \bigg( \dfrac{\mathrm{d} \hat{s}}{\mathrm{d} \hat{t}} \bigg)^2 - \, g \, k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \, \ln \hat{s} ~.\]

Here \(g\) is a yet-to-be-specified dimensionless constant, \(k_{\mathrm{B}}\) is Boltzmann’s constant, \(T_{\mathrm{ext}}\) is the temperature we seek to impose externally, and \(U\) is the potential energy of the system. The quantity \(\hat{s}\) is dimensionless, implying the corresponding inertial term \(Q\) has dimensions of energy \(\cdot\) time². However, since \(\hat{s}\) is thought of as relating to a heat bath, \(Q\) is often called the thermal mass or thermal inertia. With the Lagrangian, the conjugate momenta are determined according to the standard procedure:

\[\boldsymbol{\hat{p}}^\alpha \, = \, \dfrac{\partial \hat{\mathcal{L}}}{\partial ( \mathrm{d} \boldsymbol{\hat{q}}^\alpha / \mathrm{d} \hat{t} \mkern1mu) } \, = \, m \hat{s}^2 \, \dfrac{\mathrm{d} \boldsymbol{\hat{q}}^\alpha}{\mathrm{d} \hat{t}} \qquad\quad \text{and} \qquad\quad \hat{p}^{}_{\! s} \, = \, \dfrac{\partial \hat{\mathcal{L}}}{\partial ( \mathrm{d} \hat{s} / \mathrm{d} \hat{t} \mkern1mu) } \, = \, Q \, \dfrac{\mathrm{d} \hat{s}}{\mathrm{d} t} ~.\]

We can then construct the corresponding fictitious Hamiltonian as the Legendre transform of the Lagrangian—expressed as

\[\hat{\mathcal{H}} \big( \{ \boldsymbol{\hat{q}}^\alpha \}, \, \{ \boldsymbol{\hat{p}}^\alpha \}, \, \hat{s}, \, \hat{p}^{}_{\! s} \big) \, = \, \sum_{\alpha = 1}^N \, \boldsymbol{\hat{p}}^\alpha \boldsymbol{\cdot} \dfrac{\mathrm{d} \boldsymbol{\hat{q}}^\alpha}{\mathrm{d} \hat{t}} \, + \, \hat{p}^{}_{\! s} \, \dfrac{\mathrm{d} \hat{s}}{\mathrm{d} \hat{t}} \, - \, \hat{\mathcal{L}} ~,\]

from which we find

\[\hat{\mathcal{H}} \, = \, \sum_{\alpha = 1}^N \dfrac{\boldsymbol{\hat{p}}^\alpha \boldsymbol{\cdot} \boldsymbol{\hat{p}}^\alpha}{2 m \hat{s}^2} \, + \, U \big( \{\boldsymbol{\hat{q}}^\alpha\} \big) \, + \, \dfrac{\hat{p}^2_{\! s}}{2 Q} \, + \, g \, k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \, \ln \hat{s} ~.\]

At this point, two tasks remain. We seek to interpret all of the fictitious quantities and relate them to our physical system, and then show that microcanonical simulations of the extended system (as governed by \(\hat{\mathcal{H}}\)) cause the physical system to sample the canonical distribution. In the process of doing so, the parameter \(g\) will be determined.

The physical interpretation of fictitious quantities

It is immediately clear from the fictitious Lagrangian that we need to interpret \(\hat{s}\), as if \(\hat{s} = 1\) then \(\hat{\mathcal{L}}\) is the usual physical Lagrangian governing an \(N\)-particle system at constant energy. To this end, we examine how the fictitious Lagrangian is modified under the choice of scaling

\[\hat{s} \, \rightarrow \, b \mkern1mu \hat{s} ~, \qquad \hat{t} \, \rightarrow \, b^m \mkern1mu \hat{t} ~, \qquad \text{and} \qquad \hat{\boldsymbol{q}}^\alpha \, \rightarrow \, b^n \mkern1mu \hat{\boldsymbol{q}}^\alpha ~,\]

where the exponents \(m\) and \(n\) will be chosen such that the Lagrangian is invariant to the value of the dimensionless scale factor \(b\). Under such a change of variables, the fictitious Lagrangian transforms to

\[\hat{\mathcal{L}} \, \rightarrow \, b^{2 - 2m + 2n} \mkern1mu \sum_{\alpha = 1}^N \dfrac{m \hat{s}^2}{2} \mkern1mu \dfrac{\mathrm{d} \boldsymbol{\hat{q}}^\alpha}{\mathrm{d} \hat{t}} \! \boldsymbol{\cdot} \! \dfrac{\mathrm{d} \boldsymbol{\hat{q}}^\alpha}{\mathrm{d} \hat{t}} \, - \, U \big( \{b^n \mkern1mu \boldsymbol{\hat{q}}^\alpha\} \big) \, + \, b^{2 - 2m} \mkern1mu \dfrac{Q}{2} \mkern-1mu \bigg( \dfrac{\mathrm{d} \hat{s}}{\mathrm{d} \hat{t}} \bigg)^2 - \, g \, k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \mkern1mu \ln (b \mkern1mu \hat{s}) ~.\]

If we choose \(m = 1\) and \(n = 0\), then the Lagrangian is shifted by the constant factor \(- g \mkern1mu k_{\mathrm{B}} T_{\mathrm{ext}} \mkern1mu \ln b\), which does not affect the dynamics. Since \(\boldsymbol{\hat{q}}^\alpha\) is unaffected by the scaling, we interpret it as the real position: \(\boldsymbol{\hat{q}}^\alpha = \boldsymbol{q}^\alpha\). We also observe that \(\hat{s} \rightarrow b \mkern1mu \hat{s}\) and \(\hat{t} \rightarrow b \mkern1mu \hat{t}\), for which the ratio \(\hat{t} / \hat{s}\) is invariant to the scaling. Since \(\hat{t}\) corresponds to the real time when \(\hat{s} = 1\), we interpret \(\hat{t}\) as a scaled time, with \(\hat{s}\) the (continuously changing) value of the scaling. Since MD simulations of the fictitious system will generate \(\hat{s}\) as a function of \(\hat{t}\), we can convert between a fictitious time \(\hat{\tau}\) and real time \(\tau\) according to

\[\mathrm{d} t \, = \, \dfrac{\mathrm{d} \hat{t}}{\hat{s}} \qquad \text{and} \qquad \tau \, = \, \int_0^{\hat{\tau}} \dfrac{\mathrm{d} \hat{t}}{\hat{s} (\hat{t})} ~.\]

Sampling the canonical ensemble

Our final task is to show how solving Hamilton’s equations in the extended, fictitious system sets the temperature to a constant value in the real system. We do so by examining the partition function, which describes how often a system can be found in a particular microstate.

In the microcanonical ensemble, the partition function \(W (E)\) is simply the number of microstates in which the system has total energy \(E\)—with each state being equally probable. For the fictitious system, we have

\[\hat{W} (\hat{E}) \, = \, \dfrac{E_0}{N! \, h^{3N+1}} \, \int \mathrm{d} \hat{p}^{}_{\! s} \int \mathrm{d} \hat{s} \, \prod_{\alpha = 1}^N \bigg( \int \! \mathrm{d} \boldsymbol{\hat{p}}^\alpha \int \! \mathrm{d} \boldsymbol{\hat{q}}^\alpha \bigg) \, \bigg\{ \delta \Big( \hat{\mathcal{H}} \big( \{ \boldsymbol{\hat{q}}^\alpha \}, \, \{ \boldsymbol{\hat{p}}^\alpha \}, \, \hat{s}, \, \hat{p}^{}_{\! s} \big) \, - \, \hat{E} \mkern1mu \Big) \bigg\} ~,\]

where \(\delta (\ldots)\) is the Dirac \(\delta\)-function, \(h\) is Planck’s constant, and \(E_0\) is an arbitrary constant with units of energy. Note the factor of \(E_0 / h^{3N + 1}\) ensures \(\hat{W}\) is dimensionless, and the factor of \(N!\) accounts for the \(N\) particles being indistinguishable. At this point, Nosé recognized it was useful to recast the partition function in terms of the real positions and momenta of the particles, rather than their fictitious counterparts. We already noted \(\boldsymbol{\hat{q}}^\alpha = \boldsymbol{q}^\alpha\), and now recognize

\[\boldsymbol{\hat{p}}^\alpha \, = \, m \hat{s}^2 \, \dfrac{\mathrm{d} \boldsymbol{\hat{q}}^\alpha}{\mathrm{d} \hat{t}} \, = \, m \hat{s} \, \dfrac{\mathrm{d} \boldsymbol{q}^\alpha}{\mathrm{d} t} ~, \qquad \text{for which} \qquad \boldsymbol{p}^\alpha \, = \, \dfrac{\boldsymbol{\hat{p}}^\alpha}{\hat{s} \,} \, = \, m \mkern1mu \boldsymbol{\dot{q}}^\alpha\]

is the physical momentum. The fictitious Hamiltonian is then expressed in terms of the physical microstate \(\boldsymbol{\Gamma}\) as

\[\hat{\mathcal{H}} \big( \boldsymbol{\Gamma}, \, \hat{s}, \, \hat{p}^{}_{\! s} \big) \, = \, \mathcal{H} ( \boldsymbol{\Gamma} ) \, + \, \dfrac{\hat{p}^2_{\! s}}{2 Q} \, + \, g \, k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \, \ln \hat{s} ~,\]

where

\[\mathcal{H} ( \boldsymbol{\Gamma} ) \, = \, \sum_{\alpha = 1}^N \dfrac{\boldsymbol{p}^\alpha \boldsymbol{\cdot} \boldsymbol{p}^\alpha}{2 m} \, + \, U \big( \{\boldsymbol{q}^\alpha\} \big)\]

is the physical Hamiltonian. Moreover, we recognize \(\mathrm{d} \boldsymbol{\hat{p}}^\alpha = \hat{s}^3 \mathrm{d} \boldsymbol{p}^\alpha\) in a three-dimensional system, and express the partition function as

\[\hat{W} ( \hat{E} ) \, = \, \dfrac{E_0}{N! \, h^{3N+1}} \int \! \mathrm{d} \hat{p}^{}_{\! s} \int \! \mathrm{d} \hat{s} \int \! \mathrm{d} \boldsymbol{\Gamma} \, \hat{s}^{3 N} \, \delta \bigg( \mathcal{H} ( \boldsymbol{\Gamma} ) \, + \, \dfrac{\hat{p}^2_{\! s}}{2 Q} \, + \, g \mkern1mu k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \mkern1mu \ln \hat{s} \, - \, \hat{E} \mkern1mu \bigg) ~.\]

Note that here we introduced the shorthand \(\mathrm{d} \boldsymbol{\Gamma} := \prod_{\alpha = 1}^N ( \int \mathrm{d} \boldsymbol{q}^\alpha \int \mathrm{d} \boldsymbol{p}^\alpha )\). Importantly, the argument of the \(\delta\)-function now contains only a single term involving \(\hat{s}\). By applying the \(\delta\)-function identity

\[\delta \big( \hat{f} (\hat{s}) \big) \, = \, \dfrac{\delta ( \hat{s} - \hat{s}_0 )}{ \lvert \hat{f}{\mkern1mu}'(\hat{s}_0) \rvert }\]

for any function \(\hat{f} (\hat{s})\) with \(f (\hat{s}_0) = 0\), and introducing \(\beta_{\mathrm{ext}} := ( k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} )^{-1}\) as the inverse temperature, we find

\[\delta \bigg( \mathcal{H} ( \boldsymbol{\Gamma} ) \, + \, \dfrac{ \hat{p}^2_{\! s}}{2 Q} \, + \, g \, k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \, \ln \hat{s} \, - \, \hat{E} \mkern1mu \bigg) \, = \, \dfrac{\hat{s}_0 \mkern1mu \beta_{\mathrm{ext}}}{g} \, \delta \big( \hat{s} \, - \, \hat{s}_0 \big) ~,\]

where

\[\hat{s}_0 \, = \, \exp \bigg\{ \dfrac{\beta_{\mathrm{ext}}}{g} \, \bigg( \hat{E} \, - \, \hat{\mathcal{H}} ( \boldsymbol{\Gamma} ) \, - \, \dfrac{\hat{p}^2_{\! s}}{2 Q} \bigg) \bigg\} ~.\]

The integration over \(\hat{s}\) is now easily carried out. After some rearrangement, we obtain

\[\hat{W} (\hat{E}) \, = \, \dfrac{ \beta^{}_{\mathrm{ext}} \mkern1mu E_0 } {g \, N! \, h^{3N+1}} \, \int \! \mathrm{d} \hat{p}^{}_{\! s} \int \! \mathrm{d} \boldsymbol{\Gamma} \, \exp \Bigg\{ \dfrac{ \beta^{}_{\mathrm{ext}} (3N + 1)}{g} \bigg( \hat{E} \, - \, \mathcal{H} ( \boldsymbol{\Gamma} ) \, - \, \dfrac{\hat{p}^{2}_{\! s}}{2 Q} \bigg) \Bigg\} ~.\]

As the real system is meant to sample the canonical distribution, we choose

\[g \, = \, 3 \mkern1mu N \, + \, 1\]

to recover the Boltzmann factor. The partition function of the fictitious system can then be expressed as

\[\hat{W} (\hat{E}) \, = \, \dfrac{ \beta^{}_{\mathrm{ext}} \mkern1mu E_0 \, \exp ( \beta^{}_{\mathrm{ext}} \mkern1mu \hat{E} ) } {(3N + 1) \, N! \, h^{3N+1}} \, \int \! \mathrm{d} \hat{p}^{}_{\! s} \exp \bigg\{ - \dfrac{\beta^{}_{\mathrm{ext}} \, \hat{p}_{\! s}^2}{2 Q} \bigg\} \, \int \! \mathrm{d} \boldsymbol{\Gamma} \, \exp \Big\{ - \beta^{}_{\mathrm{ext}} \mkern1mu \mathcal{H} ( \boldsymbol{\Gamma} ) \Big\} ~.\]

The Gaussian integral over \(\hat{p}^{}_{\! s}\) can be calculated exactly, with which the fictitious microcanonical partition function is expressed as

\[\hat{W} (\hat{E}) \, = \, \dfrac{ E_0 \, \sqrt{ 2 \pi Q \beta^{}_{\mathrm{ext}} \,} \, \exp ( \beta^{}_{\mathrm{ext}} \mkern1mu \hat{E} ) }{(3N + 1) \, h} \, Q^{}_{ N V T_{\mathrm{ext}}} ~.\]

Here,

\[Q^{}_{ N V T_{\mathrm{ext}}} \, := \, \dfrac{1}{N! \, h^{3N}} \, \int \mathrm{d} \boldsymbol{\Gamma} \, \exp \Big\{ - \beta^{}_{\mathrm{ext}} \mkern1mu \mathcal{H} ( \boldsymbol{\Gamma} ) \Big\}\]

is the canonical partition function of a system of \(N\) particles in a volume \(V\), maintained at a constant temperature \(T_{\mathrm{ext}}\).

Significance

The two previous equations contain a crucial result from the analysis of Nosé, and it is useful to further highlight their importance. Let us imagine simulating the fictitious system at a constant energy \(\hat{E}\): each microstate with \(\hat{\mathcal{H}} (\boldsymbol{\hat{q}}^\alpha, \boldsymbol{\hat{p}}^\alpha, \hat{s}, \hat{p}^{}_{\! s}) \, = \, \hat{E}\) is equally likely. However, if we then only examine the physical degrees of freedom, we find the probability \(\mathcal{P}\) of observing a physical microstate \(\boldsymbol{\Gamma}\) satisfies

\[\mathcal{P} (\boldsymbol{\Gamma}) \, \sim \, \exp \Big\{ - \beta_{\mathrm{ext}} \mkern1mu \mathcal{H} (\boldsymbol{\Gamma}) \Big\} ~,\]

which is exactly the probability distribution of a system at the desired temperature \(T_{\mathrm{ext}}\). We have thus shown that when the fictitious system is in the microcanonical ensemble, the real system samples the canonical ensemble at the desired temperature \(T_{\mathrm{ext}}\).

Practical considerations

The above result concludes our presentation of Nose’s theoretical developments. In practice, standard MD methods to simulate the fictitious system in the microcanonical ensemble would generate

\[\boldsymbol{\hat{q}}^\alpha (\hat{t}), \, \qquad \boldsymbol{\hat{p}}^\alpha (\hat{t}), \, \qquad \hat{s} (\hat{t}), \, \qquad \text{and} \qquad\quad \hat{p}^{}_{\! s} (\hat{t})\]

at a set of discrete fictitious times \(\hat{t} \in \{ \hat{t}_0, \, \hat{t}_1, \, \hat{t}_2, \, \ldots, \, \hat{t}_T \}\). Generally, the time interval \(\Delta \hat{t}\) in simulations is constant, for which \(\hat{t}_j = j \Delta \hat{t}\). From the scale factor \(\hat{s}\) at discrete \(\hat{t}\), we approximate the real time according to

\[t_i \, = \, \int_0^{t_i} \mathrm{d} t \, = \, \int_0^{\hat{t}_i} \dfrac{\mathrm{d} \hat{t}}{\hat{s} (\hat{t})} \, \approx \, \sum_{j = 0}^i \dfrac{w_j \, \Delta \hat{t}}{\hat{s} (\hat{t}_j)} ~,\]

where the \(\{ w_j \}\) are the weights of some quadrature formula—for example, the trapezoidal rule. The real positions and momenta are then given by

\[\boldsymbol{q}^\alpha (t_i) \, = \, \boldsymbol{\hat{q}}^\alpha (\hat{t}_i) \qquad \text{and} \qquad \boldsymbol{p}^\alpha (t_i) \, = \, \dfrac{\boldsymbol{\hat{p}}^\alpha (\hat{t}_i)}{ \hat{s} (\hat{t}_i) } ~,\]

with which all other observables can be calculated. Importantly, the calculated \(\boldsymbol{\Gamma} (t_i)\) appropriately sample the canonical ensemble at the specified temperature \(T_{\mathrm{ext}}\).

The extension by Hoover

Nosé’s equations mark a significant advancement in constant temperature MD simulations. However, the equations are slightly inconvenient to use in practice because of the need to calculate \(t\) from \(\hat{t}\), and the resultant non-uniform (real) time intervals \(\Delta t\). In the year following Nosé’s seminal work, Hoover recognized it was possible to simulate Nosé’s equations in terms of the real (rather than fictitious) time—and do away with any time scaling altogether.³ In what follows, we present the main results of Hoover’s theoretical developments.

We begin with the equations of motion for the extended, fictional system. By substituting the fictional Hamiltonian into Hamilton’s equations of motion, we find

\[\dfrac{\mathrm{d} \boldsymbol{\hat{q}}^\alpha}{\mathrm{d} \hat{t}} \, = \, \dfrac{\boldsymbol{\hat{p}}^\alpha}{m \hat{s}^2} ~, \qquad\qquad \dfrac{\mathrm{d} \boldsymbol{\hat{p}}^\alpha}{\mathrm{d} \hat{t}} \, = \, - \, \dfrac{\partial U}{\partial \boldsymbol{\hat{q}}^\alpha} ~, \qquad\qquad \dfrac{\mathrm{d} \hat{s}}{\mathrm{d} \hat{t}} \, = \, \dfrac{\hat{p}^{}_{\! s}}{Q} ~,\]

and

\[\dfrac{\mathrm{d} \hat{p}^{}_{\! s}}{\mathrm{d} \hat{t}} \, = \, \sum_{\alpha = 1}^N \dfrac{\boldsymbol{\hat{p}}^\alpha \boldsymbol{\cdot} \boldsymbol{\hat{p}}^\alpha}{m \hat{s}^3} \, - \, \dfrac{g \, k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}}}{\hat{s}} ~.\]

Note that we once again express the dynamical equation for \(\hat{p}^{}_{\! s}\) in terms of \(g\), rather than the previously determined value of \(3N + 1\), due to a subtlety that will be discussed subsequently. At this point, Hoover recognized the utility in expressing these equations entirely in terms of real quantities. To this end, we introduce the following fictitious quantities that depend on the real time (note that we drop the ‘hat’ accent):

\[\eta \, := \, \ln \hat{s} \qquad\quad \text{and} \qquad\quad p^{}_{\mkern-3mu \eta} \, := \, \hat{p}^{}_{\! s} ~.\]

By recognizing

\[\dfrac{\mathrm{d} \hat{s}}{\mathrm{d} \hat{t}} \, = \, \dfrac{1}{\hat{s}} \, \dfrac{\mathrm{d} \hat{s}}{\mathrm{d} t} \, = \, \dfrac{\mathrm{d} \ln \hat{s}}{\mathrm{d} t} \, = \, \dfrac{\mathrm{d} \eta}{\mathrm{d} t} \, = \, \dot{\eta}\]

and defining

\[\boldsymbol{f}^\alpha \, := \, - \, \dfrac{\partial U}{\partial \boldsymbol{q}^\alpha}\]

as the force on the \(\alpha^{\mathrm{th}}\) particle, we express the governing equations in terms of the real time \(t\) as

\[\boldsymbol{\dot{q}}^\alpha \, = \, \dfrac{\boldsymbol{p}^\alpha}{m} ~, \qquad\quad \boldsymbol{\dot{p}}^\alpha \, = \, \boldsymbol{f}^\alpha \, - \, \bigg( \dfrac{p^{}_{\mkern-3mu \eta}}{Q} \bigg) \, \boldsymbol{p}^\alpha ~, \qquad\quad \dot{\eta} \, = \, \dfrac{p^{}_{\mkern-3mu \eta}}{Q} ~,\]

and

\[\dot{p}^{}_{\mkern-3mu \eta} \, = \, \sum_{\alpha = 1}^N \dfrac{\boldsymbol{p}^\alpha \boldsymbol{\cdot} \boldsymbol{p}^\alpha}{m} \, - \, g \, k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} ~.\]

In this manner, the dynamics of the fictitious system are simulated as a function of the real time \(t\), rather than the fictitious time \(\hat{t}\). Consequently, the conversion between \(\hat{t}\) and \(t\) is no longer required, and an MD simulation generates

\[\boldsymbol{q}^\alpha (t) ~, \qquad \boldsymbol{q}^\alpha (t) ~, \qquad \eta (t) ~, \qquad \text{and} \qquad p^{}_{\mkern-3mu \eta} (t)\]

at a set of discrete and uniformly spaced real times \(t_j = j \mkern1mu \Delta t\). The practical consequences of this formulation are significant, and the above equations—referred to as the Nosé–Hoover equations—are generally used in numerical implementations.

Ergodicity, time averages, and the value of \(g\)

Recall that in MD simulations, the expectation value of a generic phase variable \(B (\boldsymbol{\Gamma})\) is calculated as a time average. This value is assumed to be equal to the ensemble average with the ergodic hypothesis. Thus, for a fictitious system evolved in the scaled time \(\hat{t}\), we expect

\[\Big\langle B \big( \{ \boldsymbol{q}^\alpha \}, \{ \boldsymbol{p}^\alpha \} \big) \Big\rangle_{\! \mathrm{c}} \, = \, \Big\langle B \big( \{ \boldsymbol{\hat{q}}^\alpha \}, \{ \boldsymbol{\hat{p}}^\alpha / \hat{s} \} \big) \Big\rangle_{\! \hat{t}} \, := \, \lim_{\hat{\tau} \rightarrow \infty} \dfrac{1}{\hat{\tau}} \int_0^{\hat{\tau}} B \big( \{ \boldsymbol{\hat{q}}^\alpha (\hat{t}) \}, \{ \boldsymbol{\hat{p}}^\alpha (\hat{t}) / \hat{s} (\hat{t}) \} \big) ~ \mathrm{d} \hat{t} ~,\]

where the left-hand side denotes the expectation value in the canonical ensemble, i.e. an ensemble average. A complication arises when transitioning from the Nosé to the Nosé–Hoover formalism, however, because the times \(\hat{t}\) and \(t\) are related by a non-constant scale factor. Thus, some care must be taken to ensure that a uniform sampling in real time with the Nosé–Hoover dynamics correctly determines equilibrium averages. To this end, we begin by expressing the average value of a phase variable in real time as

\[\Big\langle B \big( \{ \boldsymbol{q}^\alpha \}, \{ \boldsymbol{p}^\alpha \} \big) \Big\rangle_{\! t} \, := \, \lim_{\tau \rightarrow \infty} \dfrac{1}{\tau} \int_0^{\tau} B \big( \{ \boldsymbol{q}^\alpha (t) \}, \{ \boldsymbol{p}^\alpha (t) \} \big) ~ \mathrm{d} t ~.\]

Next, we multiply and divide by \(\hat{\tau}\) and rearrange terms to obtain

\[\Big\langle B \big( \{ \boldsymbol{q}^\alpha \}, \{ \boldsymbol{p}^\alpha \} \big) \Big\rangle_{\! t} \, = \, \lim_{\tau \rightarrow \infty} \Bigg[ \Bigg( \dfrac{1}{\hat{\tau}} \int_0^{\tau} B \big( \{ \boldsymbol{q}^\alpha (t) \}, \{ \boldsymbol{p}^\alpha (t) \} \big) ~ \mathrm{d} t \Bigg) \bigg/ \bigg( \dfrac{\tau}{\hat{\tau}} \bigg) \Bigg] ~.\]

At this point, we change variables from \(t\) to \(\hat{t}\), for which \(\mathrm{d} t = \mathrm{d} \hat{t} / \hat{s}\) and \(\tau = \int_0^{\hat{\tau}} \mathrm{d}\hat{t} / \hat{s}\). After some algebra, the time average of a phase variable in the Nosé–Hoover dynamics can be expressed in terms of canonical ensemble averages as

\[\Big\langle B \big( \{ \boldsymbol{q}^\alpha \}, \{ \boldsymbol{p}^\alpha \} \big) \Big\rangle_{\! t} \, = \, \dfrac{ \Big\langle B \big( \{ \boldsymbol{q}^\alpha \}, \{ \boldsymbol{p}^\alpha \} \big) / \hat{s} \Big\rangle_{\! \mathrm{c}} }{ \big\langle 1 / \hat{s} \big\rangle_{\! \mathrm{c}} } ~.\]

The canonical averages are now expressed in terms of an integral over the extended phase space. Since the quantity to be averaged depends on \(\hat{s}\), we use a form of the fictitious partition function before the integral over \(\hat{s}\) is taken. We find

\[\big\langle B ( \boldsymbol{\Gamma} ) \big\rangle_{\! t} \, = \, \dfrac{\displaystyle \int \! \mathrm{d} \hat{p}^{}_{\! s} \int \! \mathrm{d} \hat{s} \int \! \mathrm{d} \boldsymbol{\Gamma} \ B ( \boldsymbol{\Gamma} ) \, \hat{s}^{3 N - 1} \, \delta \bigg( \mathcal{H} ( \boldsymbol{\Gamma} ) \, + \, \dfrac{\hat{p}^2_{\! s}}{2 Q} \, + \, g \mkern1mu k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \mkern1mu \ln \hat{s} \, - \, \hat{E} \mkern1mu \bigg) }{\displaystyle \int \! \mathrm{d} \hat{p}^{}_{\! s} \int \! \mathrm{d} \hat{s} \int \! \mathrm{d} \boldsymbol{\Gamma} \ \hat{s}^{3 N - 1} \, \delta \bigg( \mathcal{H} ( \boldsymbol{\Gamma} ) \, + \, \dfrac{\hat{p}^2_{\! s}}{2 Q} \, + \, g \mkern1mu k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \mkern1mu \ln \hat{s} \, - \, \hat{E} \mkern1mu \bigg) }\]

Importantly, in both integrands \(\hat{s}\) is now raised to the power \(3N - 1\), rather than \(3N\). Applying the \(\delta\)-function identity to both the numerator and denominator and simplifying the result yields

\[\big\langle B ( \boldsymbol{\Gamma} ) \big\rangle_{\! t} \, = \, \dfrac{\displaystyle \int \! \mathrm{d} \boldsymbol{\Gamma} \ B ( \boldsymbol{\Gamma} ) \, \exp \bigg\{ - \beta_{\mathrm{ext}} \mathcal{H} ( \boldsymbol{\Gamma} ) \cdot \dfrac{3N}{g} \mkern1mu \bigg\} }{\displaystyle \int \! \mathrm{d} \boldsymbol{\Gamma} \ \exp \bigg\{ - \beta_{\mathrm{ext}} \mathcal{H} ( \boldsymbol{\Gamma} ) \cdot \dfrac{3N}{g} \mkern1mu \bigg\} } ~.\]

Thus, to recover the canonical ensemble averages within the Nosé–Hoover formalism, we set

\[g \, = \, 3 \mkern1mu N\]

rather than the previously determined value of \(3N + 1\), which is appropriate only within the Nosé formalism.²^,⁴

Additional insights

With the Nosé–Hoover equations expressed in terms of the real time \(t\), we can draw additional insights about the dynamics of the extended system. First, the quantity \((p^{}_{\mkern-3mu \eta} / Q)\) acts as a drag coefficient: each physical particle feels an additional force \(- (p^{}_{\mkern-3mu \eta} / Q) \boldsymbol{p}^\alpha\) due to the heat bath. In this case, however, the drag coefficient dynamically evolves in time and can be either positive or negative. Moreover, recall from the equipartition thoerem that

\[k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \, = \, \dfrac{1}{N_{\mathrm{DOFs}}} \, \bigg\langle \sum_{\alpha = 1}^N \dfrac{\boldsymbol{p}^\alpha \boldsymbol{\cdot} \boldsymbol{p}^\alpha}{m} \bigg\rangle ~,\]

where \(N_{\mathrm{DOFs}}\) is the number of physical degrees of freedom. Consequently, we define an “instantaneous temperature” \(T(t)\) according to

\[k_{\mathrm{B}} \mkern1mu T(t) \, := \, \dfrac{1}{N_{\mathrm{DOFs}}} \, \sum_{\alpha = 1}^N \dfrac{\boldsymbol{p}^\alpha \boldsymbol{\cdot} \boldsymbol{p}^\alpha}{m} ~.\]

For a system of \(N\) particles with no global constraints on the dynamics, \(N_{\mathrm{DOFs}} = 3N\), for which the dynamical equation governing \(p^{}_{\mkern-3mu \eta}\) can be written as

\[\dot{p}^{}_{\mkern-3mu \eta} \, = \, N_{\mathrm{DOFs}} \Big( k_{\mathrm{B}} \mkern1mu T(t) \, - \, k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \Big) ~.\]

We thus observe a simple physical interpretation for the drag force, which affects all particle momenta:

it “cools down” a system that is too “hot,” for which \(T(t) > T_{\mathrm{ext}}\)
it “heats up” a system that is too “cold,” for which \(T(t) < T_{\mathrm{ext}}\)

We note that in situations with periodic boundary conditions in all three directions, the total linear momentum of the system is conserved. In this case, \(N_{\mathrm{DOFs}} = 3 N - 3 = g\). The latter equality arises because the integration over phase space involves \(N - 1\) vector momenta, rather than \(N\). Thus, the above equation for \(\dot{p}^{}_{\mkern-3mu \eta}\) remains valid.

We close by noting that the Nosé equations can be derived from the fictitious Hamiltonian \(\mathcal{H} (\boldsymbol{\hat{q}}^\alpha, \boldsymbol{\hat{p}}^\alpha, \hat{s}, \hat{p}^s)\). However, the same is not true for the Nosé–Hoover equations, in which the fundamental variables are the real (rather than fictitious) particle positions and momenta. Nonetheless, the value of the fictitious Hamiltonian remains a constant of motion in the extended system, and can be expressed in terms of real quantities. We thus define the quantity

\[E_{\mathrm{NH}} \, := \, \sum_{\alpha = 1}^N \dfrac{\boldsymbol{p}^\alpha \boldsymbol{\cdot} \boldsymbol{p}^\alpha}{2 m} \, + \, U \big( \{\boldsymbol{q}^\alpha\} \big) \, + \, \dfrac{p^2_{\mkern-3mu \eta}}{2 Q} \, + \, g \, k_{\mathrm{B}} \mkern1mu T_{\mathrm{ext}} \, \eta ~,\]

which has the dimensions of an energy and is conserved by the Nosé–Hoover equations.

Numerical details

Since the Nosé–Hoover equations are non-Hamiltonian, they cannot be solved numerically using the usual integration methods from constant energy MD simulations. The numerical simulation of the Nosé–Hoover equations will be discussed in a future post. In the meantime, we recommend the excellent textbook by Allen and Tildesley.⁵

References

S. Nosé. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 52, 255–268 (1984) ↩ ↩²
S. Nosé. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 81, 511–519 (1984) ↩ ↩² ↩³
W.G. Hoover. Canonical Dynamics: Equilibrium phase-space distributions. Phys. Rev. A 31, 1695–1697 (1985) ↩ ↩²
D.J. Evans and B.L. Holian. The Nosé–Hoover thermostat. J. Chem. Phys. 83, 4069–4074 (1985) ↩
M.P. Allen and D.J. Tildesley. Computer Simulation of Liquids, Chapter 3.8.2, 2^nd Edition. New York: Oxford University Press, 2017. See also the Companion Website and Github Repository. ↩