Potential Subdivision and Multiple Time Steps Integrators for NVE Simulations

The ideas developed in the preceding sections can be used to build multiple time step integrators. Multiple time step integration is based on the concept of reference system. Let us now assume that the system potential $V$ be subdivided in $n$ terms such that

$$V = V_{0} + V_{1} + ... + V_{n}.$$ (2.26)

Additionally, we suppose that the corresponding average values of the square modulus of the forces $F_{k} = \vert \partial V_{k} / \partial {\bf x}\vert$ and of their time derivatives $\dot F_{k} = \vert d /dt (\partial V_{k} / \partial {\bf x}) \vert$ satisfy the following condition:

$$F_{0}^{2} >> F_{1}^{2} >> ..>> F_{n}^{2}$$

$$\dot {F}_{0}^{2} >> \dot {F}_{1}^{2} >> ..>>\dot {F}_{n}^{2} .$$ (2.27)

These equations express the situation where different time scales of the system correspond to different pieces of the potential. Thus, the Hamiltonian of the $k$-th reference system is defined as

$$H = T + V_{0} + .. V_{k},$$ (2.28)

with a perturbation given by:

$$P = V_{k+1} + V_{k+2} .. + .. V_{n}.$$ (2.29)

For a general subdivision of the kind given in Eq. (2.26) there exist $n$ reference nested systems. In the general case of a flexible molecular systems, we have fast degrees of freedom which are governed by the stretching, bending and torsional potentials and by slow intermolecular motions driven by the intermolecular potential. As we shall discuss with greater detail in section 4, in real systems there is no clear cut condition between intra and intermolecular motions since their time scales may well overlap in many cases. The conditions Eq. (2.27) are, hence, never fully met for any of all possible potential subdivisions.

Given a potential subdivision Eq. (2.26), we now show how a multi-step scheme can be built with the methods described in section 2.2. For the sake of simplicity, we subdivide the interaction potential of a molecular system into two components only: One intra molecular, $V_{0}$, generating mostly ``fast'' motions and the other intermolecular, $V_{1}$, driving slower degrees of freedom. Generalization of the forthcoming discussion to a $n$-fold subdivision, Eq. (2.26), is then straightforward.

For the 2-fold inter/intra subdivision, the system with Hamiltonian $H = T + V_{0}$ is called the intra-molecular reference system whereas $V_{1}$ is the intermolecular perturbation to the reference system. Correspondingly, the Liouvillean may be split as

$$iL_{0} = \dot q {\partial \over \partial q} - {\partial V_{0}\over \partial q} {\partial \over \partial p}$$

$$iL_{1} = - {\partial V_{1}\over \partial q} {\partial \over \partial p} .$$ (2.30)

Here $L_{0}$ is the Liouvillean of the 0-th reference system with Hamiltonian $T + V_{0}$, while $L_{1}$ is a perturbation Liouvillean. Let us now suppose now that $\Delta t_{1}$ is a good time discretization for the time evolution of the perturbation, that is for the slowly varying intermolecular potential. The discrete $e^{iL\Delta t_{1}}\equiv e^{(iL_{0}+iL_{1})\Delta t_{1}}$ time propagator can be factorized as

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

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

where we have used Eq. (2.21) and we have defined

$$\Delta t_{0} = {\Delta t_{1} \over n}$$ (2.32)

as the time step for the ``fast'' reference system with Hamiltonian $T + V_{0}$. The propagator (2.31) is unitary and hence time reversible. The external propagators depending on the Liouvillean $L_{1}$ acting on the state vectors define a symplectic mapping, as it can be easily proved by using Eq. (2.8). The full factorized propagator is therefore symplectic as long as the inner propagator is symplectic. The Liouvillean $iL_{0} \equiv \dot q {\partial / \partial q} - {\partial V_{0}/ \partial q} {\partial / \partial p}$ can be factorized according to the Verlet symplectic and reversible breakup described in the preceding section, but with an Hamiltonian $T + V_{0}$. Inserting the result into Eq. (2.31) and using the definition (2.30), the resulting double time step propagator is then

$$e^{iL\Delta t_{1}} = e^{{-\partial V_{1} \over \partial q} {\part... ...{-\partial V_{1} \over \partial q} {\partial \over \partial p} \Delta t_{1}/2 }$$ (2.33)

This propagator is unfolded straightforwardly using the rule (2.23) generating the following symplectic and reversible integrator from step $t=0$ to $t = \Delta t_1$:

$$\begin{array}{lll} p & \left ( {\Delta t_1 \over 2} \right )&... ... 2}) + F_{1} (n\Delta t_{0}) {\Delta t_{1} \over 2} \end{array}$$ (2.34)

Note that the slowly varying forces $F_{1}$ are felt only at the beginning and the end of the macro-step^2.2 $\Delta t_{1}$. In the inner $n$ steps loop the system moves only according to the Hamiltonian of the reference system $H = T + V_{0}$. When using the potential breakup, the inner reference system is rigorously conservative and the total energy of the reference system (i.e. $T + V_0 + \dots + V_k$) is conserved during the $P$ micro-steps.^2.3

The integration algorithm given an arbitrary subdivision of the interaction potential is now straightforward. For the general subdivision (2.26) the corresponding Liouvillean split is

$$iL_{0} = \dot q {\partial \over \partial q} - {\partial V_{i}\over \partia... ...; ~~ iL_{n} = - {\partial V_{k}\over \partial q} {\partial \over \partial p}.~~$$ (2.35)

We write the discrete time operator for the Liouville operator $iL = L_{0} + ... L_{n}$ and use repeatedly the Hermitian approximant and Trotter formula to get a hierarchy of nested reference systems propagator, viz.

$$e^{i(\sum_{i=0}^{n} L_{i}) \Delta t_{n}} = e^{iL_{n} {\Delta t_{n}\over 2}} \left ( e^{i(\sum_{i=0}^{n-1} L_{i})\Delta t_{n-1}} \right )^{P_{n-1}} e^{iL_{n} {\Delta t_{n}\over 2}}$$ (2.36)

$$~~~~~~~~~ \Delta t_{n} = \Delta t_{n-1} P_{n-1} ~~~~~~~~~~~~$$

$$e^{i(\sum_{i=0}^{n-1} L_{i}) \Delta t_{n-1}} = e^{iL_{n-1} {\Delta t_{n-1}\over 2}} \left ( e^{i(\sum_{i=0}^{n-2} L_{i})\Delta t_{n-2}} \right )^{P_{n-2}} e^{iL_{n-1} {\Delta t_{n-1}\over 2}}$$ (2.37)

$$~~~~~~~~~ \Delta t_{n-1} = \Delta t_{n-2} P_{n-2} ~~~~~~~~$$

$$....$$

$$e^{i(L_{0} + L_{1} +L_2) \Delta t_{2}} = e^{iL_{2} {\Delta t_{2}\over 2}} \left ( e^{i \left ( L_{1} + L_0 \right )\Delta t_{1}} \right )^{P_{1}} e^{iL_{2} {\Delta t_{2}\over 2}}$$ (2.38)

$$~~~~~~~~~ \Delta t_{2} = \Delta t_{1} P_{1} ~~~~~~~~~~~~~~~$$

$$e^{i(L_{0} + L_{1}) \Delta t_{1}} = e^{iL_{1} {\Delta t_{1}\over 2}} \left ( e^{iL_{0}\Delta t_{0}} \right )^{P_{0}} e^{iL_{1} {\Delta t_{1}\over 2}}$$ (2.39)

$$~~~~~~~~~ \Delta t_{1} = \Delta t_{0} P_{0} ~~~~~~~~~~~~~~~$$

where $\Delta t_{i}$ is the generic integration time steps selected according to the time scale of the $i$-th force $F_i$. We now substitute Eq. (2.39) into Eq. (2.38) and so on climbing the whole hierarchy until Eq. (2.36). The resulting multiple time steps symplectic and reversible propagator is then

$$\begin{array}{l} e^{iL\Delta t_{n}} = e^{F_{n}{\partial \over... ...}{\partial \over \partial p}{\Delta t_{n} \over 2}} \end{array}$$ (2.40)

The integration algorithm that can be derived from the above propagator was first proposed by Tuckerman, Martyna and Berne and called r-RESPA, reversible reference system propagation algorithm [20]