Cayley modification for strongly stable path-integral and ring-polymer molecular dynamics

TL;DR

This paper introduces a Cayley transform-based modification to standard path-integral molecular dynamics algorithms, significantly improving numerical stability and ergodicity without increasing computational cost, thus enabling larger timesteps.

Contribution

The authors propose a Cayley transform-based modification that replaces the standard harmonic normal mode evolution in path-integral MD, enhancing stability and ergodicity.

Findings

Removes numerical artifacts in RPMD and TRPMD

Enables larger MD timesteps without loss of accuracy

Provides a simple, universally applicable modification

Abstract

Path-integral-based molecular dynamics (MD) simulations are widely used for the calculation of numerically exact quantum Boltzmann properties and approximate dynamical quantities. A nearly universal feature of MD numerical integration schemes for equations of motion based on imaginary-time path integrals is the use of harmonic normal modes for the exact evolution of the free ring-polymer positions and momenta. In this work, we demonstrate that this standard practice creates numerical artifacts. In the context of conservative (i.e., microcanonical) equations of motion, it leads to numerical instability. In the context of thermostatted (i.e., canonical) equations of motion, it leads to non-ergodicity of the sampling. These pathologies are generally proven to arise at integration timesteps that depend only on the system temperature and the number of ring-polymer beads, and they are numerically demonstrated for the cases of conventional ring-polymer molecular dynamics (RPMD) and thermostatted RPMD (TRPMD). Furthermore, it is demonstrated that these numerical artifacts are removed via replacement of the exact free ring-polymer evolution with a second-order approximation based on the Cayley transform. The Cayley modification introduced here can immediately be employed with almost every existing integration scheme for path-integral-based molecular dynamics - including path-integral MD (PIMD), RPMD, TRPMD, and centroid MD - providing strong symplectic stability and ergodicity to the numerical integration, at no penalty in terms of computational cost, algorithmic complexity, or accuracy of the overall MD timestep. Furthermore, it is shown that the improved numerical stability of the Cayley modification allows for the use of larger MD timesteps. We suspect that the Cayley modification will therefore find useful application in many future path-integral-based MD simulations.

Equations

$$e^{\Delta tL}=e^{\frac{\Delta t}{2}{L}_{\gamma }}{e}^{\frac{\Delta t}{2}{L}_V}{e}^{\Delta t{L}_0}{e}^{\frac{\Delta t}{2}{L}_V}{e}^{\frac{\Delta t}{2}{L}_{\gamma }}+\mathcal{O}(\Delta t^3)$$

$$e^{\Delta tL}=e^{\frac{\Delta t}{2}{L}_V}{e}^{\Delta t{L}_0}{e}^{\frac{\Delta t}{2}{L}_V}+\mathcal{O}(\Delta t^3).$$

$$\hat{H}=\frac{{\hat{p}}^2}{2m}+V(\hat{q})$$

$$Q=\mathrm{Tr}[e^{-\beta \hat{H}}],$$

$$Q_n=\frac{m^n}{(2\pi \hslash {)}^n}\int d^n\mathbf{q}\int d^n\mathbf{v}{e}^{-\beta H_n(\mathbf{q},\mathbf{v})},$$

$$H_n(\mathbf{q},\mathbf{v})=H_n^0(\mathbf{q},\mathbf{v})+V_n^{\text{ext}}(\mathbf{q}),$$

$$V_n^{\text{ext}}(\mathbf{q})=\frac{1}{n}\sum _{j=0}^{n-1}V(q_j)$$

$$H_n^0(\mathbf{q},\mathbf{v})=\frac{m_n}{2}\sum _{j=0}^{n-1}\left[v_j^2+{\omega }_n^2(q_{j+1}-q_j{)}^2\right],$$

$${\dot{q}}_j$$

$${\dot{v}}_j$$

$${\tilde{C}}_{qq}(t)=\frac{1}{Q}\mathrm{Tr}[e^{-\beta \hat{H}}\tilde{q}(0)\hat{q}(t)]$$

$$\tilde{q}=\frac{1}{\beta }{\int }_0^{\beta }{e}^{\lambda \hat{H}}\hat{q}{e}^{-\lambda \hat{H}}d\lambda$$

$${\tilde{C}}_{qq}(t)=\frac{1}{Q_n}\int d^n\mathbf{q}\int d^n\mathbf{v}{e}^{-\beta H_n(\mathbf{q},\mathbf{v})}\bar{q}(0)\bar{q}(t)$$

$$\bar{q}(t)=\frac{1}{n}\sum _{j=0}^{n-1}{q}_j(t),$$

$$\left[\begin{array}{c}\dot{\mathbf{q}}\\ \dot{\mathbf{v}}\end{array}\right]=\mathbf{A}\left[\begin{array}{c}\mathbf{q}\\ \mathbf{v}\end{array}\right]+\left[\begin{array}{c}0\\ \mathbf{F}(\mathbf{q})/m_n\end{array}\right],\text{where}\mathbf{A}=\left[\begin{array}{cc}0& \mathbf{I}\\ \mathbf{L}& 0\end{array}\right],$$

$$\mathbf{L}={\omega }_n^2\left[\begin{array}{cccccc}-2& 1& 0& \cdots & 0& 1\\ 1& -2& 1& 0& \cdots & 0\\ & \ddots & \ddots & \ddots & & \\ & & \ddots & \ddots & \ddots & \\ 0& \cdots & 0& 1& -2& 1\\ 1& 0& \cdots & 0& 1& -2\end{array}\right].$$

$$\mathbf{L}=-\mathbf{U}\mathbf{\Omega }{\mathbf{U}}^{\mathrm{T}},\text{where Ω=diag(ω02,…,ωn−12)}$$

$${\omega }_j^2=\{\begin{array}{ll}4{\omega }_n^2{\sin }^2\left(\frac{\pi j}{2n}\right)& \text{if j is even},\\ 4{\omega }_n^2{\sin }^2\left(\frac{\pi (j+1)}{2n}\right)& \text{else},\end{array}$$

$$\left[\begin{array}{c}\dot{q}\\ \dot{v}\end{array}\right]=\mathbf{A}\left[\begin{array}{c}q\\ v\end{array}\right]\text{where}\mathbf{A}=\left[\begin{array}{cc}0& 1\\ -{\omega }^2& 0\end{array}\right],$$

$${\mathbf{M}}_{\Delta t}=\left[\begin{array}{cc}1-\frac{\Delta t^2{\omega }^2}{2}& \Delta t\\ -\frac{1}{2}\Delta t{\omega }^2(2-\frac{\Delta t^2{\omega }^2}{2})& 1-\frac{\Delta t^2{\omega }^2}{2}\end{array}\right].$$

$$\Delta t=\frac{\pi k}{\omega }\text{for all k≥1}.$$

$$\mathrm{cay}(\Delta t\mathbf{A})\equiv (\mathbf{I}-(1/2)\Delta t\mathbf{A}{)}^{-1}(\mathbf{I}+(1/2)\Delta t\mathbf{A}).$$

$$\mathrm{cay}(\Delta t\mathbf{A}{)}^T\mathbf{J}\mathrm{cay}(\Delta t\mathbf{A})=\mathbf{J}\text{where}\mathbf{J}=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]$$

$$\left[\begin{array}{c}\dot{\mathbf{q}}\\ \dot{\mathbf{v}}\end{array}\right]=\mathbf{A}\left[\begin{array}{c}\mathbf{q}\\ \mathbf{v}\end{array}\right].$$

$$\Delta t=\frac{\pi k}{{\omega }_j}\text{for all k≥1 and 1≤j≤n−1}.$$

$$\mathrm{cay}(\Delta t\mathbf{A})=\left[\begin{array}{cc}\mathbf{U}& 0\\ 0& \mathbf{U}\end{array}\right]\mathrm{cay}\left(\Delta t\left[\begin{array}{cc}0& \mathbf{I}\\ -\mathbf{\Omega }& 0\end{array}\right]\right)\left[\begin{array}{cc}{\mathbf{U}}^{\mathrm{T}}& 0\\ 0& {\mathbf{U}}^{\mathrm{T}}\end{array}\right].$$

$$\begin{array}[]{lc}\textbf{Velocity half-step:}&\bm{v}\leftarrow\bm{v}+\frac{\Delta t}{2}\frac{\bm{F}}{m_{n}}\\ \textbf{Free ring-polymer step:\quad}&(\bm{q},\bm{v})\leftarrow\textrm{FRP}(\bm{q},\bm{v};\Delta t)\\ \textbf{Force evaluation:}&\bm{F}=-\nabla V^{\textrm{ext}}_{n}(\bm{q})\\ \textbf{Velocity half-step:}&\bm{v}\leftarrow\bm{v}+\frac{\Delta t}{2}\frac{\bm{F}}{m_{n}}\end{array}$$

$$\mathbf{\varrho }=\mathbf{U}\mathbf{q}\text{and}\mathbf{\phi }=\mathbf{U}\mathbf{v}$$

$$\left(\begin{array}{c}{\varrho }_j(t+\Delta t)\\ {\phi }_j(t+\Delta t)\end{array}\right)=\exp (\Delta t{\mathbf{A}}_j)\left(\begin{array}{c}{\varrho }_j(t)\\ {\phi }_j(t)\end{array}\right)$$

$${\mathbf{A}}_j=\left[\begin{array}{cc}0& 1\\ -{\omega }_j^2& 0\end{array}\right],$$

$$\left(\begin{array}{c}{\varrho }_j(t+\Delta t)\\ {\phi }_j(t+\Delta t)\end{array}\right)=\mathrm{cay}(\Delta t{\mathbf{A}}_j)\left(\begin{array}{c}{\varrho }_j(t)\\ {\phi }_j(t)\end{array}\right),$$

Full text

I Introduction

Feynman’s path-integral formulation of quantum statistical mechanics Feynman and Hibbs (1965) offers powerful and widely used strategies for including nuclear quantum effects in complex chemical systems. These strategies are based on the observation that the quantum Boltzmann statistical mechanics of a quantum system is exactly reproduced by the classical Boltzmann statistical mechanics of an isomorphic ring-polymer system.Chandler and Wolynes (1981) For the numerically exact calculation of quantum Boltzmann statistical properties, the classical Boltzmann distribution of the ring-polymer system can be sampled using Monte CarloCeperley (1995) (i.e., path-integral Monte Carlo, or PIMC) or molecular dynamicsParrinello and Rahman (1984) (PIMD). For the approximate calculation of dynamical quantities, such as reaction rates,Craig and Manolopoulos (2005a, b); Geva, Shi, and Voth (2001) diffusion coefficients,Miller III and Manolopoulos (2005a, b); Hone and Voth (2004) and absorption spectra,Miller III and Manolopoulos (2005b); Habershon, Fanourgakis, and Manolopoulos (2008); Kaczmarek, Shiga, and Marx (2009); Witt et al. (2009); Paesani and Voth (2010) the Newtonian dynamics of the classical isomorphic system can be numerically integrated as a model for the real-time quantum dynamics, as in ring-polymer molecular dynamics (RPMD)Craig and Manolopoulos (2004); Habershon et al. (2013) and centroid molecular dynamics (CMD).Cao and Voth (1994); Voth (1996) These and related methods have enjoyed broad applicability in recent years for exploring nuclear quantum effects in the domains that span physical, bio-, geo-, and materials chemistry. Markland and Ceriotti (2018)

For PIMD and RPMD calculations, considerable effort has been dedicated to the development and refinement of numerical integration schemes. This work falls into two distinct categories. In the first, the RPMD equations of motion are preconditioned by modifying the ring polymer mass matrix; this causes the integrated trajectories to differ from those of the RPMD model,Martyna, Hughes, and Tuckerman (1999); Minary, Martyna, and Tuckerman (2003); Beskos et al. (2008, 2011); Lu and Zhou (2018); Zhang et al. (2017); Liu, Li, and Liu (2016) but it can lead to efficient and strongly stableBeskos et al. (2008, 2011); Lu and Zhou (2018) sampling of the quantum Boltzmann distribution. In the second category, no modification is made to the ring-polymer mass matrix (i.e., the “physical” masses of the ring-polymer beads are employed).Braams and Manolopoulos (2006); Ceriotti et al. (2010); Ceriotti, Manolopoulos, and Parrinello (2011); Rossi, Ceriotti, and Manolopoulos (2014); Rossi, Kapil, and Ceriotti (2018)

Within the second category, it is common to apply a thermostat to the internal ring-polymer motions, with two primary aims: to more efficiently sample the quantum Boltzmann distribution,Ceriotti et al. (2010); Ceriotti, Manolopoulos, and Parrinello (2011); Zhang et al. (2017) or to avoid the “spurious resonance” artifact of the microcanonical (i.e., un-thermostatted) RPMD equations of motion in which internal ring-polymer modes mechanically couple to physical modes of the system.Rossi, Ceriotti, and Manolopoulos (2014); Rossi, Kapil, and Ceriotti (2018) PIMD and RPMD integration schemes in the second category (which preserve the RPMD model dynamics) typically employ a Trotter-like factorization of the time evolution operator.Tuckerman et al. (1993); Miller III and Manolopoulos (2005a, b); Ceriotti et al. (2010); Ceriotti, Manolopoulos, and Parrinello (2011); Rossi, Ceriotti, and Manolopoulos (2014); Rossi, Kapil, and Ceriotti (2018) For the example of thermostatted RPMD (TRPMD)Rossi, Ceriotti, and Manolopoulos (2014) using the generalized Langevin equation (GLE) thermostat,Ceriotti et al. (2010) the numerical integration is performed usingBussi, Donadio, and Parrinello (2007)

[TABLE]

where the Liouvillian $L=L_{V}+L_{0}+L_{\gamma}$ includes contributions from the physical potential, $L_{V}$, the purely harmonic free ring-polymer motion, $L_{0}$, and the friction and thermal noise, $L_{\gamma}$; note that the standard microcanonical RPMD numerical integration scheme is then recovered in the limit of zero coupling to the thermostat, such that Miller III and Manolopoulos (2005a)

[TABLE]

Standard practice in these RPMD and PIMD integration schemes is to exactly evolve the harmonic free ring-polymer dynamics associated with $\exp(\Delta tL_{0})$ using the uncoupled free ring-polymer normal modes.Tuckerman et al. (1993); Miller III and Manolopoulos (2005a, b); Ceriotti et al. (2010)

The first major conclusion of the current work is that any PIMD, RPMD, CMD,Stern and Berne (2001); Poulsen, Nyman, and Rossky (2003); Ceriotti et al. (2016); Tuckerman et al. (1993); Craig and Manolopoulos (2004); Habershon et al. (2013); Cao and Voth (1994); Parrinello and Rahman (1984); Hone, Rossky, and Voth (2006); Trenins, Willatt, and Althorpe (2019); Benson, Trenins, and Althorpe (2019); Ceriotti et al. (2010); Ceriotti, Manolopoulos, and Parrinello (2011); Rossi, Ceriotti, and Manolopoulos (2014); Rossi, Kapil, and Ceriotti (2018) or other integration scheme that involves the exact integration of the free ring polymer (i.e., involves the ubiquitous $\exp(\Delta tL_{0})$ step in terms of the ring-polymer normal modes) will exhibit provable numerical deficiencies, including resonance instabilities and non-ergodicity. For the case of the standard microcanonical RPMD integration scheme in Eq. 2, which is a symplectic map, exact evolution of the free ring-polymer step leads to the provable loss of strong symplectic stability and the demonstrable appearance of resonance instabilities in the integrated trajectories. For thermostatted RPMD and PIMD integration schemes that involve a free ring-polymer step,Ceriotti et al. (2010); Ceriotti, Manolopoulos, and Parrinello (2011); Rossi, Ceriotti, and Manolopoulos (2014); Rossi, Kapil, and Ceriotti (2018) exact evolution of that step leads to the provable and numerically demonstrable non-ergodicity.

The second major conclusion of the current work is that these numerical artifacts can be eliminated by simply replacing the exact evolution of the free ring polymer step with an approximation based on the Cayley transform: an alternative to exact free ring-polymer evolution that is no more costly, no more complicated, and no less accurate in the context of the full integration timestep. In particular, we show that this Cayley modification eliminates the resonance instabilities that occurs when trajectories are evolved using standard microcanonical RPMD integrators, and we show that it restores ergodicity to thermostatted RPMD and PIMD trajectories. Furthermore, we show that the improved numerical properties of the Cayley modification generally allows for larger RPMD and PIMD integration timesteps to be employed.

The paper is organized as follows. In section II we articulate the numerical instability problem in the context of standard RPMD numerical integration and introduce the Cayley modification as the solution. Section III numerically illustrates the instability of standard RPMD numerical integration and shows that the Cayley modification removes this problem. Finally, in section IV we generalize these findings to thermostatted trajectories.

II Theory

The theory introduced in this paper adapts and advances previous mathematical results on the numerical approximation of general second order Langevin stochastic partial differential equations with space-time white noise.Bou-Rabee (2017)

II.1 RPMD

We consider a quantum particle in 1D with Hamiltonian operator given by

[TABLE]

where $\hat{q}$, $\hat{p}$, and $m$ represent the particle position, momentum, and mass, respectively, and $V(\hat{q})$ is a potential energy surface. All results presented here are easily generalized to multiple dimensional quantum systems.

The thermal equilibrium properties of the system are described by the quantum mechanical Boltzmann partition function,

[TABLE]

where $\beta=(k_{B}T)^{-1}$ is the inverse temperature. Using a path-integral discretization, $Q$ can be approximated by a classical partition function $Q_{n}$ of a ring-polymer with $n$ beads,Parrinello and Rahman (1984)

[TABLE]

where $\bm{q}=(q_{0},\dots,q_{n-1})$ is the vector of bead positions, and $\bm{v}$ is the corresponding vector of velocities. The ring-polymer Hamiltonian is given by

[TABLE]

which includes contributions from the physical potential

[TABLE]

and the free ring-polymer Hamiltonian

[TABLE]

where $m_{n}=m/n$, $\omega_{n}=n/(\hbar\beta)$ and $q_{n}=q_{0}$.

If we let $n=1$ in Eq. 5, the classical partition function of the system (governed by a classical Hamiltonian, Eq. 6 with $n=1$) is recovered, i.e. $Q_{1}=Q_{cl}$. In the limit $n\rightarrow\infty$, the path-integral approximation converges to the exact quantum Boltzmann statistics for the system, such that $Q_{\infty}=Q$. The thermal ensemble of ring-polymer configurations associated with Eq. 5 can be sampled using either molecular dynamics (leading to PIMD methods) or Monte Carlo (leading to PIMC methods).

The classical equations of motion associated with the ring-polymer Hamiltonian in Eq. 6,

[TABLE]

yield the RPMD model for the real-time dynamics of the system.Craig and Manolopoulos (2004); Habershon et al. (2013) RPMD provides a means of approximately calculating Kubo-transformed thermal time-correlation functions, such as the position autocorrelation function

[TABLE]

where the Kubo-transformed position operator $\tilde{q}$ is

[TABLE]

and the time-evolved operator $\hat{q}(t)$ is $e^{i\hat{H}t/\hbar}\hat{q}e^{-i\hat{H}t/\hbar}$.

Specifically, the RPMD approximation to Eq. 10 is

[TABLE]

where $\bar{q}$ is the bead-averaged position

[TABLE]

and the pair $(\bm{q}(t),\bm{v}(t))$ are evolved by the RPMD equations of motion in Eq. 9 with initial conditions drawn from the classical Boltzmann-Gibbs measure.

The RPMD equations of motion can be compactly rewritten as

[TABLE]

$\bm{F}(\bm{q})=-\nabla V^{\textrm{ext}}_{n}(\bm{q})$, $\bm{I}$ is an $n\times n$ identity matrix, $\bm{0}$ is an array of zeros, and $\bm{L}$ is the $n\times n$ Toeplitz matrix

[TABLE]

We recognize $\bm{L}$ as the 1D discrete Laplacian endowed with periodic boundary conditions; it is negative semi-definite with spectral radius that scales as $n^{2}$, and since $\bm{L}$ is circulant, it can be diagonalized by the $n\times n$ real discrete Fourier transform (DFT) matrix. In particular, the spectral decomposition of $\bm{L}$ can be written as

[TABLE]

is a diagonal matrix of eigenvalues ordered in descending order and given by

[TABLE]

and $\bm{U}$ is an $n\times n$ matrix whose columns are the corresponding orthonormal eigenvectors.

In nontrivial applications, the RPMD equations of motion in Eq. 14 cannot be solved analytically. It is then necessary to employ approximate numerical integration of the equations of motion. As we discuss next, designing good numerical integrators for Eq. 14 is complicated by the interplay between the time-evolution of the free ring-polymer (obtained by setting $\bm{F}=0$ in Eq. 14) and the contributions from the physical forces $\bm{F}$.

II.2 Cayley removes instabilities in a free ring-polymer mode

RPMD is an example of highly oscillatory Hamiltonian dynamics. Petzold, Jay, and Yen (1997) To understand why numerical integration of such systems is tricky and why the Cayley modification is needed, it helps to consider the equations of motion for a particular normal mode of the free ring polymer with Matsubara frequency $\omega>0$:

[TABLE]

which are also the equations of motion for a linear oscillator with natural frequency $\omega$. If $\omega$ is large, Eq. 18 is highly oscillatory. Solving Eq. 18 amounts to approximating the matrix exponential $\operatorname{exp}(\Delta t\bm{A})$ where $\Delta t$ is a timestep size. A good $2\times 2$ matrix approximation $\bm{M}_{\Delta t}$ should satisfy:

(P1) Accuracy

$\|\bm{M}_{\Delta t}-\operatorname{exp}(\Delta t\bm{A})\|=O(\Delta t^{3})$.

(P2) Strong Stability

For all $\omega>0$, and for all $\Delta t$ smaller than some constant independent of $\omega$, $\bm{M}_{\Delta t}$ is a strongly stable symplectic matrix.

(P3) Time-Reversibility

For all $\omega>0$ and $\Delta t>0$, $\bm{M}_{\Delta t}$ is reversible with respect to the velocity flip matrix $\bm{R}=\begin{bmatrix}1&0\\ 0&-1\end{bmatrix}$, i.e., $\bm{R}\bm{M}_{\Delta t}\bm{R}=\bm{M}_{\Delta t}^{-1}$.

We briefly comment on each of these criteria for a good approximation. Property (P1) is a basic requirement that ensures second-order accuracy on finite-time intervals. Property (P3) is particularly useful for sampling from the stationary distribution, since a reversible map can be readily Metropolized Tierney (1998); Bou-Rabee and Vanden-Eijnden (2012); Bou-Rabee and Sanz-Serna (2018), and since time-reversibility in a volume-preserving numerical integrator leads to a doubling of the accuracy order (see Propositions 5.2 and Theorem 6.2 of Ref. Bou-Rabee and Sanz-Serna, 2018, respectively). Property (P2) is the most interesting. A symplectic matrix $\bm{S}$ is stable if all powers of the matrix $\bm{S}$ are bounded. A symplectic matrix $\bm{S}$ is strongly stable if $\bm{S}$ is stable and all sufficiently close symplectic matrices are also stable. In other words, $\bm{S}$ is strongly stable if there exists an $\epsilon>0$, such that all symplectic matrices $\bm{S}^{\epsilon}$ that are within a distance $\epsilon$ away from $\bm{S}$ are also stable. A sufficient condition for $\bm{S}$ to be strongly stable is that the eigenvalues of $\bm{S}$ are on the unit circle in the complex plane and are distinct; both the necessary and sufficient conditions for strong stability of symplectic matrices are known.Krein (1950)

Figure 1 illustrates the concept of strong stability. In particular, for different values of $t$ (as indicated in each panel), the black dots correspond to the eigenvalues of the symplectic matrix $\bm{S}=\operatorname{exp}(t\bm{A})$ with $\omega=3$, and the grey dots are the eigenvalues of a perturbation of $\bm{S}$ which preserves the symplectic nature of the matrix, specifically $\bm{S}^{\epsilon}=\operatorname{exp}((1/2)t\bm{B})\operatorname{exp}(t\bm{A})\operatorname{exp}((1/2)t\bm{B})$ where $\bm{B}=\begin{bmatrix}0&\epsilon\\ \epsilon&0\end{bmatrix}$ and $\epsilon=0.15$. For any $t$, note that the two eigenvalues of $\bm{S}$ are always on the unit circle, and hence, $\bm{S}$ is always stable, but as the figure shows, $\bm{S}$ is not always strongly stable. Indeed, in Figure 1 (a), we see that the two eigenvalues of $\bm{S}$, represented by a single black dot, are both equal to $(-1,0)$, which violates the condition for strong stability, and in this case, we see that one of the eigenvalues of $\bm{S}^{\epsilon}$ has modulus greater than unity, which implies that $\bm{S}^{\epsilon}$ is unstable. In Figure 1 (b), the two eigenvalues of $\bm{S}$ are distinct and equal to $(0,\pm 1)$, and hence, $\bm{S}$ is strongly stable. Since $\bm{S}$ is strongly stable, and $\epsilon$ is sufficiently small, $\bm{S}^{\epsilon}$ has eigenvalues that are on the unit circle, and hence, is itself stable. For a more detailed discussion of the concept of strong stability of symplectic matrices, see Section 42 of Ref. Arnold, 2013.

A natural candidate for an approximation $\bm{M}_{\Delta t}$ that satisfies these criteria is the Verlet integrator, which is ubiquitous in the classical simulation of molecular systems.Allen and Tildesley (1987); Frenkel and Smit (2002); Bou-Rabee (2014); Leimkuhler and Matthews (2015) For a single Matsubara frequency of the free ring polymer, the Verlet integrator gives

[TABLE]

However, for $\Delta t>2/\omega$, the eigenvalues of $\bm{M}_{\Delta t}$ are real and distinct, so that one of them has modulus $>1$, and therefore the powers of $\bm{M}_{\Delta t}$ grow exponentially. Thus, numerical stability requires $\Delta t<2/\omega$, and Verlet does not satisfy (P2), since this numerical stability requirement is not uniform in $\omega$.

Surprisingly, the exact solution for the normal-mode dynamics also does not satisfy (P2). To see why, note that the eigenvalues of the matrix exponential $\operatorname{exp}(\Delta t\bm{A})$ are $e^{\pm i\omega\Delta t}$ and (P2) requires that $e^{i\omega\Delta t}\neq e^{-i\omega\Delta t}$ which is violated if and only if

[TABLE]

At these timesteps, the exact solution violates strong stability. This is illustrated in Figure 2 (a), where the two eigenvalues of $\operatorname{exp}(\Delta t\bm{A})$ are plotted in the complex plane for a range of time-step sizes. Although the two eigenvalues of $\operatorname{exp}(\Delta t\bm{A})$ lie on the unit circle for all $\Delta t$, strong stability fails to hold whenever the eigenvalues are both equal to $(\pm 1,0)$.

A simple strategy to avoid these artificial resonances is to use a random timestep size $\delta t$, e.g., take as timestep size an exponential random variable $\delta t$ with mean $\Delta t$. Averaging $\operatorname{exp}(\delta t\bm{A})$ over the exponential probability density function yields $\bm{M}_{\Delta t}=\mathbb{E}(\operatorname{exp}(\delta t\bm{A}))=(\bm{I}-\Delta t\bm{A})^{-1}$, where here $\bm{I}$ is the $2\times 2$ identity matrix. Unfortunately, as can be easily verified, this matrix satisfies none of our criteria: it is neither symplectic, nor reversible, nor sufficiently accurate. However, we can easily turn this approximation into one that satisfies (P1), by simply composing $1/2$ step of this integrator with $1/2$ step of its adjoint $\bm{M}_{\Delta t}^{-1}$. This correction yields the Cayley transform of the matrix $\Delta t\bm{A}$,

[TABLE]

In fact, the Cayley transform satisfies all three of the specified criteria for a good numerical integrator. It is time-reversible since $\bm{R}\operatorname{cay}(\Delta t\bm{A})\bm{R}=(\bm{R}-(1/2)\Delta t\bm{R}\bm{A})^{-1}(\bm{R}+(1/2)\Delta t\bm{A}\bm{R})=\operatorname{cay}(\Delta t\bm{A})^{-1}$, where we used that $\bm{R}^{-1}=\bm{R}$. It is a symplectic matrix since

[TABLE]

where we used the fact that $\bm{A}$ is a Hamiltonian matrix (See Ref. Marsden and Ratiu, 1999, Section 2.5). More importantly, it is a strongly stable symplectic matrix for all $\Delta t>0$, as illustrated in Figure 2 (b); in contrast with the exponential map, for all $\omega>0$ and $\Delta t>0$ the eigenvalues of the Cayley map are $(4-\Delta t^{2}\omega^{2}\pm 4i\Delta t\omega)/(4+\Delta t^{2}\omega^{2})$, which are distinct and of unit modulus. Thus, not only is every matrix power of $\operatorname{cay}(\Delta t\bm{A})$ bounded, but the Cayley map is strongly stable uniformly in $\omega$ and $\Delta t$.

II.3 Cayley removes instabilities in microcanonical RPMD

For numerical integration of the conservative RPMD equations of motion (Eq. 9 or Eq. 14), it is standard practiceMiller III and Manolopoulos (2005a, b); Habershon et al. (2013) to employ a symmetrically split second-order integrator of the form in Eq. 2.

Furthermore, it is standard practice to exactly perform the free ring-polymer time evolution step,Habershon et al. (2013) using an exponential map of the form $\exp(\Delta tL_{0})=\operatorname{exp}(\Delta t\bm{A})$ where $\bm{A}$ is the matrix associated with the dynamics of the free ring-polymer Hamiltonian,

[TABLE]

In practice, the exact exponential map is executed by successively (i) changing from the Cartesian bead positions and velocities to the normal modes of the free ring polymer, (ii) numerically integrating each of the uncoupled normal mode equations of motion, and (iii) translating the time-evolved normal mode coordinates back into the Cartesian bead positions and velocities. Therefore, the numerical stability of standard RPMD numerical integration may be analyzed in normal mode coordinates, where the free ring-polymer equations of motion in Eq. 21 decouple into a system of $n$ independent oscillators with natural frequencies given by the eigenvalues of the matrix $\bm{L}$ in Eq. 17.

By applying Eq. 19 to each normal mode coordinate, we find that strong stability of the exact free ring-polymer time evolution is violated when

[TABLE]

Unstable pairs of $\Delta t$ and $n$ are plotted using solid lines in Fig. 3(b) for selected values of $j$ and $k$. The horizontal asymptotes in this figure reflect the fact that the eigenvalues of $\bm{L}$ converge to the eigenvalues of the continuous Laplacian endowed with periodic boundary conditions.

Unlike the exact free ring-polymer step used in standard RPMD numerical integration, the Cayley modification $\exp(\Delta tL_{0})\approxeq\operatorname{cay}(\Delta t\bm{A})$ is strongly stable for all $\Delta t>0$ uniformly in $n$. To see this, note that the Cayley transform can be equivalently computed in either bead or normal mode coordinates. More precisely, let $\bm{L}=-\bm{U}\bm{\Omega}\bm{U}^{\mathrm{T}}$ be the spectral decomposition of $\bm{L}$ given in Eq. 16. Direct computation then shows that

[TABLE]

Using this correspondence, one can invoke the preceding results on the one-dimensional oscillator, to conclude that $\operatorname{cay}(\Delta t\bm{A})$ is second-order accurate, strongly stable symplectic, and time-reversible.

Since the Cayley transform meets our criteria (P1)-(P3), and under suitable conditions on the force $\bm{F}$, the Cayley modification to the RPMD numerical integrator is provably stable and second-order accurate on finite-time intervals with a stability requirement that is uniform with respect to the number of ring polymer beads. On the other hand, standard RPMD integrators may display artificial resonance instabilities because the free RP step is not always strongly stable. These instabilities often manifest as exponential growth in energy when strong stability is lost, as will be discussed in Section III.

We emphasize that the improved numerical stability of the Cayley modification comes at zero cost in terms of algorithmic complexity or computational expense, and it preserves the same order of accuracy for the overall timestep. Use of this improved integration algorithm simply involves replacing the exact normal mode free ring-polymer step in the standard RPMD integrator with the Cayley modification.

II.4 Algorithmic comparison: Standard vs. Cayley

For complete clarity, we now present a side-by-side comparison of the full RPMD timestep (Eq. 2) with the free ring-polymer motion $\exp(\Delta tL_{0})$ implemented using either the standard exponential map (i.e., exact normal mode evolution) or via the Cayley modification. In both cases, the full RPMD timestep associated with the splitting in Eq. 2 is implemented using the algorithm

[TABLE]

In standard RPMD numerical integration, the free ring-polymer step is performed exactly, using:

Convert bead Cartesian coordinates to normal modes using the orthogonal transformation:

[TABLE]

where $\bm{U}$ is the real DFT matrix defined in Eq. 16. 2. 2.

From $t$ to $t+\Delta t$, exactly evolve the free ring polymer in the normal mode coordinates:

[TABLE]

where

[TABLE]

for $0\leq j\leq n-1$ with $\omega_{j}$ defined in Eq. 17. 3. 3.

Convert back to bead Cartesian coordinates using the inverse of $\bm{U}$, which is just its transpose, since $\bm{U}$ is orthogonal.

In the Cayley modification, the only change is to use the following in place of Eq. 25:

[TABLE]

where $\operatorname{cay}$ is the Cayley transform given in Eq. 20.

As a final algorithmic comparison, we note that another popular means of evolving the free ring-polymer involves multiple timestepping (MTS) with the reversible reference system propagator algorithm (RESPA) Tuckerman, Berne, and Martyna (1992); Martyna, Hughes, and Tuckerman (1999); Tuckerman et al. (1993), which introduces an inner loop of short timesteps. However, it is easily shown that MTS-RESPA can exhibit the same problem of resonance instabilities as exact normal-mode evolution, due to fact that MTS-RESPA also lacks the property of strong stability. Consequently, we will not further discuss MTS algorithms in the current work, although we recognize that combining the Cayley modification with MTS in the context of the ring-polymer contraction methodMarkland and Manolopoulos (2008a, b) is straightforward and worth pursuing.

III Results for RPMD

In this section, we demonstrate the numerical integration of the microcanonical RPMD equations of motions (Eq. 14). Specifically, we compare the performance of the standard RPMD integrator, which involves exact integration of the free ring-polymer modes (Eq. 25) and our refinement in which the Cayley modification is used (Eq. 26). Results are presented for simple one-dimensional potentials, including

[TABLE]

and using a mass of $m=1$.

We begin by numerically testing the conditions for loss of strong stability (Eq. 22) for the example of the harmonic potential (Eq. 27). Figure 3(a) shows a typical example of one of the approximately $25\%$ of trajectories that fail for the standard RPMD integration scheme with $\beta=1$, $n=16$, and $\Delta t=0.1$. The unstable trajectories start out with the typical values of ring-polymer energy in Eq. 6 (i.e. they are not the “hot” initial conditions from the tail of the thermal distribution), and they diverge to exponentially large energies after relatively short propagation time when run with the standard RPMD. All of these trajectories are stable when run with the Cayley modification.

The solid lines in Fig. 3(b) indicate predicted conditions for instability (Eq. 22). These analytical predictions are overlaid with a heatmap showing the fraction of stable RPMD trajectories on the harmonic potential using the standard RPMD integration scheme; for the purposes of the current section, a trajectories is deemed to be unstable if energy conservation associated with the ring-polymer Hamiltonian (Eq. 6) is violated by more than $10\%$ within 100 time units of simulation. There are clear correlations in Fig. 3(b) between the predicted instabilities and observed simulation results.

Finally, Fig. 3(c) presents the corresponding heatmap for the Cayley-modified RPMD integration scheme. The Cayley modification preserves the conditions for strong stability, and the only numerically unstable trajectories are found for extremely large timesteps ($\Delta t>0.6$). Comparison of Figs. 3(b) and (c) reveals the clear numerical advantages of the Cayley-modified RPMD integration scheme over the standard RPMD integration scheme.

Before proceeding, we emphasize the generality of the loss of strong stability with the standard RPMD numerical integrator: Eq. 22 makes no assumption with regard to the form of the physical potential, the dimensionality of the system, or the mass of the particles; it only depends on the temperature of the system and the number of ring-polymer beads in relation to the size of the integration timestep. Considering Eq. 22 for the $k=1$ index and the highest Matsubara frequency of the ring-polymer, it is straightforward to show that the smallest possible timestep $\Delta t_{*}$ at which strong stability is violated is given by

[TABLE]

We thus arrive at a highly practical expression for the “maximum safe timestep” that depends only on $\beta$ and $n$, such that all smaller timesteps avoid the loss of strong stability associated with Eq. 22. In Fig. 3(b), this result is plotted (dotted, black line) and seen to follow the convex hull of smallest timesteps created by the other curves. In passing, we note that if $\beta$ corresponds to room temperature and $n=64$, then the maximum safe timestep is $0.6$ fs, which is strikingly consistent with the conventional $0.5$ fs timestep employed in many PIMD simulations of liquid water.

Figure 4 confirms that the numerical instabilities of the standard RPMD integrator also manifest for anharmonic potentials. For both the weakly anharmonic (Eq. 28) and quartic (Eq. 29) potentials, we plot the fraction of stable trajectories as a function of timestep, comparing the standard RPMD integration scheme with the Cayley modification. Also shown are the fraction of stable classical mechanical trajectories (i.e., the 1-bead limit of RPMD) when integrated using the Verlet algorithm. Indeed, the standard RPMD integration scheme exhibits clear numerical instabilities at particular timesteps (which depend on the choice of $\beta$ and $n$), whereas the Cayley-modified integration scheme (like the classical integration scheme) avoids these pronounced instabilities.

For the results in Fig. 4, the maximum safe timestep is $\Delta t_{*}\approx 0.029.$ Note that the standard RPMD integration scheme on the weakly anharmonic potential does not exhibit significant loss of stability at this timestep, due to the fact that the unstable ring-polymer mode apparently does not sufficiently couple to the other modes on the timescale of the trajectories. However, the expected artifact at this timestep is indeed observed for the quartic potential. These results illustrate that the degree to which the resonance instabilities of standard RPMD integration manifest will depend on the application, but regardless of the system, these resonance instabilities can be removed using the Cayley modification. Finally, panel (c) in this figure compares the accuracy of the standard and Cayley-modified RPMD integration schemes for the case of the quartic oscillator, revealing that even with time-steps that three-fold exceed the maximum safe timestep of the standard integration scheme, the Cayley-modified scheme shows negligible loss of accuracy in the trajectories.

Figure 5 explores the degree to which the Cayley modification enables the use of larger timesteps in comparison to the standard RPMD integration scheme. Defining the “critical timestep” as the largest value of $\Delta t$ for which 980 out of 1000 trajectories are stable, we compare this quantity for standard and Cayley-modified RPMD numerical integration as a function of the number of ring-polymer beads; the trends in the figure are insensitive to the precise definition of the critical timestep. Also shown is the maximum safe timestep for the standard RPMD integration scheme (Eq. 30). The improved stability of the Cayley-modified integration scheme is seen to consistently allow for the use of larger RPMD timesteps. The numerical behavior of the standard RPMD integration scheme closely tracks the predictions of the maximum safe timestep, although as seen previously, the resonance instabilities do not always manifest on the timescale of the simulated trajectories. Interestingly, for small $n$ in the quartic-oscillator simulations, the standard RPMD integration scheme actually underperforms the prediction of the maximum safe timestep, given that it exhibits large energy fluctuations ($>10\%$) without fully encountering a resonance instability. In summary, using the maximum safe timestep for the standard RPMD integration scheme as a reference, the figure indicates that in these systems, the Cayley modification allows for substantial improvements in the allowed timestep size (three-fold or more for large $n$).

IV Results for TRPMD

Thermostatted RPMD (TRPMD) involves thermalization of the internal ring-polymer modes during RPMD dynamics, with the aims of improving sampling of the Boltzmann distributionCeriotti et al. (2010) or avoiding the “spurious resonance” artifact that can appear in RPMD simulations of vibrational spectra.Habershon, Fanourgakis, and Manolopoulos (2008); Rossi, Ceriotti, and Manolopoulos (2014) Following Refs. Ceriotti et al., 2010 and Rossi, Ceriotti, and Manolopoulos, 2014, we implement TRPMD using the splitting in Eq. 1, where $L_{T}$ corresponds to

[TABLE]

$\dot{\bm{W}}(t)$ is a white-noise vector (since $\bm{W}$ is an $n$-dimensional standard Brownian motion), and $\bm{\gamma}$ is an $n\times n$ friction matrix defined such that $\bm{U}^{\mathrm{T}}\bm{\gamma}\bm{U}$ is a diagonal matrix whose $k$th diagonal entry is equal to $\omega_{k}$ (Eq. 17). In normal mode coordinates (cf. Eq. 24), this thermostat is implemented by adding the following at the beginning and end of the full integration step outlined in Eq. 23:

[TABLE]

where $\xi_{j}$ is a standard normal variate.

IV.1 Cayley removes non-ergodicity in TRPMD

Given that it helps to avoid spurious resonances,Rossi, Ceriotti, and Manolopoulos (2014); Rossi, Kapil, and Ceriotti (2018) one might expect that a Langevin thermostat can also eliminate the instabilities we have observed in standard RPMD integrators. This turns out to be only partly true. Here, we show that (i) lack of strong stability in the free RP step induces non-ergodicity in standard TRPMD integrators, and (ii) the Cayley modification eliminates these non-ergodicity issues.

For this purpose, we revisit the simple case of a single free ring-polymer mode, as in Section II.2. Consider Eq. 18 with a Langevin thermostat,

[TABLE]

where $\gamma\geq 0$ is a friction factor and $\dot{W}(t)$ is a scalar white noise. The solution $(q(t),v(t))$ of Eq. 32 is a bivariate Gaussian with mean vector and covariance matrix given respectively by

[TABLE]

In the limit as $t\to\infty$, the probability distribution of $(q(t),v(t))$ converges to the classical Boltzmann-Gibbs measure, which in this case, is a bivariate normal distribution with mean vector and covariance matrix given respectively by

[TABLE]

In this situation, the standard TRPMD splitting in Eq. 1 inputs $(q_{0},v_{0})$ and outputs $(q_{1},v_{1})$ defined as

[TABLE]

where $\xi_{0},\eta_{0}$ are independent standard normal random variables, $\bm{E}=\operatorname{exp}(\Delta t\bm{A})$, and $\bm{O}$ is the $2\times 2$ matrix

[TABLE]

Moreover, the numerical solution after $N$ integration steps is a Gaussian vector with mean vector and covariance matrix given respectively by

[TABLE]

where

[TABLE]

From Eq. 19, if $\Delta t=k\pi/\omega$ for any $k\geq 1$, then $\bm{E}$ is not strongly stable. At these timesteps, the eigenvalues of the matrix $\bm{O}\bm{E}\bm{O}$ are given by $\lambda_{+}=(-1)^{k}$ and $\lambda_{-}=(-1)^{k}\exp(-k\pi\gamma/\omega)$. By the Cayley-Hamilton theorem for $2\times 2$ matrices,Andreescu (2016) we have the following representation of the $N$th power of $\bm{O}\bm{E}\bm{O}$

[TABLE]

Since $|\lambda_{+}|=1$, it follows from this representation that $\bm{\mu}_{N}$ does not converge to $\bm{\mu}$ in Eq. 34, since $\bm{\mu}_{N}$ clearly depends on the initial condition. Similarly, the covariance matrix $\bm{\Sigma}_{N}$ fails to converge to $\bm{\Sigma}$.

If we modify the above by replacing every instance of $\bm{E}$ with $\bm{C}=\operatorname{cay}(\Delta t\bm{A})$, the modified splitting is ergodic. More precisely, provided that the timestep is sufficiently small such that

[TABLE]

then the eigenvalues of $\bm{O}\bm{C}\bm{O}$ are a complex conjugate pair with complex modulus $|\lambda_{\pm}|=\exp(-\gamma\Delta t/2)$. Hence, the matrix $\bm{O}\bm{C}\bm{O}$ is asymptotically stable. Under condition 37, the Cayley-modified scheme converges to the exact classical Boltzmann-Gibbs measure, in this example.

These results carry over to TRPMD, where the free ring-polymer equations of motion in Eq. 21 decouple into a system of $n$ independent oscillators with natural frequencies given by the eigenvalues of the matrix $\bm{L}$ in Eq. 17. Although the analysis of TRPMD in this section was performed for the specific case of the splitting in Eq. 1 (i.e., the Bussi-Parrinello or OBABO splitting), we have confirmed that the same problem of non-ergodicity arises in the BAOAB splittingLeimkuhler and Matthews (2013) and can likewise be fixed via the Cayley modification.

IV.2 TRPMD numerical results

Figure 6 presents TRPMD results on the harmonic potential (Eq. 27) using $n=6$ and $\beta=1$. For a single TRPMD trajectory, we histogram the distribution of the normal mode coordinates that are sampled, employing the smallest timestep for which numerical instability is observed in the microcanonical case for this number of beads (see Fig. 3b); specifically, we use $\Delta t=0.26$, which corresponds to the instability condition in Eq. 22 for the case of $n=6$, $j=5$, and $k=1$. Using both standard and Cayley-modified TRPMD integration, the trajectory is sampled at every timestep for a total of 770 timesteps.

The centroid mode (panel a) follows harmonic motion, that is decoupled from the other degrees of freedom. With both integrators, the lower-frequency ($j=1-4$) internal ring-polymer modes are efficiently sampled and converge to the correct Gaussian distribution (panels c-f). However, the $j=5$ mode behaves qualitatively differently, as predicted by Eq. 22, with the standard TRPMD integrator showing clear non-ergodicity. The Cayley modification leads to ergodic sampling of all ring-polymer modes.

The lower frequency internal modes can also be afflicted with non-ergodicity at larger timesteps in this system. For the next-smallest unstable timestep in Fig. 3 ($\Delta t=0.3$, which corresponds to the instability condition in Eq. 22 with $j=3,4$, and $k=1$), the simulations were repeated. As predicted by the instability condition, modes 3 and 4 are found to be non-ergodic if sampled using the standard TRPMD integrator (Fig. 7); again, ergodicity is recovered using the Cayley modification. The same non-ergodicity problems appear for anharmonic potentials using the standard TRPMD integrator and can easily be avoided with use of the Cayley modification.

We emphasize that the TRPMD results presented here employ a white-noise thermostat; there are additional non-ergodicity problems for coloured-noise TPRMD which the Cayley modification is not expected to improve, since they likely arise from the attenuation of the thermostat across particular frequency bands.Ceriotti, Manolopoulos, and Parrinello (2011); Ceriotti and Markland (2013)

V Summary

Strong stability is a relevant – and under-appreciated – concept for path-integral-based molecular dynamics methods. Without strong stability, numerical integration schemes are prone to numerical instabilities in the microcanonical case and non-ergodicity in the canonical case. Fortunately, one can easily imbue existing integration schemes, including those for PIMD, RPMD, TRPMD, and many CMD methods, with strongly stability via the Cayley modification introduced here. This can be done without downside in terms of the computational cost, algorithmic complexity, or accuracy of the numerical integration scheme. The numerical results presented here suggest that this will have practical benefits for simulation studies, including improved stability, improved sampling efficiency, and improved efficiency via the use of larger MD timesteps.

While the Cayley transformation is familiar in the chemical physics literature in the context of the Crank-Nicolson propagatorCrank and Nicolson (1947) for wavepacket dynamics,Press et al. (2007); Judson et al. (1991) and real-time path integralsMa, Kouri, and Hoffman (1993) it has not to our knowledge been utilized for molecular dynamics, due to an under-appreciation of the property of strong stability. We conclude by noting that path-integral-based MD methods are far from unique in the physical sciences in exhibiting highly oscillatory dynamics, with other notable examples including Markov-Chain-Monte-Carlo-based Bayesian statistical inversion,Kaipio and Somersalo (2005); Dashti and Stuart (2017); Borggaard, Glatt-Holtz, and Krometis (2018) transition path sampling, Reznikoff and Vanden-Eijnden (2005); Pinski and Stuart (2010); Bolhuis et al. (2002); Miller III and Predescu (2007) stochastic wave equations,Newhall and Vanden-Eijnden (2017) Drude-oscillator models for many-body polarizability and dispersion,Jones et al. (2013); Tkatchenko et al. (2012); Odbadrakh and Jordan (2016) and Carr-Parrinello molecular dynamics.Car and Parrinello (1985) We anticipate that the Cayley modification introduced here may have similar advantages in these and other areas of application.

Acknowledgements.

We thank Jesús Sanz-Serna, Xuecheng Tao, and Eric Vanden-Eijnden for helpful discussions. N. B.-R. was supported in part by the National Science Foundation under Award No. DMS-1816378. R. K. and T. F. M. acknowledge support from the Department of Energy under Award No. DE-FOA-0001912 and the Office of Naval Research under Award No. N00014-10-1-0884.