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. [9]. What we present here is essentially the first two parts of the second section of Reference [9], 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}}$$ (122)

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}}$$ (123)

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

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

which we solve directly to yield

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

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)$$ (126)

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$$ (127)

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)$$ (128)

Consider the action of the partial Liouville operator

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

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]$$ (130)

$$= \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]$$ (131)

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

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\righta... ...ght)\exp\left(iL_t2/P\right)\exp\left(iL_1t/2P\right)\right]^P$$ (133)

When $P$ is large but finite:

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

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)... ...t/2\right)\Gamma\left(t\right) = \Gamma\left(t+\Delta t\right)$$ (135)

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}}$$ (136)

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

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}{\pa... ...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... ...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{\... ... 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{\... ...\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) + ... ...left(\Delta t\right)^2}{2} \frac{{\bf F}\left[{\bf r}\left(0\right)\right]}{m},$$ (138)

$$\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\}$$ (139)

This is the velocity-Verlet algorithm, seen previously in Eqs 119-121. 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}}$$ (140)

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

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]$$ (142)

$${\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].$$ (143)

This is termed the position Verlet algorithm [9]. 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.