Isothermal-isotension ensemble

An extended ensemble, which was implemented in this project, is the isothermal-isotension ensemble with constant number of atoms (NtT). We study the behavior of vanadium containing self-interstitials and vacancies, at temperatures close to the melting point. In this case, one expects changes of the volume and shape of the sample due to the thermal expansion and point defects. Therefore, in order to describe the system accordingly we have to implement an algorithm, which allows for fluctuation of the shape and volume of the computational cell, and at the same time controls the pressure of the system.

Parinello and Rahman in 1980 [66] invented a new method, in which the shape and volume of the computational cell are dynamical variables. This technique is very helpful in study of the phase transitions. The statistical ensemble designed by Parrinello and Rahman, combined with Nose-Hoover thermostat is identified as an isothermal-isotension (NtT) ensemble.

The realization of this method is in the following way. First of all, we introduce scaled coordinates $\{s^i_\alpha \}$ in addition to the usual ones $\{r^i_\alpha \}$ :

$$s^i_\alpha=\sum_{\beta=1}^{N=3}H_{\alpha \beta} {r^i_\beta}$$ (11.1)

where the Greek indices $\alpha, \beta, ....$ are the coordinate indices $\{\alpha=1,2,3~ or~\alpha=x,y,z\}$, the Latin indices $i$ $\{~i=1,N \}$ count the atoms, and $H_{\alpha \beta}$ is a transformation matrix.

The volume of the computational box is given by:

$$V=det(H_{\alpha \beta })$$ (11.2)

One can introduce a metric tensor $G_{\alpha \beta}$ by using the $H_{\alpha \beta}$ matrix:

$$G_{\alpha \beta}=\sum_{\gamma=1}^{N=3}H_{\gamma\alpha }H_{\gamma \beta}$$ (11.3)

or

$$G=H^TH$$ (11.4)

where $T$ is stands for transpose. A distance between any two atoms are $i$ and $j$ is simply given by:

$$r_{ij}^2=\sum_{\alpha}\left ( r^i_\alpha-r^j_\alpha \right )... ...ha \right )G_{\alpha\beta} \left( s^i_\beta-s^j_\beta \right)$$ (11.5)

where $\alpha$ and $\beta$ are the coordinate indices.

The Lagrangian of the system has to be extended from 3N coordinate degrees of freedom to 3N+9, including the new 9 degrees of freedom of the real matrix $H_{\alpha \beta}$ :

$$L=T-U \longrightarrow L'=T+T_H -U -U_H$$ (11.6)

where the extra 'kinetic' energy $T_H$ and the corresponding 'potential' $V_H$ energy are added to the original Lagrangian to account for the new decrees of freedom.
The 'kinetic' energy term is given by:

$$T_H=\frac{1}{2}Q\sum_{\alpha=1}^{3}\sum_{\beta=1}^{3}\dot H_{\alpha \beta}\dot H_{\alpha \beta}$$ (11.7)

where $Q$ is a parameter, its physical meaning will be explained latter.
The 'potential' $U_H$ energy term is:

$$U_H=PV$$ (11.8)

here $V$ is volume of the system $V=det(H_{\alpha \beta})$ and $P$ its pressure, which can be obtained according to the virial theorem:

$$P=Nk_bT+\frac{1}{3}<\sum_{i,j}\vec r_{ij} \vec F_{ij}>$$ (11.9)

where $\vec F_{ij}$ is the force acting on atom $i$ due to the atom $j$, and $\vec r_{ij}=\vec r_i -\vec r_j$ is the vector which connects the atoms $i$ and $j$, and $< ... >$ stands for ensemble averaging.

The modified Lagrangian is written in the terms of new variables:

$$L'=\frac{1}{2}\sum_{i=1}^{3} \sum_{\alpha,\beta}^3m\dot s^i_... ...sum_{\beta=1}^{3}\dot H_{\alpha \beta}\dot H_{\alpha \beta}-PV$$ (11.10)

here, for the sake of simplicity we omit for a while the terms corresponding to the Noose-Hoover thermostat. The equations of motion are obtained in the usual way:

$$\frac {\partial L'}{\partial s^i_\alpha}-\frac{d}{dt}\frac {\partial L'}{\partial \dot s^i_\alpha}=0$$ (11.11)

$$\frac {\partial L'}{\partial h^i_\alpha}-\frac{d}{dt}\frac {\partial L'}{\partial \dot h^i_\alpha}=0$$ (11.12)

Therefore we obtain:

$$m \ddot s^i_\alpha=\sum_{\beta}(H^{-1}) _{\alpha \beta}f^i_\... ...(G^{-1})_{\alpha \beta}(\dot G)_{\beta \gamma}\dot s^i_\gamma$$ (11.13)

$$Q\ddot H_{\alpha \beta}=V\sum_{\gamma}(\sigma_{\alpha \gamma}-P\delta _{\alpha \gamma})(H^{-1}){\beta \gamma}$$ (11.14)

where $f^i_\beta=\sum_{j:~ j \neq i}F^{ij}_\beta$ and $m_1=m_2=...=m_i \equiv m$,
the $\sigma_{\alpha \gamma}$ is the stress tensor which describes anisotropy of the solid.

The stress tensor $\sigma_{\alpha \beta}$ is defined in the following way [67]:

$$\sigma_{\alpha \beta}=\frac{1}{V}\left (\sum_{i=1}^N \frac {... ...rt r_{ij}^\alpha r_{ij}^\beta }{\vert\vec r_{ij}\vert} \right)$$ (11.15)

here $v_i^\alpha$ is $\alpha$th component of the atom $i$ velocity,
and $r_i^\alpha$ is $\alpha$th component of the vector connecting atoms $i$ and $j$.
The hydrostatic pressure components are the diagonal elements of the stress tensor:

$$p_x=<\sigma_{xx}>;~~ p_y=<\sigma_{yy}>;~~ p_z=<\sigma_{zz}>~~$$ (11.16)

The stress tensor can be rewritten in the terms of the scaled variables $\{s^i_\alpha \}$:

$$\sigma_{\alpha \beta}=\frac{1}{V}\left (\sum_{i=1}^N m (\su... ...H_{\alpha \gamma} (s^i_\delta-s^j_\delta) F_{ij,\beta} \right)$$ (11.17)

here $F_{ij,\beta}$ is the $\beta$th component of the force acting on the atom $i$ due to the atom $j$.

The system is driven by dynamic imbalance between the applied external stress and the internally generated stress. The $Q$ parameter corresponds to the relaxation time of recovery from the imbalance between the external and the internal stress. It can be shown that the value of $Q$ does not influence the ensemble average [67].

The equations of motion generated by the Lagrangian can be effectively solved numerically. One has take into account that the Parinello-Rahman model is represented by a system of coupled ordinary differential equations for scaled coordinates $\{s^i_\alpha \}$ and for the $\{ H_{\alpha \beta } \}$ elements of the $H$ matrix. These equations have to be solved simultaneously.

This can be done using the predictor-corrector (PG) method, which includes the following steps:
1.) New scaled positions $\{ s^i_\alpha (t+\delta t) \}$ and velocities $\{ \dot s^i_\alpha (t+\delta t)\}$ are predicted on the basis of the old ones,
including the accelerations $\{ s^i_\alpha (t),\dot s^i_\alpha (t),\ddot s^i_\alpha (t) ,\ddot s^i_\alpha (t-\delta t) ,\ddot s^i_\alpha (t-2\delta t) \}$ :

$$s^i_\alpha(t+\delta t)=s^i_\alpha(t)+ \dot s^i_\alpha(t)\de... ... t)^2 \sum_{k=1}^3 a^{PG}_k \ddot s^i_\alpha(t+(1-k)\delta t)$$ (11.18)

$$\dot s^i_\alpha(t+\delta t)\delta t=s^i_\alpha(t+\delta t)-s... ... t)^2 \sum_{k=1}^3 b^{PG}_k \ddot s^i_\alpha(t+(1-k)\delta t)$$ (11.19)

here $\{ a^{PG}_k,~ b^{PG}_k \}$ are the (PG) prediction coefficients.
2.) The force acting on the atom $i$ is calculated using the predicted coordinates $\{ s^i_\alpha (t+\delta t) \}$:

$$\vec F_i= -\vec \nabla _iU_{pot}(r(t+\delta t)),~~ r^i_\alpha=\sum_{\beta=1}^3H_{\alpha \beta}s^i_\beta$$ (11.20)

3.) The temperature of the system and the stress tensor are evaluated at $t+\delta t$:

$$T=\frac {2}{3}\sum_{i=1}^N \sum_{\alpha=1}^3m (\dot r_\alpha^i)^2$$ (11.21)

$$\sigma_{\alpha \beta}=\frac{1}{V}\left (\sum_{i=1}^N m (\su... ...H_{\alpha \gamma} (s^i_\delta-s^j_\delta) F_{ij,\beta} \right)$$ (11.22)

4.) The values of the elements of the matrices $H_{\alpha \beta },\dot H_{\alpha \beta }$ are predicted:

$$H_{\alpha \beta}(t+\delta t)=H_{\alpha \beta}(t)+ \dot H_{\a... ...\sum_{k=1}^3 a^{PG}_k \ddot H_{\alpha \beta}(t+(1-k)\delta t)$$ (11.23)

$$\dot H_{\alpha \beta}(t+\delta t)\delta t=H_{\alpha \beta}(t... ...\sum_{k=1}^3 b^{PG}_k \ddot H_{\alpha \beta}(t+(1-k)\delta t)$$ (11.24)

5.) The values of the matrix $\ddot H_{\alpha \beta}$ is evaluated at $t+\delta t$ using the predicted values of the matrix $H_{\alpha \beta}$ according to:

$$\ddot H_{\alpha \beta}=\frac{V}{Q}\sum_{\gamma=1}^3(\sigma_{\alpha \gamma}-P\delta _{\alpha \gamma})(H^{-1}){\beta \gamma}$$ (11.25)

6.) The acceleration of of the atom $i$ is calculated:

$$\ddot s^i_\alpha(t+\delta t)=\sum_{\beta}(H^{-1})_{\alpha \b... ...m(G^{-1})_{\alpha \beta}(\dot G)_{\beta \gamma}\dot s^i_\gamma$$ (11.26)

if we take into account the Noose-Hoover thermostat:

$$\ddot s^i_\alpha(t+\delta t)=\sum_{\beta}(H^{-1})_{\alpha \b... ... s^i_\gamma - \frac{1}{M_s} \psi (t+\delta t) \dot s^i_\alpha$$ (11.27)

$$\psi(t+\delta t)=\psi(t)+ \frac{\delta t}{M_sNk_bT}\left [ ... ...2m (\dot r_i)^2}{3}-Nk_bT \right ]\dot s^i_\alpha,~~~\psi(0)=0$$ (11.28)

7.) Now the values of $H_{\alpha \beta}$ and $\dot H_{\alpha \beta}$ are corrected:

$$H_{\alpha \beta}(t+\delta t)=H_{\alpha \beta}(t)+ \dot H_{\a... ...\sum_{k=1}^3 c^{PG}_k \ddot H_{\alpha \beta}(t+(2-k)\delta t)$$ (11.29)

$$\dot H_{\alpha \beta}(t+\delta t)\delta t=H_{\alpha \beta} (... ...\sum_{k=1}^3 d^{PG}_k \ddot H_{\alpha \beta}(t+(2-k)\delta t)$$ (11.30)

here $\{ c^{PG}_k,~ d^{PG}_k \}$ are the (PG) correction coefficients.
8.) Finally, the cycle is completed by correction of positions $\{ s^i_\alpha (t+\delta t) \}$ and velocities $\{ \dot s^i_\alpha (t+\delta t)\}$ of all atoms:

$$s^i_\alpha(t+\delta t)=s^i_\alpha(t)+ \dot s^i_\alpha(t)\del... ... t)^2 \sum_{k=1}^3 c^{PG}_k \ddot s^i_\alpha(t+(2-k)\delta t)$$ (11.31)

$$\dot s^i_\alpha(t+\delta t)\delta t=s^i_\alpha(t+\delta t)- ... ...t)^2 \sum_{k=1}^3 d^{PG}_k \ddot s^i_\alpha(t+(2-k)\delta t)$$ (11.32)

9.) Go to the step 1.
The structural phase transitions can be studied using this algorithm. We used this algorithm to find the volume and the shape of the sample under the zero external pressure, close to the melting point $\bf T_m$.