next up previous
Next: Calculation of the shear Up: POINT DEFECTS, LATTICE STRUCTURE Previous: Predictor - corrector method

Nose-Hoover algorithm

The simplest extension of the NVE ensemble is the canonical one (NVT), where the number of particles, the volume and the temperature are fixed to the prescribed values. The temperature $T$, in contrast to the number of particles $N$ and volume $V$, is an intensive parameter. The extensive counterpart is kinetic energy, related to $T$ through:
\begin{displaymath}
T=\frac{2}{3}\frac{E_{kin}}{Nk_b}=\frac{1}{3Nk_b}\sum_{i=1}^{N}\frac{p_i^2}{m}
\end{displaymath} (9.1)

There are different methods which control the temperature: the differential, proportional, stochastic and the integral thermostat [59]. The integral thermostat method was selected for our MD simulations. The integral thermostat method (which sometimes called the extended system method or the Nose-Hoover algorithm [60]) introduces additional degrees of freedom into the system's Hamiltonian, for which equation of motion can be derived. These equations for the additional degrees of freedom are integrated together with "usual" equations for spatial coordinates and momenta.

The idea of the method proposed by Nose [59,96] was to reduce the effect of an external system, acting as heat reservoir, to an additional degree of freedom. This heat reservoir controls the temperature of the given system, i.e. the temperature fluctuates around target value. Actually, the thermal interaction between the heat reservoir and the system results in exchange of the kinetic energy between them. Nose introduced two sets of variables: real $\{p_i,q_i\}$ and virtual $\{\pi_i,\rho_i\}$ ones. The virtual variables are derived from the so called Sundman's transformation [62]:

\begin{displaymath}
s=\frac {d \tau}{dt}
\end{displaymath} (9.2)

where $\tau$ is the virtual time, $t$ is the real time and $s$ is a resulting scaling factor, which also treated as a dynamical variable. The transformation from the virtual variables $\{\pi_i,\rho_i\}$ to the real ones $\{p_i,q_i\}$ is performed according to:
\begin{displaymath}
p_i=\pi_i
\end{displaymath} (9.3)


\begin{displaymath}
q_i=\rho_i
\end{displaymath} (9.4)

The introduction of the effective mass $ M_s$ connects also a momentum to the additional degree of freedom $\pi_s$. The resulting Hamiltonian, expressed in terms of the virtual coordinates can be written as:
\begin{displaymath}
H^*=\sum_{i=1}^{N}\frac{\pi_i^2}{2ms}+U(\rho_1,\rho_2,.....,\rho_N)+\frac{\pi_s^2}{2M_s}+gk_bTln(s)
\end{displaymath} (9.5)

where $g=3N+1$ is the number degrees of freedom of the extended system ($N$ particles + 1 the new degree of freedom). It was shown, that this Hamiltonian $H^*$ leads to a density of probability in phase space, corresponding to the canonical ensemble[61].

The equations of motion obtained from the Hamiltonian are:

\begin{displaymath}
\frac{d\vec \rho_i}{d\tau}=\frac{\partial H^*}{\partial \vec \pi _i}=\frac{\vec \pi_i}{ms^2}
\end{displaymath} (9.6)


\begin{displaymath}
\frac{d\vec \pi_i}{d\tau}=-\frac{\partial H^*}{\partial \vec \rho _i}=-\frac{\partial U}{\partial \vec \rho_i}
\end{displaymath} (9.7)


\begin{displaymath}
\frac{ds}{d\tau}=\frac{\partial H^*}{\partial \pi _s}=\frac{\pi_s}{M_s}
\end{displaymath} (9.8)


\begin{displaymath}
\frac{d\pi_s}{d\tau}=-\frac{\partial H^*}{\partial s}=\frac{gk_bT}{s}+\sum_{i=1}^{N}\frac{\pi_i^2}{2ms}
\end{displaymath} (9.9)

If one transforms these equations back into the real variables $\{p_i,q_i\}$ and introduces a new variable $\zeta$:

\begin{displaymath}
\zeta=s\frac{ds}{dt}=s\frac{ds}{d\tau}\frac{d\tau}{dt}
=s\frac{\partial H^*}{d\pi_s}\frac{d\tau}{dt}=s^2\frac{\pi_s}{M_s}
\end{displaymath} (9.10)

then one obtains (according to Hoover [63]):
\begin{displaymath}
\frac{d\vec q_i}{d t}=\frac{\vec p_i}{m_i}
\end{displaymath} (9.11)


\begin{displaymath}
\frac{d\vec p_i}{d\tau}=-\frac{\partial U}{\partial \vec q _i}-\zeta \vec p_i
\end{displaymath} (9.12)


\begin{displaymath}
\frac{\partial ln(s)}{\partial t}=\zeta
\end{displaymath} (9.13)


\begin{displaymath}
\frac{d \zeta}{dt}=\frac{1}{M_s}\left ( \sum_{i=1}^{N}\frac{ p_i^2 }{2m_i}-gk_bT\right),~~~~ p_i\equiv \vert\vec p_i\vert
\end{displaymath} (9.14)

These equations describe the Noose-Hoover thermostat [64]. The parameter $ M_s$ is a thermal inertia parameter, which determines rate of the heat transfer. The value of this parameter must be set carefully, because if it is chosen to be too small the phase space of the system will not be canonical [65], and it is chosen to be too large the temperature control will not be efficient. In our simulation, for example, we set $M_s=1/6$. By means of the Noose-Hoover thermostat we can impose time averaged value of temperature to be equal to the prescribed value.
next up previous
Next: Calculation of the shear Up: POINT DEFECTS, LATTICE STRUCTURE Previous: Predictor - corrector method
2003-01-15