Latest revision as of 14:49, 29 April 2013

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{T}}=\mathbf{G}(\mathbf{T})$ \begin{align} \dot{q}&=p\\ \dot{p}&=-q-z p\\ \dot{z}&=p^2 - T(q) - xz\\ \dot{x}&=z^2 - T(q) \end{align} with $T=1+ \varepsilon \arctan(q)$. The terms $p^2 - T(q)$ and $z^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)}_{=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 \\ \frac{-\varepsilon}{q^2+1} & 2p & -x & -z \\ \frac{-\varepsilon}{q^2+1} &0 &2z &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}\delta \dot q\\ \delta \dot p \\ \delta \dot z\\ \delta \dot T\\ \end{pmatrix}=\begin{pmatrix} 0 & 1 & 0 &0 \\ -1 & -z & -p & 0 \\ \frac{-\varepsilon}{q^2+1} & 2p & -x & -z \\ \frac{-\varepsilon}{q^2+1} &0 &2z &0 \end{pmatrix}\begin{pmatrix}\delta q\\ \delta p \\ \delta z\\ \delta 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. Summing of $i$ iteration steps we find the mean Lyapunov Exponent \begin{align} \lambda_j=\frac{1}{n \Delta t} \sum_i^n \ln \frac{ \mathbf{\delta}_j(t_i)}{\mathbf{\hat{\delta}}_j(t_{i-1})}=\frac{1}{n \Delta t} \sum_i^n \mathbf{\delta}_j(t_i) \end{align} where the last step is justified since we always start with a unit vector $|\mathbf{\hat{\delta}}_j(t_{i-1})|=1$ with length one.

Computing procedure

• Use the Runge Kutta method (4th order) to integrate the nonlinear system for intitial conditions $q(0)=z(0)=x(0)=0$ and $p(0)=1$, iteration steps of $\Delta t =0.001$ and iterate to get solutions for $t_0 \leq t \leq t_1$
• Plug in the values $\mathbf{T}(t)$ into the Jacobian $\mathbf{J}$ to get the linearised system at each point of the fiducial trajectory $\mathbf{G}(\mathbf{T})$
• choose four initial orthonormal pertubation vectors $\delta_1^{(0)}(t_0)=[1,0,0,0]$,$\delta_2^{(0)}(t_0)=[0,1,0,0]$,$\delta_3^{(0)}(t_0)=[0,0,1,0]$,$\delta_4^{(0)}(t_0)=[0,0,0,1]$ and integrate the linear system for each pertubation vector with the 4th order Runge Kutta method

The resulting vectors $\delta_1^{(0)}(t_1)$,$\delta_2^{(0)}(t_1)$,$\delta_3^{(0)}(t_1)$,$\delta_4^{(0)}(t_1)$ at time $t_1$ won't be orthogonal, but they can be orthogonalized by the Gram-Schmidt algorithm.

• keep the first vector

\begin{align}\delta_1^{(0)}(t_1)'=\delta_1^{(0)}(t_1)\end{align}

make the second vector

\begin{align}\delta_2^{(0)}(t_1)'=\delta_2^{(0)}(t_1)-\delta_2^{(0)}(t_1) \cdot \delta_1^{(0)}(t_1) \hat{\delta}_1^{(0)}(t_1)\end{align}

the thrid vector

\begin{align}\delta_3^{(0)}(t_1)'=\delta_3^{(0)}(t_1)-\delta_3^{(0)}(t_1) \cdot \delta_2^{(0)}(t_1) \hat{\delta}_2^{(0)}(t_1)-\delta_3^{(0)}(t_1) \cdot \delta_1^{(0)}(t_1) \hat{\delta}_1^{(0)}(t_1)\end{align}

and the fourth vector

\begin{align}\delta_4^{(0)}(t_1)'=\delta_4^{(0)}(t_1)-\delta_4^{(0)}(t_1) \cdot \delta_3^{(0)}(t_1) \hat{\delta}_3^{(0)}(t_1)-\delta_4^{(0)}(t_1) \cdot \delta_2^{(0)}(t_1) \hat{\delta}_2^{(0)}(t_1)-\delta_4^{(0)}(t_1) \cdot \delta_1^{(0)}(t_1) \hat{\delta}_1^{(0)}(t_1)\end{align}

• The first approximation of the i-th Lyapunov Exponent is

\begin{align} \lambda_i^{(1)} = \frac{1}{ \Delta t} \ln |\delta_i^{(0)}(t_1)'| \end{align}

the second approximation is the average

\begin{align} \lambda_i^{(2)} = \frac{1}{2 \Delta t} \ln |\delta_i^{(0)}(t_1)'|+\ln |\delta_i^{(1)}(t_2)'| \end{align}

and the $N$-th approximation is

\begin{align} \lambda_i^{(N)} = \frac{1}{n \Delta t} \sum_{k=1}^N \ln |\delta_i^{(k-1)}(t_k)'| \end{align}

after each approximation, the orthogonalized vectors are normalized and used as the new orthonormal pertubation vectors.

• For large $N$ this sum converges to the Lyapunov Exponent of the nonlinear system. The existence of the limit is ensured by the theorem of Oseledec

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]
[4] Wm. G. Hoover, C. G. Hoover, H. A. Posch - Dynamical Instabilities, Manifolds and Local Lyapunov Spectra Far From Equilibrium. Computational Methods in Science and Technology 2001 [2]
[5] Giancarlo Benettin et al. - Lyapunov Characteristic Exponents for Smooth Dynamical Systems and for Hamitonian Systems; A method for computing all of them [3]