Thermostats

Standard molecular dynamics generates the microcanonical or $NVE$ ensemble where in addition the total momentum is conserved [13]. The temperature is not a control variable and cannot be preselected and fixed. But it is evident that also within molecular dynamics the possibility to control the average temperature (as obtained from the average kinetic energy of the nuclei and the energy equipartition theorem) is welcome for physical reasons. A deterministic algorithm of achieving temperature control in the spirit of extended system dynamics [30] by a sort of dynamical friction mechanism was devised by Nosé and Hoover [26,27], see e.g. Refs. [12,26,13] for reviews of this technique. Thereby, the canonical or $NVT$ ensemble is generated in the case of ergodic dynamics.

It is well-known that the standard Nosé-Hoover thermostat method suffers from non-ergodicity problems for certain classes of Hamiltonians, such as the harmonic oscillator [27]. A closely related technique, the so-called Nosé-Hoover-chain thermostat [32], cures that problem and assures ergodic sampling of phase space even for the pathological harmonic oscillator. This is achieved by thermostatting the original thermostat by another thermostat, which in turn is thermostatted and so on. In addition to restoring ergodicity even with only a few thermostats in the chain, this technique is found to be much more efficient in imposing the desired temperature. The underlying equations of motion read

$$M_I {\bf\ddot R}_I = - {\partial U ({\bf R}^N) \over \partial {\bf R}_I } - M_I {\dot \xi }_1 {\bf\dot R}_I$$ (14)

$$Q_1^{\rm n} {\ddot \xi}_1 = \left[ \sum_I M_I {\bf\dot R}_I^2 - g k_{\rm B} T \right] - Q_1^{\rm n} {\dot \xi}_1 {\dot \xi}_2$$

$$Q_k^{\rm n} {\ddot \xi}_k = \left[ Q_{k-1}^{\rm n} {\dot \xi}_{k-1}^2 - k_{\rm B} T \right] -... ...n} {\dot \xi}_k {\dot \xi}_{k+1} \left( 1-\delta_{kK} \right) \enspace \enspace k=2, \dots , K \enspace .$$

By inspection of Eq. (14) it becomes intuitively clear how the thermostat works: ${\dot \xi }_1$ can be considered as a dynamical friction coefficient. The resulting ``dissipative dynamics'' leads to non-Hamiltonian flow, but the friction term can acquire positive or negative sign according to its equation of motion. This leads to damping or acceleration of the nuclei and thus to cooling or heating if the instantaneous kinetic energy of the nuclei is higher or lower than $k_{\rm B}T$ which is preset. As a result, this extended system dynamics can be shown to produce a canonical ensemble in the subspace of the nuclear coordinates and momenta. In spite of being non-Hamiltonian, Nosé-Hoover (-chain) dynamics is also distinguished by conserving an energy quantity of the extended system; see Eq. (16).

The desired average physical temperature is given by $T$ and $g$ denotes the number of dynamical degrees of freedom to which the nuclear thermostat chain is coupled (i.e. constraints imposed on the nuclei have to be subtracted). It is found that this choice requires a very accurate integration of the resulting equations of motion (for instance by using a high-order Suzuki-Yoshida integrator [33]). The integration of these equations of motion is discussed in detail in Ref. [33] using the velocity Verlet algorithm. One of the advantages of the velocity Verlet integrator is that it can be easily used together with higher order schemes for the thermostats. Multiple time step techniques can be used and time reversibility of the overall algorithm is preserved.

       Integrate Thermostats for dt/2
       V(:) := V(:) + dt/(2*M(:))*F(:)
       R(:) := R(:) + dt*V(:)
       Calculate new forces F(:)
       V(:) := V(:) + dt/(2*M(:))*F(:)
       Integrate Thermostats for dt/2

The choice of the ``mass parameters'' assigned to the thermostat degrees of freedom should be made such that the overlap of their power spectra and the ones of the thermostatted subsystems is maximal [33]. The relations

$$Q_1^{\rm n} = {{gk_{\rm B}T} \over {\omega_{\rm n}^2}} \enspace , \hspace*{1.0cm} Q_k^{\rm n} = {{k_{\rm B}T} \over {\omega_{\rm n}^2}} \enspace ,$$ (15)

assures this if $\omega_{\rm n}$ is a typical phonon or vibrational frequency of the nuclear subsystem (say of the order of 2000 to 4000 cm$^{-1}$ ). There is a conserved energy quantity in the case of thermostatted molecular dynamics. This constant of motion reads

$$E_{\rm cons}^{\rm NVT} = \sum_I {1 \over 2} M_i {\bf\dot R}_i^2 + U({\bf R}^N ) + \sum_{k=... ...{\rm n}_k {\dot \xi}_k^2 + \sum_{k=2}^K k_{\rm B} T \xi_k + g k_{\rm B} T \xi_1$$ (16)

for Nosé-Hoover-chain thermostatted molecular dynamics.