The Parrinello-Rahman-Nosé Hamiltonian and the Equations of Motion

In order to derive the multiple time step integration algorithm using the Liouville formalism described in the preceding sections we must switch to the Hamiltonian formalism. Thus, we evaluate the conjugate momenta of the coordinates $S_{i \alpha}$, $l_{ik\alpha}$, $h_{\alpha \beta}$ and $s$ by taking the derivatives of the Lagrangian in Eq. (3.6) with respect to corresponding velocities, i.e.

$${\bf T}_{i} = M_{i}{\bf G } s^{2} \dot {\bf S}_i$$ (3.7)

$${\bf p}_{ik} = m_{ik} s^{2} \dot {\bf l}_{ik}$$ (3.8)

$${\bf P}_{h} = s^{2} W \dot {\bf h}$$ (3.9)

$$p_s = Q\dot s.$$ (3.10)

Where we have defined the symmetric matrix

$${\bf G } = {\bf h}^{t} {\bf h}$$ (3.11)

The Hamiltonian of the system is obtained using the usual Legendre transformation [72]

$$H(p,q) = \sum \dot q p - {\cal L}(q,\dot q) .$$ (3.12)

One obtains

$$H = {1 \over 2} \sum_{i}^{N} {\dot{{\bf T}_{i}} {\bf G}^{-1} \dot{{\b... ...{tr \left ( {\bf P}_h^t {\bf P}_h\right )\over s^{2} W} + {p_{s}^{2} \over 2 Q}$$

$$+ V + P\Omega + {g \ln s \over \beta}$$ (3.13)

In the extended systems formulation we always deal with real and virtual variables. The virtual variables in the Hamiltonian (3.13) are the scaled coordinates and momenta while the unscaled variables (e.g ${\bf R}_{i} = {\bf h} {\bf S}_{i}$ or $p_{ik\alpha}' = p_{ik\alpha}/s$ are the real counterpart. The variable $s$ in the Nosé formulation plays the role of a time scaling [86,80,94]. The above Hamiltonian is given in terms of virtual variables and in term of a virtual time and is indeed a true Hamiltonian function and has corresponding equation of motions that can be obtained applying Eq. (2.3) with ${\bf x} \equiv S_{i \alpha},l_{ik\alpha},h_{\alpha \beta},s,T_{i \alpha},p_{ik\alpha},\pi_{\alpha \beta},p_{s}$ in a standard fashion. Nonetheless, the equations of motions in terms of these virtual variable are inadequate for several reasons since for example one would deal with a fluctuating time step [86,94]. It is therefore convenient to work in terms of real momenta and real time. The real momenta are related to the virtual counterpart through the relations

$$T_{i \alpha} \rightarrow T_{i \alpha}/ s$$ (3.14)

$$p_{ik\alpha} \rightarrow p_{ik\alpha}/ s$$ (3.15)

$$\left( {\bf P}_h\right )_{\alpha\beta} \rightarrow \left( {\bf P}_h\right )_{\alpha\beta} / s$$ (3.16)

$$p_{s} \rightarrow p_{s} / s$$ (3.17)

It is also convenient [26] to introduce new center of mass momenta as

$${\bf P}_{i} \equiv {\bf G^{-1}}{\bf T}_{i} .$$ (3.18)

such that the corresponding velocities may be obtained directly without the knowledge of the ``coordinates'' ${\bf h}$ in ${\bf G}$^3.4, namely

$$\dot {\bf S}_{i} = {{\bf P}_{i} \over M}.$$ (3.19)

Finally, a real time formulation and a new dynamical variable $\eta$ are adopted:

$$t \rightarrow t/s$$ (3.20)

$$\eta \equiv \mathrm{ln} s$$ (3.21)

The equations of motions for the newly adopted set of dynamical variables are easily obtained from the true Hamiltonian in Eq. (3.13) and then using Eqs. (3.14-3.22) to rewrite the resulting equations in terms of the new momenta. In so doing, we obtain:

$$\dot {\bf l}_{ik} = \frac{{\bf p}_{ik}}{m_{ik}}, ~~~~~ \qquad \dot {\bf S}_{i}= \frac... ... \Dot{\bf h}= \frac{{\bf P}_h}{W}, ~~~~~ \qquad \Dot{\eta} = \frac{p_{\eta}}{Q}$$ (3.22)

$$\dot {\bf p}_{ik} = {\bf f}_{ik}^c - \frac{p_{\eta}}{Q}{\bf p}_{ik},$$ (3.23)

$$\dot {\bf P}_{i} = {\bf h}^{-1} {\bf F}_{i}- \overleftrightarrow{\bf G}^{-1}\dot \overleftrightarrow{\bf G}{\bf P}_{i}- \frac{p_{\eta}}{Q}{\bf P}_{i},$$ (3.24)

$$\dot {\overleftrightarrow{\bf P}}_h = \left ( {\boldmath\overleftrightarrow{\mathcal{V}}}+ \boldmath\ov... ...{K}}- {\bf h}^{-1} P_{ext} \det {\bf h}\right ) - \frac{p_{\eta}}{Q} {\bf P}_h,$$ (3.25)

$$\Dot{p}_{\eta} = {\mathcal{F}}_{\eta},$$ (3.26)

It can be verified that the conserved quantity ${\cal H}$ is associated with the above equations of motion, namely

$${\cal H} = \frac{1}{2} \sum_{i=1}^{N} \frac{ {\bf P}_{i}^{t} \overleftrighta... ...k}}{m_{ik}} + \frac{1}{2} \frac{ tr \left ( {\bf P}_h^t {\bf P}_h\right )}{W} +$$

$$+ \frac{1}{2} \frac{ p_{\eta}p_{\eta}}{Q} + V + P_{ext} \det {\bf h}+ g k_B T \eta.$$ (3.27)

The atomic force ${\bf f}_{ik}^{c} = {\partial V \over \partial r_{ik\alpha}} - \frac{ m_{ik}}{M_i} {\bf F}_{i}$ includes a constraint force contribution which guarantees that the center of mass in the intramolecular frame of the $l_{ik\alpha}$ coordinates remains at the origin. ${\boldmath\overleftrightarrow{\mathcal{V}}}$ and $\boldmath\overleftrightarrow{\mathcal{K}}$ are the virial and ideal gas contribution to the internal pressure tensor ${\bf P}_{int} = {\boldmath\overleftrightarrow{\mathcal{V}}}+ \boldmath\overleftrightarrow{\mathcal{K}}$ and they are defined as^3.5

$${\boldmath\overleftrightarrow{\mathcal{V}}} = \sum_{i=1}^{N} {\bf F}_{i}{\bf S}_{i}^{t}$$

$$\boldmath\overleftrightarrow{\mathcal{K}} = \sum_{i=1}^{N} M_i \left ( {\bf h}\dot {\bf S}_{i}\right ) \dot {\bf S}_{i}^{t}.$$ (3.28)

Finally ${\cal F}_{\eta}$ is the force driving the Nosé thermostat

$${\cal F}_{\eta} = \frac{1}{2} \sum_{i=1}^{N} M_i \dot {\bf S}_{i}... ...i=1}^{N} \sum_{k=1}^{n_i} \frac{ {\bf p}_{ik}^t {\bf p}_{ik}}{m_{ik}} - gk_B T,$$ (3.29)

with $g$ equal to the number of all degrees of freedom $N_{f}$ including those of the barostat^3.6.

Eqs. (3.14-3.21) define a generalized coordinates transformation of the kind of Eq. (2.4). This transformation is non canonical, i.e. the Jacobian matrix of the transformation from the virtual coordinates does not obey Eq. (2.8). This means that ${\cal H}$ in terms of the new coordinates Eq. (3.28) is ``only'' a constant of motion, but is no longer a true Hamiltonian: application of Eq. (2.1) does not lead to Eqs. (3.23-3.27). Simulations using the real variables are not Hamiltonian in nature in the sense that the phase space of the real variables is compressible [96] and that Liouville theorem is not satisfied [91]. This ``strangeness'' in the dynamics of the real variables in the extended systems does not of course imply that the sampling of the configurational real space is incorrect. To show this, it suffices to evaluate the partition function for a microcanonical distribution of the kind $\delta ({\cal H} - E)$, with ${\cal H}$ being given by Eq. (3.28). The Jacobian of the transformation of Eqs. (3.14-3.22) must be included in the integration with respect to the real coordinates when evaluating the partition function for the extended system. If the equations of motion in terms of the transformed coordinates are known, this Jacobian, ${\mathcal J}$, can be readily computed from the relation [73]:

$$\frac{d {\mathcal{J}}}{d t} = - {\mathcal{J}} \left ( \frac{\partial}{\partial {\bf y} } \cdot \Dot{\bf y} \right ).$$ (3.30)

Where ${\bf y}$ has the usual meaning of phase space vector containing all independent coordinates and momenta of the systems. Inserting the equations of motion of Eq. (3.27) into Eq. (3.31) and integrating by separation of variables yields

$${\cal J} = e^{ N_{f} \eta } \left [ \det {\bf h}\right ]^{6 N}.$$ (3.31)

Using (3.32) and integrating out the thermostat degrees of freedom, the partition function can be easily shown [91,97] to be equivalent to that that of $N{\bf P}T$ ensemble, i.e.

$$\Delta_{N{\bf P}T} \propto \int d{\bf h} e^{-\beta P_{text} \det( {\bf h})} Q ({\bf h})$$ (3.32)

with $Q({\bf h})$ being the canonical distribution of a system with cell of shape and size define by the columns of ${\bf h}$.^3.7