Nosé Hoover Thermostat
From bio-physics-wiki
Molecular Dynamics (MD) Simulations can be used to calculate the position and momenta of a many particle system. The initial condition for particles of such simulations are often taken from a Maxwell distribution. However, since in this systems have constant energy, particle number and Volume, they represent Microcanonical Ensembles. Under experimental conditions it is easier to control the temperatur of a system, to keep it constant. In the so called Nosé Hoover Thermostat the hamiltonian does not only contain terms arising from the particles, but also an additional term. It can be shown that such a system represents a Canonical Ensemble$^{1}$ (T,N,V const.). The additional term simulates the energy transfer, thus keeping the termperature constant. This is why such equations are called thermostated or in this case the Nosé Hoover Thermostat. The Hamiltonian of the system reads \begin{align} H=\sum_i \frac{\mathbf{p}_i^2}{2ms^2}+\frac{1}{2} \sum_{ij,i \not = j} U(\mathbf{r}_i-\mathbf{r}_j) + \frac{p_s^2}{2Q}+gkT ln(s) \end{align} from the Hamiltionian we can derive the ODE system \begin{align} \frac{d\mathbf{r}_i}{dt}&=\frac{\partial H}{\partial \mathbf{p}_i}=\frac{\mathbf{p}_i}{ms^2}\\ \frac{ds}{dt}&=\frac{\partial H}{\partial p_s}=\frac{p_s}{Q}\\ \frac{d \mathbf{p}_i}{dt}&=\frac{\partial H}{\partial \mathbf{r}_i}=- \nabla _i U(R)\\ \frac{dp_s}{dt}&=\frac{\partial H}{\partial s}=\left( \frac{\sum_i \mathbf{p}_i^2}{ms^2}-gk_BT \right)/s \end{align}
According to Hoover$^{2}$ the equations for a one dimensional harmonic osciallator can be written as \begin{align} \dot{q}&=p\\ \dot{p}&=-q-\xi p\\ \dot{\xi}&=(p^2-1)/Q\\ \end{align}
Example
qpzx-Modul $\dot{\mathbf{x}}=\mathbf{G}(\mathbf{T})$ \begin{align} \dot{q}&=p\\ \dot{p}&=-q-z p\\ \dot{z}&=p^2 - T(q) - xz\\ \dot{T}&=1+ \varepsilon \arctan(q) \end{align} Fixed points: setting $\mathbf{G}(\mathbf{T})$ zero, $p$ must be zero, thus $q$ must be zero too, but $arctan(0)=\pm \infty$ diverges. ???
Linearising the system$^3$ about a fixpoint $\mathbf{T}_0$ gives \begin{align} \mathbf{G}(\mathbf{T}) \approx \underbrace{\mathbf{G}(\mathbf{T}_0)}_{=0} + \underbrace{\frac{\partial G_i(\mathbf{T}_0)}{\partial T_j}}_{\mathbf{J}} \cdot \delta T_j + o(\| \delta T_j \|) \end{align} The Jacobian of the linearised system reads \begin{align} \mathbf{J}=\begin{pmatrix} 0 & 1 & 0 &0 \\ -1 & -z & -p & 0 \\ 0 & 2p & -x & -1 \\ \frac{\varepsilon}{q^2+1} &0 &0 &0 \end{pmatrix} \end{align} At the fixed point $(q, p, z, T)=(xxxx)$ we get the linear System \begin{align} \delta \dot{T}_j=\mathbf{J} \cdot \delta T_j =\begin{pmatrix} \dot q\\ \dot p \\ \dot z\\ \dot T\\ \end{pmatrix}=\begin{pmatrix} 0 & 1 & 0 &0 \\ -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 \\ \varepsilon &0 &0 &0 \end{pmatrix}\begin{pmatrix} q\\ p \\ z\\ T\\ \end{pmatrix} \end{align} For initial conditions $q(0)=z(0)=x(0)=0$ and $p(0)=1$ we integrate numerically
Furhter Reading:
[1] J.M. Thijssen - Computational Physics
[2] Willian G. Hoover - Canonical dynamics: Equilibrium phase-space distibutions
[3] Denis J. Evans and Gary P. Morriss - Statistical Mechanics of NonEquilibrium Liquids (Page 69) [1]