Molecular Dynamics in Different Ensembles

• Abstract
• Molecular Dynamics in NVE
• Molecular dynamics in NPH

Abstract

The most important concept of molecular dynamics is that the ensemble average of a state function can be represented by the average of certain time-dependent trajectories. The procedure of molecular dynamics is to use classical mechanics to generate such trajectories, i.e. the equations of motion(EoM) of a series of $N$ particles in volume $V$ are solved numerically.

In this paper, Anderson expanded the molecular dynamics application from the traditional micro canonical ensemble (NVE) to three other ensembles: canonical (NVT) ensemble, isothermal-isobaric (NPT) ensemble and the isoenthalpic-isobaric (NPH) ensemble. All these four ensembles are equivalent for the calculation of the thermodynamics quantities. For intensive properties, $F_{NVE}$, $F_{NVT}$, $F_{NPT}$ and $F_{NPH}$ are equivalent except for the order of $N^{-1}$. And if $F$ is extensive property, all four should be the same.

Molecular Dynamics in NVE

The object of interest in Anderson’s paper is atomic fluid with a small number (between 50 and 1000) of $N$ particles, with coordinates $\textbf{r}_1, \textbf{r}_2,…,\textbf{r}_N$ in a cubic volume $V$. In real world, this physical system is more like a droplet where severe surface effect happens. A periodic boundary condition is employed to obtain the results for a bulk liquid: for particle $i$ at $\textbf{r}_i$, a series of virtual particles exists at position $\textbf{r}_i+\textbf{n}V^{1/3}$, where $\textbf{n}$ is a vector with integer components.

The Lagrangian for this system is:

$$\begin{equation} \mathcal{L}1(\textbf{r}^N,\dot{\textbf{r}}^N)= \frac{m}{2}\sum\limits_i\dot{\textbf{r}}i\cdot\dot{\textbf{r}}i- \sum\limits{i<j} u(\textbf{r}_{ij})\label{lagr1} \end{equation}\tag{1}$$

where $\sum\limits_{i<j} u(\textbf{r}_{ij})$ is the total potential energy of the atoms

The Hamiltonian is:

$$\begin{equation} \mathcal{H}1(\textbf{r}^N,\textbf{p}^N,V)=\sum\limits_i \dot{\textbf{r}}i\cdot\textbf{p}i-\mathcal{L}1 = (2m)^{-1}\sum\limits_i\textbf{p}i\cdot\textbf{p}i+\sum\limits{i<j}u(\textbf{r}{ij}) \end{equation}\tag{2}$$

Then we can build Lagrangian EoM for this system:

$$\begin{align} \frac{d\textbf{r}i}{dt} &= \frac{\partial \mathcal{H}}{\partial \textbf{p}i} = \frac{\textbf{p}i}{m}\label{EoM1}\tag{3}\\ \frac{d\textbf{p}i}{dt} &= -\frac{\partial \mathcal{H}}{\partial \textbf{r}i} = -\sum\limits{j(\ne i)}\hat{r}{ij}u^\prime(\textbf{r}{ij})\label{EoM2}\tag{4} \end{align}$$

In statistical thermodynamics, a thermodynamic property $F$ of the system can be expressed as a function $F(\textbf{r}^N,\textbf{p}^N,V)$. The micro-canonical ensemble average is denoted as $F_{NVE}(N, V, E)$:

$$\begin{equation} F_{NVE}(N, V, E)=[N!\Omega(N,V,E)]^{-1}\int_Vd\textbf{r}^N\int d\textbf{p}^N\delta[\mathcal{H}1(\textbf{r}^N,\textbf{p}^N,V)-E]F(\textbf{r}^N,\textbf{p}^N,V)\label{FNVE} \end{equation}\tag{5}$$

where

$$\begin{equation} \Omega(N,V,E) = (N!)^{-1} \int_V d \textbf{r}^N\int d\textbf{p}^N\delta(\mathcal{H}_1(\textbf{r}^N,\textbf{p}^N,V)-E)\label{PartitionNVE} \end{equation}\tag{6}$$

is the microcanonical ensemble partition function.

Then we can define the trajectory average to calculate the thermodynamics properties. For a function $F(\textbf{r}^N,\textbf{p}^N,V)$, the initial state $\textbf{r}^N(0)$ and $\textbf{p}^N(0)$ is chosen by initialize lattice coordinates and temperatures. The trajectory average of function $F(\textbf{r}^N,\textbf{p}^N,V)$ is defined as:

$$\begin{equation} \bar{F} = \lim_{T\rightarrow\infty}\int\limits_0^Tdt F(\textbf{r}^N,\textbf{p}^N,V)\label{tave} \end{equation}\tag{7}$$

When $N$ and $V$ are fixed, implement EoM in Eq.\ref{EoM}. The hypothesis is made that the trajectory spends equal time in every state(equal volume with the same value of energy). So that

$$\begin{equation} \bar{F}=F_{NVE}(N,V,E) \end{equation}\tag{8}$$

This is the theoretical foundation of the molecular dynamics. And Anderson made modification on the usual molecular dynamics to generate trajectories that can be used in $NPH$, $NVT$ and $NPT​$ ensembles.

Molecular dynamics in NPH

At constant pressure, the volume is no longer a fixed value with fluctuation present. In this case, the trajectory average will be the case of $NPH​$. Some alteration need to be made to accommodate to this CPMD situation.

Firstly, we replace the coordinates $\textbf{r}_1, \textbf{r}_2,…,\textbf{r}_N$ of atoms with a scaled coordinates$\textbf{$\rho$}_1, \textbf{$\rho$}_2,…,\textbf{$\rho$}_N$:

$$\begin{equation} \rho i=r i/V^{1/3},\ i = 1,2,...,N \end{equation}\tag{9}$$

So that the coordinates in the new system are no longer a fixed number but a quantities vary with volume $V​$ between 0 and 1. Thus we can call it scaled system.

Then we introduce such Lagrangian for this system:

$$\begin{equation} \mathcal{L}2(\rho,\dot{\rho}^N,Q,\dot{Q}) =\frac{1}{2} m Q^{2/3} \sum\limits{i=1}^N\dot{\rho}i\cdot\dot{\rho}i-\sum\limits{i<j=1}^N u(Q^{1/3}\rho{ij})+\frac{1}{2}M\dot{Q}^2-\alpha Q\label{lagr2} \end{equation}\tag{10}$$

Here $Q$ is a new “volume” and $\alpha$ is the “pressure”. Compare Eq.$\ref{lagr2}$ with Eq.$\ref{lagr1}$, we can find that the first two terms are just the same. The third term in Eq.$\ref{lagr2}$ is the kinetic energy of $Q$, and $\alpha Q$ is the potential energy of $Q$. Anderson gives the last two term a physical interpretation. Suppose the as-simulated fluid in a volume-variable container, whose volume can be changed by a piston. Thus $Q$ will be the coordinate of the piston and $\alpha V$ is the potential by an external pressure $\alpha$.

The generalized momentum conjugate to $\rho_i$ is :

$$\begin{equation} \pi_i=\frac{\partial \mathcal{L}2}{\partial \dot{\rho}i}=mQ^{2/3}\dot{\rho}_i \end{equation}\tag{11}$$

and which for $Q$ is:

$$\begin{equation} \Pi = \frac{\partial \mathcal{L}2}{\partial \dot{\rho}i} = m \dot{Q} \end{equation}\tag{12}$$

Thus the Hamiltonian is:

$$\begin{align} \mathcal{H}2(\rho^N,\pi^N,Q,\Pi)&=\sum\limits{i=1}^{N}\dot{\rho}i\cdot\pi_i+\dot{Q}\Pi-\mathcal{L}2(\rho^N,\dot{\rho}^N,Q,\dot{Q}) \nonumber \tag{13}\\ &=(2mQ^{2/3})^{-1}\sum\limits{i=1}^N\pi_i\cdot\pi_i+\sum\limits{i<j=1}^Nu(Q^{1/3}\rho_{ij})+(2M)^{-1}\Pi^2+\alpha Q\tag{14} \end{align}$$

So the EoM are:

$$\begin{align} \dot{\rho}i &=\frac{\partial \mathcal{H}}{\partial \pi_i} = \frac{\pi_i}{mQ^{2/3}}\label{typo}\tag{15}\\ \dot{\pi} &=-\frac{\partial \mathcal{H}}{\partial\rho_i} = -Q^{1/3}\sum\limits{j(\ne i)=1}\frac{u^\prime\rho{ij}}{|\rho{ij}|}\tag{16}\\ \dot{Q} &= \frac{\partial \mathcal{H}}{\partial \Pi} = \frac{\Pi}{W}\tag{17}\\ \dot{\Pi} &= \frac{\partial \mathcal{H}}{\partial Q} = -\frac{1}{3Q}(-\frac{1}{mQ^{2/3}}\sum\limits{i=1}^N\pi_i\cdot\pi_i+Q^{1/3}\sum\limits{i<j}\rho{ij}u^\prime(Q^{1/3\rho{ij}})+3\alpha Q)\tag{18} \end{align}$$

Actually, there’s one typo in Anderson’s paper in Eq.$\ref{typo}$, where $\dot{p}$ should be $\dot{\rho}$.

And there exist one-to-one correspondences between parameters in the scaled system and original system in phase space

$$\begin{align} V &= Q\label{corespondences1}\tag{19}\\ r_i &= Q^{1/3}\rho_i\label{corespondences2}\tag{20}\\ p_i &= \frac{\pi_i}{Q^{1/3}}\label{corespondences3}\tag{21} \end{align}$$

With EoM, we can generate the trajectory of the new system by $\rho(t)$, $\pi(t)$, $Q(t)$, $\Pi(t)$. And using the correspondences in phase space, we can generate the trajectories in original system by the trajectories in the scaled system. Add time dependence to Eq.$\ref{corespondences1}$-$\ref{corespondences3}$, we get:

$$\begin{align} V(t) &= Q(t)\tag{22}\\ r_i(t) &= Q^{1/3}(t)\rho_i(t)\tag{23}\\ p_i(t) &= \frac{\pi_i(t)}{Q^{1/3(t)}}\tag{24} \end{align}$$

Substitute them into the EoM of the scaled system:

$$\begin{align} \frac{dr_i}{dt} &= \frac{p_i}{m} + \frac{1}{3}r_i\frac{d\ln V}{dt}\label{EoMO1}\tag{25}\\ \frac{dp_i}{dt} &= -\sum\limits{j(\ne i)=1}^N\hat{r}{ij}u^\prime(r{ij})-\frac{1}{3}p_i\frac{d\ln V}{dt}\label{EoMO2}\tag{26}\\ \frac{Md^2V}{dt^2} &=-\alpha+(\frac{2}{3}\sum{i=1}^{N}\frac{p_i\cdot p_j}{2m}-\frac{1}{3}\sum\limits{i<j=1}^Nr{ij}u^\prime(r_{ij}))/V\label{EoMO3}\tag{27} \end{align}$$

Compared these equations with Eq.$\ref{EoM1}$-$\ref{EoM2}$, when mass of the piston $M$ is relatively large, so that $dV/dt=0$, i.e. volume is invariant, Eq. $\ref{EoMO1}$-$\ref{EoMO3}$ reduced to the original system.

Of course, we can calculate the thermodynamics properties in this new isoenthalpic-isobaric system by using the time average of the trajectories generated in this system, just like in $NVE$ system, i.e.:

$$\begin{equation} \bar{F} = F_{NPH}(N,\alpha,H)\label{conclusion} \end{equation}\tag{28}$$

The trajectory average $G(\rho^N,\pi^N,Q,\Pi)$ of the scaled system has similar form as Eq.$\ref{tave}$:

$$\begin{equation} \bar{G} = \lim_{T\rightarrow\infty} T^{-1}\int\limits_0^T dt G(\rho^N(t),\pi^N(t),Q(t),\Pi(t))\label{tavec2} \end{equation}\tag{29}$$

For any function $F(r^N, p^N, V)$, there is a corresponding $G(\rho^N, \pi^N, Q, \Pi)$ that

$$\begin{equation} G(\rho^N, \pi^N, Q, \Pi) = F(r^N, p^N, V) \end{equation}\tag{30}$$

Using Eq.$\ref{tave}$ and $\ref{tavec2}$, we can get:

$$\begin{equation} \bar{G} = \bar{F} \end{equation}\tag{31}$$

This time average $\bar{G}$ is equal to an ensemble average of $G$ over an $NE$ ensemble, i.e., an ensemble with fixed energy and fixed number of particles. This is just like a scaled system in microcanonical ensemble, as the coordinates $\rho_i$ are a dimensionless unit volume. So that:

$$\begin{equation} \bar{G} = G_{NE}(N,E) \end{equation}\tag{32}$$

where

$$\begin{equation} G_{NE}(N,E) = [N!\Omega_2(N,E)]^{-1}\int_1d\rho^N\int d\pi^N\int_0^{\infty} dQ\int_{-\infty}^{\infty}d\Pi \delta[\mathcal{H}_2(\rho^N,\pi^N,Q,\Pi)-E]G(\rho^N,\pi^N,Q,\Pi)\label{GNE} \end{equation}\tag{33}$$

the partition function

$$\begin{equation} \Omega_2(N,E)=(N!)^{-1}\int_1d\rho^N\int d\pi^N\int_0^{\infty} dQ\int_{-\infty}^{\infty}d\Pi \delta[\mathcal{H}2(\rho^N,\pi^N,Q,\Pi)-E]\label{ParNE} \end{equation}\tag{34}$$

Here the coordinate $\rho_i$ is integrated over the unit cube. Value $E$ comes from Eq.$\ref{tavec2}$, where $E$ is the constant energy of the trajectory calculated.

Similar to Eq.$\ref{FNVE}$, in $NPH$ ensemble, $F_{NPH}(N,\alpha,H)$ is defined as:

$$\begin{equation} F_{NPH}(N, P, H)=[N!\Gamma(N,P,H)]^{-1}\int_0^\infty dV\int_Vd\textbf{r}^N\int d\textbf{p}^N\delta[\mathcal{H}1(\textbf{r}^N,\textbf{p}^N,V)+PV-H]F(\textbf{r}^N,\textbf{p}^N,V)\label{FNPH} \end{equation}\tag{35}$$

where

$$\begin{equation} \Gamma(N,P,H)=(N!)^{-1}\int_0^\infty dV\int_Vd\textbf{r}^N\int d\textbf{p}^N\delta[\mathcal{H}_1(\textbf{r}^N,\textbf{p}^N,V)+PV-H]\label{ParNPH} \end{equation}\tag{36}$$

Compare Eq.$\ref{GNE}$ with Eq.$\ref{FNPH}$, we can find resemblance in the two equations, besides one integral over $d\Pi$. This term does not appear in the original system, but relatively important in the $NPH$ term. For $\bar{F} = G_{NE}$, using Eq.$\ref{FNPH}$ and Eq.$\ref{ParNPH}$, we can get

$$\begin{equation} \bar{F} = \frac{\int{-\infty}^{\infty}d\Pi\ \Gamma(N,\alpha,E-\Pi^2/2M)F{NPH}(N,\alpha,E-\Pi^2/2M)} {\int_{-\infty}^{\infty}d\Pi\ \Gamma(N,\alpha,E-\Pi^2/2M)} \end{equation}\tag{37}$$

The numerator can be expended in Taylor series of $\Pi^2/2M$:

$$\begin{equation} F_{NPH}(N,\alpha,E-\Pi^2/2M) = F_{NPH}(N,\alpha,E-\Pi^2/2M) - \frac{\Pi^2}{2M}\frac{\partial F_{NPH}(N,\alpha,E)}{\partial H} + O(N^{-2}) \end{equation}\tag{38}$$

If $F$ were an extensive property, the correction term would be the order of $N^{-1}$

Thus $\bar{F}$ can be expressed as:

$$\begin{equation} \bar{F} = F_{NPH}(N,\alpha,E)-\frac{\partial F_{NPH}(N,\alpha,E)}{\partial H}\frac{\int_{-\infty}^{\infty}d\Pi\ \Gamma(N,\alpha,E-\Pi^2/2M)(\Pi^2/2M)} {\int_{-\infty}^{\infty}d\Pi\ \Gamma(N,\alpha,E-\Pi^2/2M)} \end{equation}\tag{39}$$

The ratio of the two integral is $\overline{\Pi^2/2M}$. Therefore,

$$\begin{align} \bar{F} &= F_{NPH}(N,\alpha,E) - \overline{\Pi^2/2M}\partial F_{NPH}(N,\alpha,E)/\partial H +O(N^{-2})\nonumber\tag{40}\\ &= F_{NPH}(N,\alpha,E-\overline{\Pi^2/2M})+O(N^{-2}) \tag{41} \end{align}$$

So Eq.$\ref{conclusion}$ is proved. Enthalpy $H$ for this ensemble is $E-\overline{\Pi^2/2M}$, which is the energy of the trajectory in the scaled system, minus the average kinetic energy associated with the motion of Q.

---

• Next
  Norm-Conservation Pseudopotentials(NCPP)

---