The Liouville Operator Formalism to Generating MD Integration Schemes

In this section, we present an elegant formalism for deriving MD integrators, as discussed by Tuckerman et al. [6]. What we present here is essentially the first two parts of the second section of Reference [6], including some of my own elaboration and some of that presented in section 4.3 of F&S.

Imagine a quantity $f$ which is a function of particle positions ${\bf r}^N$ and momenta ${\bf p}^N$. Its time derivative is given by

$$\dot{f} = \dot{\bf r}\frac{\partial f}{\partial{\bf r}} + \dot{\bf p}\frac{\partial f}{\partial {\bf p}}$$ (128)

We can write down a formal solution to this equation. First, define the Liouville operator as

$$iL = \dot{\bf r}\frac{\partial}{\partial{\bf r}} + \dot{\bf p}\frac{\partial}{\partial {\bf p}}$$ (129)

As Tuckerman points out, the $i$ is there by convention and ensures that the operator is Hermitian. We can re-express Eq. 128 as

$$\dot{f} = iLf$$ (130)

which we solve directly to yield

$$f\left(t\right) = \exp\left(iLt\right)f\left(0\right).$$ (131)

If $f$ is itself a vector quantity identical to the set of positions and momenta, $\Gamma$, we have a way to express, formally, the evolution of the system:

$$\Gamma\left(t\right) = U\left(t\right)\Gamma\left(0\right)$$ (132)

where $U(t) = \exp\left(iLt\right)$ is the classical propagator. The idea with numerical integration is that we find a way to represent the propagator as a discrete algorithm for constructing the system at some time $t + \Delta t$ given the system at time $t$.

Let's build our discrete integrator by decomposing the operator:

$$iL = iL_1 + iL_2$$ (133)

This does not necessarily lead to two independent propagators, because the two components do not commute; that is:

$$\exp\left[\left( iL_1 + iL_2\right)t\right] \ne \exp\left(iL_1t\right)\exp\left(iL_2t\right)$$ (134)

Consider the action of the partial Liouville operator

$$iL_1 \equiv \dot{\bf r}\left(0\right)\frac{\partial}{\partial{\bf r}},$$ (135)

which gives

$$f\left[{\bf p}^N\left(t\right),{\bf r}^N\left(t\right)\right] = \exp\left[t\dot{\bf r}\left(0\right)\frac{\partial}{\partial{\bf r}}\right] f\left[{\bf p}^N\left(0\right),{\bf r}^N\left(0\right)\right]$$ (136)

$$= \sum_{n=0}^{\infty} \frac{\left(\dot{\bf r}\left(0\right)t\right)... ...\partial{\bf r}^n}f\left[{\bf p}^N\left(0\right),{\bf r}^N\left(0\right)\right]$$ (137)

$$= f\left[{\bf p}^N\left(0\right),{\bf r}^N\left(0\right)+\dot{\bf r}\left(0\right)\right]$$ (138)

The last line is the collapse of the Taylor expansion of the line immediately above it. So, the effect of this operator fragment is a simple shift of coordinates given some initial velocities. This is an interesting fact: we can consider first-order integration as a Taylor expansion.

The next step of Tuckerman was to apply the Trotter identity:

$$\exp\left[\left( iL_1 + iL_2\right)t\right] = \lim_{P\rightarrow\... ...\left(iL_1t/2P\right)\exp\left(iL_t2/P\right)\exp\left(iL_1t/2P\right)\right]^P$$ (139)

When $P$ is large but finite:

$$\exp\left[\left(iL_1 + iL_2\right)t\right] = \left[\exp\left(iL_1... ...exp\left(iL_1t/2P\right)\right]^P\exp\left[\mathscr{O}\left(1/P^2\right)\right]$$ (140)

Now, we define a finite timestep as $\Delta t = t/P$ and we have then a discrete operator that, when applied to a configuration at time $t$, will produce the configuration at time $t + \Delta t$:

$$\exp\left(iL_1\Delta t/2\right)\exp\left(iL_2\Delta t\right)\exp\left(iL_1\Delta t/2\right)\Gamma\left(t\right) = \Gamma\left(t+\Delta t\right)$$ (141)

By performing this operation sequentially $P$ times, we recover a discretized version of the formal solution to generate $\Gamma\left(t\right)$ given $\Gamma\left(0\right)$.

Now we explicitly consider the decomposition:

$$iL_1 = \dot{\bf p}\left(0\right)\frac{\partial}{\partial{\bf p}}$$ (142)

$$iL_2 = \dot{\bf r}\left(0\right)\frac{\partial}{\partial{\bf r}}$$ (143)

We can perform one $\Delta t$'s worth of update using the following operation on $f$:

$$\exp\left(\dot{\bf p}\left(0\right)\left(\frac{\partial}{\partial... ...ta t}{2}\right) f\left[{\bf p}^N\left(t\right),{\bf r}^N\left(t\right)\right]$$

The action of the rightmost operator, $\exp\left(\dot{\bf p}\left(0\right)\left(\frac{\partial}{\partial{\bf p}}\right)\frac{\Delta t}{2}\right)$:

$$f\left[{\bf p}^N\left(t\right),{\bf r}^N\left(t\right)\right] \ri... ... p}\left(\Delta t\right)\right]^N,\left[{\bf r}\left(t\right)\right]^N\right\}$$

The action of the next rightmost, $\exp\left(\dot{\bf r}\left(0\right)\left(\frac{\partial}{\partial{\bf r}}\right)\Delta t\right)$:

$$f\left\{\left[{\bf p}\left(t\right)+\frac{\Delta t}{2}\dot{\bf p}... ...t\right)+ \Delta t \dot{\bf r}\left(\frac{\Delta t}{2}\right)\right]^N\right\}$$

Then, the action of the final operator:

$$f\left\{\left[{\bf p}\left(t\right)+\frac{\Delta t}{2}\dot{\bf p}... ...ight)+ \Delta t \dot{\bf r}\left(\frac{\Delta t}{2}\right)\right]^N\right\}\\$$

Noting that ${\bf p} = m\dot{\bf r}$ and ${\bf F} = m{\bf a} = \dot{\bf p}$, we can summarize the effect of this three-step update of the positions and velocities as

$${\bf r}\left(\Delta t\right) = {\bf r}\left(0\right) + \Delta t \dot{\bf r}\left(0\right) + \frac{\left(\Delta t\right)^2}{2} \frac{{\bf F}\left[{\bf r}\left(0\right)\right]}{m},$$ (144)

$$\dot{\bf r}\left(\Delta t\right) = \dot{\bf r}\left(0\right) + \frac{\Delta t}{2m} \left\{{\bf F}\le... ...left(0\right)\right] + {\bf F}\left[{\bf r}\left(\Delta t\right)\right]\right\}$$ (145)

This is the velocity-Verlet algorithm, seen previously in Eqs 125-127. Interestingly, we can also reverse the order of the decomposition; i.e.,

$$iL_1 = \dot{\bf r}\left(0\right)\frac{\partial}{\partial{\bf r}}$$ (146)

$$iL_2 = \dot{\bf p}\left(0\right)\frac{\partial}{\partial{\bf p}}$$ (147)

The update algorithm that arises is

$$\dot{\bf r}\left(\Delta t\right) = \dot{\bf r}\left(0\right) + \Delta t {\bf F}\left[{\bf r}\left(0\right) + \frac{\Delta t}{2m}\dot{\bf r}\left(0\right)\right]$$ (148)

$${\bf r}\left(\Delta t\right) = {\bf r}\left(0\right) + \frac{\Delta t}{2}\left[\dot{\bf x}\left(0\right) +\dot{\bf x}\left(\Delta t\right)\right].$$ (149)

This is termed the position Verlet algorithm [6]. Tuckerman et al. showed that this new algorithm results in a slightly lower drift in total energy in MD simulation of a simple Lennard-Jones fluid, when the time-step is greater than about 0.004.