Nosé- Hoover Thermostat

In the Nosé-Hoover algorithm [18] Newton's equations of motion are modified to read:

$\displaystyle {d \mbox{$\underline{r}$}(t) \over d t}$	$\textstyle =$	$\displaystyle \mbox{$\underline{v}$}(t)$
$\displaystyle {d \mbox{$\underline{v}$}(t) \over d t}$	$\textstyle =$	$\displaystyle {\mbox{$\underline{f}$}(t)\over m} - \chi(t) \mbox{$\underline{v}$}(t)$	(2.181)

The friction coefficient, $\chi$ , is controlled by the first order differential equation

where $\tau_T$ is a specified time constant (normally in the range [0.5, 2] ps). In DL_POLY_2
$\chi$ is stored at half timesteps as it has dimensions of (1/time). The integration takes place as:

$\displaystyle \chi(t+{1\over2}\Delta t)$	$\textstyle \leftarrow$	$\displaystyle \chi(t-{1\over 2}\Delta t) + {\Delta t \over \tau_T^2}\left( {{\cal T}\over T_{\rm ext}} - 1\right)$
$\displaystyle \chi(t)$	$\textstyle \leftarrow$	$\displaystyle {1\over 2} \left[\chi(t-{1\over 2}\Delta t)+\chi(t+{1\over 2}\Delta t)\right]$
$\displaystyle \mbox{$\underline{v}$}(t + {1\over 2} \Delta t)$	$\textstyle \leftarrow$	$\displaystyle \mbox{$\underline{v}$}(t-{1\over 2} \Delta t) + \Delta t\left[ {\mbox{$\underline{f}$}(t)\over m} - \chi(t)\mbox{$\underline{v}$}(t)\right]$
$\displaystyle \mbox{$\underline{v}$}(t)$	$\textstyle \leftarrow$	$\displaystyle {1\over 2} \left[\mbox{$\underline{v}$}(t-{1\over 2}\Delta t) +\mbox{$\underline{v}$}(t+{1\over 2}\Delta t)\right]$
$\displaystyle \mbox{$\underline{r}$}(t+\Delta t)$	$\textstyle \leftarrow$	$\displaystyle \mbox{$\underline{r}$}(t) + \Delta t \; \mbox{$\underline{v}$}(t+{1\over 2}\Delta t)$	(2.183)

Since $\mbox{$\underline{v}$}(t)$ is required to calculate ${\cal T}$ and itself, the algorithm requires several iterations to obtain self consistency. In DL_POLY_2 the number of iterations is set to 3 (4 if the system has bond constraints). The iteration procedure is started with the standard Verlet leapfrog prediction of $\mbox{$\underline{v}$}(t)$ and ${\cal T}$ . The conserved quantity is derived from the extended Hamiltonian for the system which, to within a constant, is the Helmholtz free energy:

If bond constraints are present an extra iteration is required due to the call to the SHAKE routine. Note that the SHAKE corrections need only be applied during the first iteration as subsequent to this the velocities, relative to the bond c.o.m. velocity, will be orthogonal to the bond vectors. The velocity scaling imposed by the thermostat is isotropic so does not destroy this orthogonality. The algorithm is implemented in the DL_POLY routines NVT_H0 and NVT_H1, the latter being for systems with bond constraints.