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. In the so called Nosé Hoover Thermostat the hamiltonian does not only contain terms arising from the particles, but also an additional term. 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. It can be shown that such a system represents a Canonical Ensemble$^{1}$ (T,N,V const.). 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{x}&=p^2 - T(q) \end{align} with $T=1+ \varepsilon \arctan(q)$. The terms $p^2 - T(q)$ represent the thermosats, respectively. Notice, there are two thermostating terms, one would not be sufficent to produce dissipative behaviour.

Linearising the system$^3$ about an arbitrary point $\mathbf{T}_0$ gives \begin{align} \mathbf{G}(\mathbf{T}) \approx \underbrace{\mathbf{G}(\mathbf{T}_0)}_{=const.} + \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} Numerical integration for initial conditions $q(0)=z(0)=x(0)=0$ and $p(0)=1$ gives the values $\mathbf{x}(t)$ for $t_1,t_2,\dots$, for each point $t_i$ in the phase space we can find a linearised 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 & -z & -p & 0 \\ 0 & 2p & -x & -1 \\ \frac{\varepsilon}{q^2+1} &0 &0 &0 \end{pmatrix}\begin{pmatrix} q\\ p \\ z\\ T\\ \end{pmatrix} \end{align} To study the convergence or divergence of neigbouring trajectories we choose a set of basis (or tangent) vectors $\mathbf{\hat{\delta}}_1,\mathbf{\hat{\delta}}_2,\mathbf{\hat{\delta}}_3,\mathbf{\hat{\delta}}_4$ and integrate the linear system, starting from the point $\mathbf{x}(t_i)+\mathbf{\hat{\delta}}_j$ for each vector $\mathbf{\hat{\delta}}_j$. After the integration we have to apply Gram-Schmidt orthonomalisation. Then we can calculate the Lyapunov Exponent for one iteration step. \begin{align} \lambda_j=\frac{1}{n \Delta t} \sum_i \ln \frac{\mathbf{\hat{\delta}}_{ji}}{\mathbf{\hat{\delta}}}_j(i-1)} \end{align}

---

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]