Liouvillean Split and Multiple Time Step Algorithm for the $N{\bf P}T$ Ensemble

We have seen in section 2.2 that the knowledge of the Liouvillean allows us to straightforwardly derive a multi-step integration algorithm. Thus, for simulation in the $N{\bf P}T$ ensemble, the Liouvillean $i \mathrm{L} = \Dot{\bf y} \nabla_{\bf y}$ is readily available from the equations of motion in (3.23-3.27). For sake of simplicity, to build our $N{\bf P}T$ multiple time step integrator we assume that the system potential contains only a fast intramolecular $V_{0}$ term and a slow intermolecular term $V_{1}$, as discussed in Sec. 2.3. Generalization to multiple intra and inter-molecular components is straightforward.

We define the following components of the $N{\bf P}T$ Liouvillean

$$i {\mathrm{L}}_x = - \sum_i {\bf P}_{i}\frac{p_{\eta}}{Q} {\mathbf \nabla}_{{\bf P}_... ... \frac{p_{\eta}}{Q} \left ( {\mathbf \nabla}_{{\bf P}_h} \right )_{\alpha\beta}$$ (3.56)

$$i {\mathrm{L}}_y = {\mathcal{F}}_{\eta} \nabla_{p_{\eta}}$$ (3.57)

$$i {\mathrm{L}}_z = \sum_i \overleftrightarrow{\bf G}^{-1}\dot \overleftrightarrow{\bf G}{\bf P}_{i}{\mathbf \nabla}_{{\bf P}_{i}}$$ (3.58)

$$i {\mathrm{L}}_u = \sum_{\alpha\beta} \left ( \boldmath\overleftrightarrow{\mathcal{... ...ght )_{\alpha\beta} \left ( {\mathbf \nabla}_{{\bf P}_h} \right )_{\alpha\beta}$$ (3.59)

$$i {\mathrm{L}}_{s} = \sum_{i} {\bf J}_{i}{\mathbf \nabla}_{{\bf P}_{i}} + \sum_{ik} {\... ...\alpha\beta} \left ( {\mathbf \nabla}_{{\bf P}_h} \right )_{\alpha\beta} \qquad$$ (3.60)

$$i \mathrm{G}_0 = \sum_i \frac{{\bf P}_{i}}{M_{i}} {\mathbf \nabla}_{{\bf S}_{i}} +... ...)_{\alpha\beta}}{W} \left ( {\mathbf \nabla}_{{\bf h}} \right )_{\alpha\beta} +$$

$$+ \frac{p_{\eta}}{Q} \nabla_{\eta} - \nabla_{{\bf l}_{ik}} V_{0} \nabla_{{\bf p}_{ik}},$$ (3.61)

where in Eq. (3.61) the scaled forces ${\bf F}_{i}$ have been replaced by its real space counterparts, i.e. ${\bf J}_{i}= {\bf h}^{-1} {\bf F}_{i}$.

The atomic scaling version of this Liouvillean breakup is derived on the basis of Eqs. (3.38). One obtains

$$i {\mathrm{L}}_x = - \sum_{ik} {\bf p}_{ik}\frac{p_{\eta}}{Q} {\mathbf \nabla}_{{\bf... ... \frac{p_{\eta}}{Q} \left ( {\mathbf \nabla}_{{\bf P}_h} \right )_{\alpha\beta}$$ (3.62)

$$i {\mathrm{L}}_y = {\mathcal{F}}_{\eta} \nabla_{p_{\eta}}$$ (3.63)

$$i {\mathrm{L}}_z = \sum_{ik} \overleftrightarrow{\bf G}^{-1}\dot \overleftrightarrow{\bf G}{\bf p}_{ik}{\mathbf \nabla}_{{\bf p}_{ik}}$$ (3.64)

$$i {\mathrm{L}}_u = \sum_{\alpha\beta} \left ( \boldmath\overleftrightarrow{\mathcal{... ...ght )_{\alpha\beta} \left ( {\mathbf \nabla}_{{\bf P}_h} \right )_{\alpha\beta}$$ (3.65)

$$i {\mathrm{L}}_{s} = \sum_{ik} {\bf j}_{ik} {\mathbf \nabla}_{{\bf p}_{ik}} + \sum_{\a... ...\alpha\beta} \left ( {\mathbf \nabla}_{{\bf P}_h} \right )_{\alpha\beta} \qquad$$ (3.66)

$$i \mathrm{G}_0 = \sum_{ik} \frac{{\bf p}_{ik}}{m_{ik}} {\mathbf \nabla}_{{\bf l}_{... ...)_{\alpha\beta}}{W} \left ( {\mathbf \nabla}_{{\bf h}} \right )_{\alpha\beta} +$$

$$+ \frac{p_{\eta}}{Q} \nabla_{\eta} - \nabla_{lik} V_{0} \nabla_{{\bf p}_{ik}} ,$$ (3.67)

where ${\bf j}_{ik} = {\bf h}^{-1} {\bf f}_{ik}$ and ${\cal V}, {\cal K}, {\cal F}_{\eta}$ are given in Eqs. (3.39,3.40). For the time scale breakup in the $N{\bf P}T$ ensemble we have the complication of the extra degrees of freedom whose time scale dynamics can be controlled by varying the parameter $Q$ and $W$. Large values of $Q$ and $W$ slow down the time dynamics of the barostat and thermostat coordinates. The potential $V$ determines the time scale of the $iG_{0}$ term (the fast component) and of the $iL_{s}$ contribution (the slow component). All other sub-Liouvilleans either handle the coupling of the true coordinates to the extra degrees of freedom ($iL_{x}$ expresses the coupling of all momenta (including barostat momenta) to the thermostat momentum, while $iL_{z}$ is a coupling term between the center of mass momenta and the barostat momentum), or drive the evolution of the extra coordinates of the barostat and thermostat ($iL_{y}$ and $iL_{u}$). The time scale dynamics of these terms depends not only on the potential subdivision and on the parameters $W$ and $Q$, but also on the type of scaling [27]. When the molecular scaling is adopted the dynamics of the virial term ${\boldmath\overleftrightarrow{\mathcal{V}}}$ contains contributions only from the intermolecular potential since the barostat is coupled only to the center of mass coordinates (see Eq. (3.29)). Indeed, the net force acting on the molecular center of mass is independent on the intramolecular potential, since the latter is invariant under rigid translation of the molecules. When atomic scaling or group (i.e. sub-molecular) scaling is adopted, the virial ${\boldmath\overleftrightarrow{\mathcal{V}}}$ (see Eq. (3.39)) depends also on the fast intramolecular such as stretching motions. In this case the time scale of the barostat coordinate is no longer slow, unless the parameter W is changed. For standard values of W, selected to obtain an efficient sampling of the $N{\bf P}T$ phase space [99,80], the barostat dependent Liouvilleans, Eqs. (3.60,3.59), have time scale dynamics comparable to that of the intramolecular Liouvillean $iG_{0}$ and therefore must be associated with this term^3.9.

Thus, the molecular split of the Liouvillean is hence given by

$$iL_{1} = iL_{x} + iL_y + iL_z + iL_u + iL_s$$

$$iL_{0} = iG_{0}$$ (3.68)

whereas the atomic split is

$$iL_{1} = iL_{x} + iL_y + iL_s$$

$$iL_{0} = iG_{0} + iL_z + iL_u$$ (3.69)

For both scaling, a simple Hermitian factorization of the total time propagator $e^{iLt}$ yields the double time discrete propagator

$$e^{iL_{1}+iL_{0}} = e^{iL_{1} \Delta t_{1} /2} (e^{iL_{0} \Delta t_{0}})^{n} e^{iL_{1} \Delta t_{1} /2}$$ (3.70)

where $\Delta t_{0}$, the small time step, must be selected according to the intramolecular time scale whereas $\Delta t_{1}$, the large time step, must be selected according to the time scale of the intermolecular motions. We already know that the propagator (3.71) cannot generate a symplectic. The alert reader may also have noticed that in this case the symmetric form of the multiple time step propagator Eq. (3.71) does not imply necessarily time reversibility. Some operators appearing in the definition of $L_{1}$ (e.g. $iL_z$ and $iL_{s}$) for the molecular scaling and in the definitions of $iL_{1}$ and $iL_{0}$ for the atomic scaling are in fact non commuting. We have seen in section 2.2 that first order approximation of non commuting propagators yields time irreversible algorithms. We can render the propagator in Eq. (3.71) time reversible by using second order symmetric approximant (i.e. Trotter approximation) for any two non commuting operators. For example in the case of the molecular scaling, when we propagate in Eq. (3.71) the slow propagator $e^{iL_{1}\Delta t/2}$ for half time step, we may use the following second order $O(\Delta t^{3})$ split

$$e^{iL_{1} \frac{\Delta t_{1}}{2}} \simeq e^{ i {\mathrm{L}}_y \frac{\Delta t_{1}}{4}} e^{ i {\mathrm{L}}_z... ...u \right ) \frac{\Delta t_{1}}{2}} e^{ i {\mathrm{L}}_x \frac{\Delta t_{1}}{4}}$$ (3.71)

An alternative simpler and equally accurate approach when dealing with non commuting operators is simply to preserve the unitarity by reversing the order of the operators in the first order factorization of the right and left operators of Eq. (3.71) without resorting to locally second order $O(\Delta t^3)$ approximation like in Eq. (3.72). Again for the molecular scaling, this is easily done by using the approximant

$$\left ( e^{iL_{1} \frac{\Delta t_{1}}{2}} \right)_{\rm left} = e^... ...elta t_{1}}{2}}e^ {iL_u \frac{\Delta t_{1}}{2}}e^ {iL_s \frac{\Delta t_{1}}{2}}$$ (3.72)

for the left propagator, and

$$\left (e^{iL_{1} \frac{\Delta t_{1}}{2}} \right )_{\rm right} = e... ...ta t_{1}}{2}} e^ {iL_y \frac{\Delta t_{1}}{2}} e^ {iL_x \frac{\Delta t_{1}}{2}}$$ (3.73)

for the rightmost propagator. Note that

$$\left ( e^{iL_{1} \frac{\Delta t_{1}}{2}} \right)_{\rm left}^{-1} = \left ( e^{iL_{1} \frac{\Delta t_{1}}{2}} \right)_{\rm right}$$ (3.74)

Inserting these approximations into (3.71) the overall integrator is found to be time reversible and second order. Time reversible integrators are in fact always even order and hence at least second order [71,68]. Therefore the overall molecular and atomic (or group) discrete time propagators are given by

$$e^{iL_{\rm mol}\Delta t_{1}} = e^ {iL_s \frac{\Delta t_{1}}{2}}e^ {iL_u \frac{\Delta t_{1}}{2}} ... ..._{1}}{2}} e^ {iL_x \frac{\Delta t_{1}}{2}} (e^{iG_{0} \Delta t_{0}})^{n} \times$$

$$\times e^{iL_x \frac{\Delta t_{1}}{2}} e^{iL_y \frac{\Delta t_{1}}{2}} e... ...elta t_{1}}{2}}e^ {iL_u \frac{\Delta t_{1}}{2}}e^ {iL_s \frac{\Delta t_{1}}{2}}$$ (3.75)

$$e^{iL_{\rm atom}\Delta t_{1}} = e^ {iL_s \frac{\Delta t_{1}}{2}} e^ {iL_y \frac{\Delta t_{1}}{2}} e^ {iL_x \frac{\Delta t_{1}}{2}} \times$$

$$\times (e^ {iL_u \Delta t_{0}/2} e^ {iL_z \Delta t_{0}/2} e^{iG_{0} \Delta t_{0}} e^ {iL_z \Delta t_{0}/2} e^ {iL_u \Delta t_{0}/2})^{n} \times$$

$$\times e^{iL_x \frac{\Delta t_{1}}{2}} e^{iL_y \frac{\Delta t_{1}}{2}} e^ {iL_s \frac{\Delta t_{1}}{2}}$$ (3.76)

The propagator $e^{iG_{0}\Delta t_{0}}$, defined in Eq. (3.62), is further split according to the usual velocity Verlet breakup of Eq. (2.22). Note that in case of molecular scaling the ``slow'' coordinates $({\bf S},~{\bf h},~\eta)$ move with constant velocity during the $n$ small times steps since there is no ``fast'' force acting on them in the inner integration. The explicit integration algorithm may be easily derived for the two propagators in Eqs. (3.76) and (3.77) using the rule in Eq. (2.23) and its generalization:

$$e^{a y \nabla_y} f \left ( y \right ) = f \left ( y e^{a} \right )$$

$$e^{\overleftrightarrow{\bf a} \mathbf{y}\nabla_{\mathbf{y}}} f \left ( \mathbf{y} \right ) = f \left ( \overleftrightarrow{e^{\mathbf{a}}} \mathbf{y} \right )$$ (3.77)

Where $a$ and $\overleftrightarrow{\bf a}$ are a scalar and a matrix, respectively. The exponential matrix $\overleftrightarrow{e^{\mathbf{a}}}$ on the right hand side of Eq. (3.78) is obtained by diagonalization of $\overleftrightarrow{\bf a}$.

As stated before the dynamics generated by Eqs. (3.23-3.27) or (3.35-3.38) in the $N{\bf P}T$ ensemble is not Hamiltonian and hence we cannot speak of symplectic integrators [100] for the $t$-flow's defined by Eqs. (3.76, 3.77). The symplectic condition Eq. (2.8) is violated at the level of the transformation (3.14-3.22) which is not canonical. However, the algorithms generated by Eqs. (3.76,3.77) are time reversible and second order like the velocity Verlet. Several recent studies have shown [26,25,27] that these integrators for the non microcanonical ensembles are also stable for long time trajectories, as in case of the symplectic integrators for the NVE ensemble.