Next: Constant-pressure integration with isotropic Up: Mechanics Previous: Extended Systems for other Contents

Nosé-Hoover Thermostatting

The Nosé-Hoover equations of motion are given by[45]

$\displaystyle \dot{{\bf r}}$	$\displaystyle =$	$\displaystyle {\bf v},$	(4.17)
$\displaystyle \dot{{\bf v}}$	$\displaystyle =$	$\displaystyle \frac{{\bf f}}{m} - \chi {\bf v} ,$	(4.18)
$\displaystyle \dot{\mathsf{A}}$	$\displaystyle =$	$\displaystyle \mathsf{A} \cdot$ skew $\displaystyle \left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right),$	(4.19)
$\displaystyle \dot{{\bf j}}$	$\displaystyle =$	$\displaystyle {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j} \right) -$ rot $\displaystyle \left(\mathsf{A}^{T} \cdot \frac{\partial V}{\partial \mathsf{A}} \right) - \chi {\bf j}.$	(4.20)

$\chi$ is an ``extra'' variable included in the extended system, and it is propagated using the first order equation of motion

$\displaystyle \dot{\chi} = \frac{1}{\tau_{T}^2} \left( \frac{T}{T_{\mathrm{target}}} - 1 \right).$

(4.21)

The instantaneous temperature is proportional to the total kinetic energy (both translational and orientational) and is given by

$\displaystyle T = \frac{2 K}{f k_B}$

(4.22)

Here,

is the total number of degrees of freedom in the system,

$\displaystyle f = 3 N + 2 N_{\mathrm{linear}} + 3 N_{\mathrm{non-linear}} - N_{\mathrm{constraints}},$

(4.23)

and

is the total kinetic energy,

$\displaystyle K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + ... ...1}{2} {\bf j}_i^T \cdot \overleftrightarrow{\mathsf{I}}_i^{-1} \cdot {\bf j}_i.$

(4.24)

$N_{\mathrm{linear}}$ is the number of linear rotors (i.e. with two non-zero moments of inertia), and $N_{\mathrm{non-linear}}$ is the number of non-linear rotors (i.e. with three non-zero moments of inertia).

In eq.(4.21), $\tau_T$ is the time constant for relaxation of the temperature to the target value. To set values for $\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use the tauThermostat and targetTemperature keywords in the meta-data file. The units for tauThermostat are fs, and the units for the targetTemperature are degrees K. The integration of the equations of motion is carried out in a velocity-Verlet style 2 part algorithm:

moveA:

$\displaystyle T(t)$	$\displaystyle \leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,$
$\displaystyle {\bf v}\left(t + h / 2\right)$	$\displaystyle \leftarrow {\bf v}(t) + \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) \chi(t)\right),$
$\displaystyle {\bf r}(t + h)$	$\displaystyle \leftarrow {\bf r}(t) + h {\bf v}\left(t + h / 2 \right),$
$\displaystyle {\bf j}\left(t + h / 2 \right)$	$\displaystyle \leftarrow {\bf j}(t) + \frac{h}{2} \left( {\bf\tau}^b(t) - {\bf j}(t) \chi(t) \right) ,$
$\displaystyle \mathsf{A}(t + h)$	$\displaystyle \leftarrow \mathrm{rotate} \left(h * {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) ,$
$\displaystyle \chi\left(t + h / 2 \right)$	$\displaystyle \leftarrow \chi(t) + \frac{h}{2 \tau_T^2} \left( \frac{T(t)} {T_{\mathrm{target}}} - 1 \right) .$

Here $\mathrm{rotate}(h * {\bf j} \overleftrightarrow{\mathsf{I}}^{-1})$ is the same symplectic Trotter factorization of the three rotation operations that was discussed in the section on the DLM integrator. Note that this operation modifies both the rotation matrix $\mathsf{A}$ and the angular momentum ${\bf j}$ . moveA propagates velocities by a half time step, and positional degrees of freedom by a full time step. The new positions (and orientations) are then used to calculate a new set of forces and torques in exactly the same way they are calculated in the doForces portion of the DLM integrator.

Once the forces and torques have been obtained at the new time step, the temperature, velocities, and the extended system variable can be advanced to the same time value.

moveB:

$\displaystyle T(t + h)$	$\displaystyle \leftarrow \left\{{\bf v}(t + h)\right\}, \left\{{\bf j}(t + h)\right\},$
$\displaystyle \chi\left(t + h \right)$	$\displaystyle \leftarrow \chi\left(t + h / 2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)} {T_{\mathrm{target}}} - 1 \right),$
$\displaystyle {\bf v}\left(t + h \right)$	$\displaystyle \leftarrow {\bf v}\left(t + h / 2 \right) + \frac{h}{2} \left( \frac{{\bf f}(t + h)}{m} - {\bf v}(t + h) \chi(t h)\right) ,$
$\displaystyle {\bf j}\left(t + h \right)$	$\displaystyle \leftarrow {\bf j}\left(t + h / 2 \right) + \frac{h}{2} \left( {\bf\tau}^b(t + h) - {\bf j}(t + h) \chi(t + h) \right) .$

Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to calculate as well as $\chi(t + h)$ , they indirectly depend on their own values at time . moveB is therefore done in an iterative fashion until $\chi(t + h)$ becomes self-consistent. The relative tolerance for the self-consistency check defaults to a value of 10 $^{-6}$ , but OOPSE will terminate the iteration after 4 loops even if the consistency check has not been satisfied.

The Nosé-Hoover algorithm is known to conserve a Hamiltonian for the extended system that is, to within a constant, identical to the Helmholtz free energy,[46]

$\displaystyle H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left( \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime \right).$

(4.25)

Poor choices of

or $\tau_T$ can result in non-conservation of $H_{\mathrm{NVT}}$ , so the conserved quantity is maintained in the last column of the .stat file to allow checks on the quality of the integration.

Bond constraints are applied at the end of both the moveA and moveB portions of the algorithm. Details on the constraint algorithms are given in section 4.7.1.

Next: Constant-pressure integration with isotropic Up: Mechanics Previous: Extended Systems for other Contents

Updated on January 16, 2006