Non-Reversibility of Molecular Dynamics Trajectories

TL;DR

This paper investigates the causes of non-reversibility in molecular dynamics trajectories, showing that rounding errors lead to irreversibility but do not significantly impact Hybrid Monte Carlo methods, supported by data from gauge theories.

Contribution

It provides an analysis of rounding error amplification in molecular dynamics and argues that non-reversibility is not a major concern for Hybrid Monte Carlo simulations.

Findings

Rounding errors cause non-reversibility in trajectories.

Non-reversibility does not significantly affect Hybrid Monte Carlo.

Data from SU(3) gauge theory and QCD supports the analysis.

Abstract

We study the non-reversibility of molecular dynamics trajectories arising from the amplification of rounding errors. We analyse the causes of such behaviour and give arguments, indicating that this does not pose a significant problem for Hybrid Monte Carlo computations. We present data for pure SU(3) gauge theory and for QCD with dynamical fermions on small lattices to illustrate and to support some of our ideas.

Equations

$$\Vert \Delta \delta U{\Vert }^2\equiv \sum _{x,\mu }\sum _{a,b}|U_{x,\mu }^{a,b}-U_{x,\mu }^{\prime a,b}{|}^2,$$

Full text

The theory of the Hybrid Monte Carlo (HMC) algorithm [1] assumes the exact reversibility of its molecular dynamics (MD) trajectories. Leapfrog integration guarantees this unless the initial conjugate gradient (CG) vector is chosen in time asymmetric way or finite precision arithmetic is used. While the first condition is easily ensured in practice by using a fixed starting vector for every CG inversion, all numerical computations carried out using floating point arithmetic are subject to rounding errors.

These rounding errors are normally not considered dangerous unless they are exponentially amplified. Indeed, without such an amplification, the time cost of reducing the error to some preset value grows only logarithmically with the number of arithmetic operations involved in the computation. This is a very small correction to the growth of the cost of the HMC algorithm as the volume and correlation length of the system are increased.

Exponential amplification will occur whenever nearby MD trajectories diverge from one another exponentially, i.e., when the MD evolution becomes unstable. There are two distinct mechanisms leading to such an instability [2]. First, this is typical for nonlinear equations in the chaotic regime. In fact, the existence of a positive leading Liapunov exponent for the MD equations of pure $SU(2)$ lattice gauge theory was proposed in Ref. [3]. The second possibility is that the result of the discrete integration scheme diverges exponentially from the true solution. This instability should grow with the number of integration steps and is thus expected to have characteristic time scale shorter than the one associated with intrinsic chaos. Our numerical results confirm this.

The integration instability can be analysed in the context of free field theory [2]. In fact, the behaviour of a single mode with frequency $\omega$ already reveals all the essential features. One can show that the instability accompanied by the exponentially decaying acceptance rate ($P_{\hbox{\rm\tiny acc}}\sim e^{-\nu\tau}$) occur when $\omega\delta\tau\geq 2$. Here $\tau$ and $\delta\tau$ are the trajectory length and the integration step size respectively. In Fig. 1 the ’$\sigma=0$’ line shows the characteristic exponent $\nu$ as a function of $\delta\tau$ with $\omega$ fixed to unity. Note the sharp “wall” arising at $\delta\tau=2$, where the instability sets in.

Qualitatively similar behaviour is observed for the case of many stable modes [2] . The onset of instability is determined by the highest frequency mode and occurs when $\omega_{max}\delta\tau=2$. In order to keep the acceptance rate constant for free field theory as the lattice volume $V\to\infty$, we must decrease $\delta\tau$ so that $V\delta\tau^{4}$ stays fixed. Consequently, the instabilities go away as we approach the thermodynamic limit. In this sense the leapfrog instability is a finite volume effect.

In interacting field theory the notion of independent modes loses its meaning. On the other hand, accepting the standard assumption that it can be still useful to think in these terms for asymptotically free field theories at short distances, it is quite plausible to expect similar scenario there too. The forces acting on the highest frequency mode due to the other modes will fluctuate in some complicated way however, and so

we expect that the “wall” at $\omega_{\max}\delta\tau=2$ will get smeared out. This is illustrated for the simple model of a harmonic oscillator whose frequency is randomly chosen from a Gaussian distribution with mean $\omega$ and standard deviation $\sigma$ before each MD step. The numerical results shown in Fig. 1 confirm that the “wall” in this model does indeed spread out.

Equipped with the above qualitative picture, we have studied reversibility numerically for $SU(3)$ gauge theory both in the pure gauge and dynamical fermion cases [2] (see also related work [4]). We evolved a typical equilibrium configuration ${U}$ using leapfrog equations for some time $\tau$, then reversed the momenta and evolved it again for the same amount of time to get the configuration ${U^{\prime}}$. Deviations from reversibility were measured by

[TABLE]

but we also recorded the change of energy at the end of the trajectory (${\delta H}$) and at the end of the reversed trajectory ($\Delta\delta H$).

In Fig. 2 we collected a typical set of data from one pure gauge configuration. The top and bottom graphs clearly show the integration instability “wall” at $\delta\tau\approx 0.6$, which has spread out as expected. At the same time the middle graph indicates that as we reach the “wall” ${\delta H}=O(10^{3})$, so the integration instabilities are of no practical importance for this system. Note however the case of the unreasonably long trajectory ($\tau=40$), where the reversibility is lost while ${\delta H}$ is very small implying a good acceptance rate.

When ploted as a function of $\tau$, all of our data show a clear exponential instability in $\|\Delta\delta U\|$. We extracted a characteristic exponent $\nu$ ($\|\Delta\delta U\|\sim e^{\nu\tau}$) and show the results in Fig. 3. Note the same qualitative behaviour we observed for the toy model in Fig. 1 except that the integration instability “wall” appears at different values of $\delta\tau$. This probably just reflects the different highest frequencies of these systems. In case of full QCD, the pseudofermions produce a force of the order of the inverse lightest fermionic mass thus giving the highest relevant frequency when simulating close to $\kappa_{c}$. This is reflected in the bottom graph where the integration instability appears at very small $\delta\tau$.

Notice also that the characteristic exponent does not approach zero for small $\delta\tau$, which confirms the existence of chaotic continuous time dynamics. Unlike the integration instability, the intrinsic chaos can not be controlled by adjusting $\delta\tau$. Moreover, accepting the standard hypothesis that the trajectory length should be scaled proportionally to the correlation length in order to reduce the critical slowing down, non-reversibility might cause problems when simulating closer to the continuum limit.

However, our numerical analysis indicates a strong $\beta-$dependence of the exponent $\nu$, characterizing the intrinsic chaos. Indeed, Fig. 4 shows this for $SU(3)$ pure gauge theory on $4^{4}$ and $8^{4}$ lattices. These results can be qualitatively understood if we hypothesize that chaos is not only a property of this continuous time evolution, but is also a property of the underlying continuum field theory. This would suggest that $\nu$ scales like a physical quantity. At small $\beta$ the lattice theory is in the strong coupling regime and does not obey the asymptotic scaling behaviour. At large $\beta$ the system is in a tiny box and is thus in the deconfined phase. The finite temperature phase transition at $N_{T}=4$ occurs near $\beta=5.7$, suggesting that the scaling region is in the vicinity of this value for our lattices. We have fitted our $8^{4}$ data at $\beta=5.4,5.5,5.6,5.7$ to the one loop asymptotic scaling form $\nu=ce^{-\beta/12\beta_{0}}$, with $\beta_{0}=(11-{2\over 3}n_{f})/16\pi^{2}$ and with constant $c$ being the only free parameter.

The resulting fit, shown in Fig. 4 is surprisingly good, suggesting that our hypothesis might indeed be correct. This would mean that the characteristic exponent is constant when measured in “physical” units, that is $\nu\xi$ would be constant as $\xi\to\infty$. If this is the case, then tuning the HMC algorithm by varying the trajectory length proportionally to the correlation length does not lead to any change in the amplification of rounding errors as we simulate closer to the continuum limit.