Convergence with respect to imaginary-time discretization (i.e., the number of ring-polymer beads) is an essential part of any path-integral-based molecular dynamics (MD) calculation. However, an unfortunate property of existing non-preconditioned numerical integration schemes for path-integral molecular dynamics—including essentially all existing ring-polymer molecular dynamics (RPMD) and thermostatted RPMD (T-RPMD) methods—is that for a given MD time step, the overlap between the exact ring-polymer Boltzmann–Gibbs distribution and that sampled using MD becomes zero in the infinite-bead limit. This has clear implications for hybrid Metropolis Monte Carlo/MD sampling schemes, and it also causes the divergence with bead number of the primitive path-integral kinetic-energy expectation value when using standard RPMD or T-RPMD. We show that these and other problems can be avoided through the introduction of “dimension-free” numerical integration schemes for which the sampled ring-polymer position distribution has non-zero overlap with the exact distribution in the infinite-bead limit for the case of a harmonic potential. Most notably, we introduce the BCOCB integration scheme, which achieves dimension freedom via a particular symmetric splitting of the integration time step and a novel implementation of the Cayley modification [R. Korol *et al.*, J. Chem. Phys. **151**, 124103 (2019)] for the free ring-polymer half-steps. More generally, we show that dimension freedom can be achieved via mollification of the forces from the external physical potential. The dimension-free path-integral numerical integration schemes introduced here yield finite error bounds for a given MD time step, even as the number of beads is taken to infinity; these conclusions are proven for the case of a harmonic potential and borne out numerically for anharmonic systems that include liquid water. The numerical results for BCOCB are particularly striking, allowing for nearly three-fold increases in the stable time step for liquid water with respect to the Bussi–Parrinello (OBABO) and Leimkuhler (BAOAB) integrators, while introducing negligible errors in the calculated statistical properties and absorption spectrum. Importantly, the dimension-free, non-preconditioned integration schemes introduced here preserve ergodicity and global second-order accuracy, and they remain simple, black-box methods that avoid additional computational costs, tunable parameters, or system-specific implementations.

## I. INTRODUCTION

Considerable effort has been dedicated to the development of numerical integration schemes for imaginary-time path-integral molecular dynamics (PIMD).^{1} In comparison with standard classical molecular dynamics, PIMD numerical integration faces the additional challenge of the highly oscillatory dynamics of the ring-polymer internal modes. Work on PIMD numerical integration generally falls into two distinct categories. In the first category, the PIMD equations of motion are *preconditioned* by modifying the ring-polymer mass matrix;^{2–10} this approach, which includes the widely used staging algorithms,^{11} causes the integrated trajectories to differ from those of the ring-polymer molecular dynamics (RPMD) model for real-time dynamics,^{12,13} but it can lead to efficient^{4–6} sampling of the quantum Boltzmann distribution.^{14,15} In the second category, no modification is made to the ring-polymer mass matrix, i.e., the equations of motion are *non-preconditioned*.^{13,16–21}

With the aim of providing useful models for real-time quantum dynamics, as well as simple and efficient algorithms for equilibrium thermal sampling, the current work focuses on non-preconditioned PIMD numerical integration, notable examples of which include RPMD^{12,13} and its thermostatted variant T-RPMD.^{20} Numerical integration schemes for the latter methods typically employ symmetric factorizations of the time-evolution operator of the form^{11,16–25}

where the operator $L=A+B+O$ includes contributions from the purely harmonic free ring-polymer motion $A$, the external potential $B$, and a thermostat $O$. Note that the standard microcanonical RPMD numerical integration scheme is recovered in the limit of zero coupling to the thermostat and that Eq. (1) yields the “OBABO” scheme of Bussi and Parrinello^{22} when *a* = 1 and the “BAOAB” scheme of Leimkuhler^{25} when *a* = 0.

In our previous work,^{26} we emphasized that earlier PIMD numerical integration schemes had overlooked a fundamental aspect of the $exp((\Delta t/2)A)$ sub-step of the time evolution in Eq. (1). Standard practice in these integration schemes has been to exactly evolve the harmonic free ring-polymer dynamics associated with $exp((\Delta t/2)A)$ using the uncoupled free ring-polymer normal modes,^{11,16–18} which was shown to lack the property of strong stability in the numerical integration, leading to resonance instabilities for microcanonical RPMD and loss of ergodicity for T-RPMD.^{26} Use of the Cayley modification to the free ring-polymer motion was shown to impart strong stability to the time evolution, thereby improving numerical stability for microcanonical RPMD and restoring ergodicity for T-RPMD.^{26}

In the current study, we focus on the accuracy of both statistical and dynamical properties of the OBABO and BAOAB schemes, as well as the corresponding integrators obtained when the exact free ring-polymer step is replaced by the strongly stable Cayley modification (OBCBO and BCOCB, respectively). Particular attention is paid to the effect of finite-time step error with these integrators in the limit of large bead numbers. Of these four integrators, it is found that only BCOCB is “dimension-free,” in the sense that the sampled ring-polymer position distribution has non-zero overlap with the exact distribution in the infinite-bead limit for the case of a harmonic potential. It is further shown that the OBCBO scheme can be made dimension-free via the technique of force mollification. It is shown that the newly introduced BCOCB integrator yields better accuracy than all other considered non-preconditioned PIMD integrators and allows for substantially larger time steps in the calculation of both statistical and dynamical properties. Importantly, these gains are made without loss of computational efficiency or algorithmic simplicity.

## II. NON-PRECONDITIONED PIMD

Consider a one-dimensional molecular system with potential energy function *V*(*q*) and mass *m*. The equations of motion for the corresponding *n*-bead ring polymer held at constant temperature *T* by a Langevin thermostat are

Here, ** W** is an

*n*-dimensional standard Brownian motion;

**(**

*q**t*) = (

*q*

_{0}(

*t*), …,

*q*

_{n−1}(

*t*)) is the vector of positions for the

*n*ring-polymer beads at time

*t*≥ 0, and $v(t)$ are the corresponding velocities;

*m*

_{n}=

*m*/

*n*and $\beta =(kBT)\u22121$; and $F(q)=\u2212\u2207Vnext(q)$, where $Vnext$ is the contribution of the external potential,

Moreover, **Ω**^{2} is the following *n* × *n* symmetric positive semi-definite matrix:

where *κ*_{n} = *n*/(*ℏβ*). Note that **Ω** can be diagonalized by an *n* × *n* orthonormal real discrete Fourier transform matrix ** U** as follows:

where *ω*_{j,n} is the *j*th Matsubara frequency^{27} given by

Finally, the matrix **Γ** in Eq. (2) is typically an *n* × *n* symmetric positive semi-definite friction matrix of the form

where *γ*_{j} is the friction factor in the *j*th normal mode.

In RPMD and T-RPMD calculations, one is often interested in the dynamics of Eq. (2) with initial conditions drawn from the stationary distribution with non-normalized density exp(−*βH*_{n}(** q**, $v$)), where

*H*

_{n}(

**, $v$) is the ring-polymer Hamiltonian defined by**

*q*and $Hn0(q,v)=(1/2)mn|v|2+qT\Omega 2q$ is the free ring-polymer Hamiltonian.

The standard method for discretizing Eq. (2) is to use a symmetric splitting method of the form of Eq. (1) that consists of a combination of three types of sub-steps: (i) exact free ring-polymer evolution of time step *τ*,

where $A=0\u2003I\u2212\Omega 2\u20030$ is the Hamiltonian matrix associated with the free ring polymer, (ii) velocity updates of time step *τ* due to forces from the external potential,

and (iii) velocity updates of time step *τ* due to the thermostat,

where ** I** is the

*n*×

*n*identity matrix and

**is an**

*ξ**n*-dimensional vector whose components are independent, standard normal random variables. The acronyms OBABO and BAOAB indicate the order in which these sub-steps are applied, as indicated in Eq. (1) with

*a*= 1 or

*a*= 0, respectively.

In previous work,^{26} we showed that the matrix exponential for the free ring-polymer evolution in Eq. (9) is not a strongly stable symplectic matrix, and as a consequence, the OBABO and BAOAB schemes can display non-ergodicity at time steps Δ*t* = *kπ*/*ω*_{j,n} for any 1 ≤ *j* ≤ *n* and *k* ≥ 1. We also identified a maximum safe time step size Δ*t*_{⋆} = *βℏπ*/(2*n*), below which the matrix exponential is strongly stable. As *n* → *∞*, this maximum safe time step goes to zero, such that no finite time step for the scheme in Eq. (1) is safe in this limit from non-ergodicity.

This non-ergodicity motivates the Cayley modification^{26} which consists of approximating the matrix exponential appearing in Eq. (9) with the Cayley transform. Specifically, for the Cayley-modified OBABO scheme (called OBCBO), we replace the exact free ring-polymer update of time step *τ* = Δ*t* with

For the Cayley-modified BAOAB scheme (called BCOCB), we replace the two exact free-ring polymer updates of half-time step *τ* = Δ*t*/2 with cay(Δ*t*** A**)

^{1/2}. While it might be expected that these half-time step updates would, instead, be replaced with cay((Δ

*t*/2)

**), such a choice leads to a loss of strong stability. Our use of the square root of the Cayley transform preserves strong stability, symplecticity, time reversibility, local third-order accuracy, and by definition satisfies cay(Δ**

*A**t*

**)**

*A*^{1/2}cay(Δ

*t*

**)**

*A*^{1/2}= cay(Δ

*t*

**). Furthermore, the square root of the Cayley transform is no more complicated to evaluate than the Cayley transform itself. Both the OBCBO and BCOCB Cayley modifications of Eq. (1) are ergodic for a fixed time step, irrespective of the number of beads; moreover, like Eq. (1), the Cayley-modified integrators exhibit locally third-order accuracy in the time step and leave invariant the free ring-polymer Boltzmann–Gibbs distribution in the special case of a constant external potential (**

*A**V*≡ const.).

^{26}

## III. BCOCB AVOIDS PATHOLOGIES IN THE INFINITE BEAD LIMIT

In this section, we show that of the OBABO, BAOAB, OBCBO, and BCOCB integration schemes, only BCOCB is dimension-free. Although the current section presents analytical results for the specific case of a harmonic external potential, these results are supported by numerical results for anharmonic external potentials in Secs. VI and VII.

To this end, consider the *j*th internal ring-polymer mode with frequency *ω*_{j,n}, in the presence of a harmonic external potential *V*(*q*) = (1/2)Λ*q*^{2} and a Langevin thermostat with friction *γ*_{j}. Expressed in terms of the normal mode coordinates, obtained from the Cartesian positions and velocities via the orthogonal transformation

where ** U** is defined in Eq. (5), the non-preconditioned PIMD equations of motion for this mode are

where $W\u0307j$ is a scalar white-noise, and we have introduced the following 2 × 2 matrices:

The solution (*ϱ*_{j}(*t*), *φ*_{j}(*t*)) of Eq. (14) is a bivariate Gaussian, and in the limit as *t* → *∞*, the probability distribution of (*ϱ*_{j}(*t*), *φ*_{j}(*t*)) converges to a centered bivariate normal distribution with the covariance matrix,

For this system, a single time step of Eq. (1) can be compactly written as

where *ξ*_{0} and *η*_{0} are independent standard normal random variables, and we have introduced the following 2 × 2 matrices:

where $P=0001$ and $Nj=ea\Delta t2Oje\Delta t2Be\Delta t2Aj$. The corresponding step for the Cayley modification is obtained by replacing exp((Δ*t*/2)*A*_{j}) in Eq. (16) with $cay(\Delta tAj)1/2$, which is given by

where

For the Cayley modification of Eq. (16), Eq. (18) still provides a sufficient condition for ergodicity, except with

Due to the lack of strong stability in the exact free ring-polymer evolution, Eq. (16) fails to meet the condition in Eq. (18) and becomes non-ergodic whenever Δ*t* = *kπ*/*ω*_{j,n}, where *k* ≥ 1;^{26} no such problem exists for the Cayley modification. Regardless, assuming that the condition in Eq. (18) holds, the numerical stationary distribution is a centered Gaussian with a 2 × 2 covariance matrix **Σ**_{j,Δt} that satisfies the linear equation,

for which the solution is

where the variance in the position and velocity marginal are $(\beta mn)\u22121sj,\Delta t2$ and $(\beta mn)\u22121rj,\Delta t2$ with

For the Cayley modification of Eq. (16),

Note that these numerical stationary distributions are independent of the friction parameter *γ*_{j}, which is a benefit of schemes based on splitting the T-RPMD dynamics into Hamiltonian and thermostat parts and using the exact Ornstein–Uhlenbeck flow in Eq. (11) to evolve the thermostat part. Moreover, comparing the exact covariance matrix in Eq. (15) with the finite-time step approximations in Eqs. (19)–(23), note that in all cases, **Σ**_{j} = lim_{Δt→0}**Σ**_{j,Δt}. These results have previously been reported for the OBABO [Eqs. (20) and (21), *a* = 1] and BAOAB [Eqs. (20) and (21), *a* = 0] schemes^{8,29} but not for the OBCBO [Eqs. (22) and (23), *a* = 1] or BCOCB [Eqs. (22) and (23), *a* = 0] schemes.

In the infinite-bead limit, the exact and numerical position-marginals can be written as an infinite product of one-dimensional centered normal distributions with variances given by $(\beta mn)\u22121sj2$ and $(\beta mn)\u22121sj,\Delta t2$, respectively. By Kakutani’s theorem,^{30,31} these two distributions have a non-zero overlap if and only if the following series converges:

For OBABO and BAOAB, due to the oscillatory cotangent term appearing in *s*_{j,Δt}, the limit $limj\u2192\u221e(1\u2212sj/sj,\Delta t)2$ does not exist, and therefore, the series does not converge. For OBCBO, the *j*th summand of this series is

which more obviously leads to a divergent series. Therefore, for OBABO, OBCBO, and BAOAB, the numerical stationary distribution has no overlap with the exact stationary distribution in the infinite-bead limit; it is in this sense that these schemes fail to exhibit the property of dimensionality freedom. Remarkably, BCOCB is exact in the position marginal and, thus, exhibits dimensionality freedom. See the Appendix for a brief summary of the properties of other symmetric splittings that were considered.

## IV. CONSEQUENCES FOR THE PRIMITIVE KINETIC ENERGY EXPECTATION VALUE

In the current section, we show that the non-overlap pathology of the OBABO, BAOAB, and OBCBO schemes causes a divergence with increasing bead number of the primitive path-integral kinetic-energy expectation value, an issue that is numerically well known for OBABO and BAOAB.^{8,29,32,33} We further show that this divergence is fully eliminated via the BCOCB scheme—as expected.

The primitive kinetic energy expectation value is given by^{34,35}

where the first equality involves a sum over the ring-polymer beads in Cartesian coordinates (with *q*_{n} = *q*_{0}), and the second equality performs the summation in terms of the ring-polymer normal modes. The divergence of this expectation value is numerically illustrated for the simple case of a harmonic oscillator [Figs. 1(a)–1(d)]; note that for larger MD time steps, the OBABO, BAOAB, and OBCBO schemes fail to reach a plateau with increasing bead number and dramatically deviate from the exact result (dashed line). The same divergence for OBABO and BAOAB has been numerically observed in many systems,^{8,29,32,33} including liquid water which we discuss later. A striking observation from Figs. 1(a)–1(d) is that the BCOCB exhibits no such divergence or error in the primitive kinetic energy expectation value at a high bead number, regardless of the employed time step.

Using Eq. (15), note that the contribution to the primitive kinetic energy expectation value from the *j*th ring-polymer mode is

such that in the infinite-bead limit,

Similarly using Eq. (19), the *j*th-mode contribution to the kinetic energy from the finite-time step numerical expectation value is

Thus, the per-mode error in kinetic energy is

where the per-mode error in the position marginal for internal mode *j* is

where *s*_{j,Δt} is given by Eq. (20) for the cases of OBABO (*a* = 1) and BAOAB (*a* = 0) and by Eq. (22) for the cases of OBCBO (*a* = 1) and BCOCB (*a* = 0). Note that this error vanishes only for the BCOCB scheme, which satisfies *ρ*_{j,Δt} = 0 for each mode *j*, irrespective of the time step Δ*t*.

Equations (29) and (30) indicate that the primitive kinetic energy estimator is a sensitive measure of the finite-time step error in the sampled ring-polymer position distribution associated with the high-frequency modes. Figure 1(e) resolves this per-mode error, *ρ*_{j,Δt}, for each internal mode in simulations that employ a total of 128 beads, including results from OBABO (red), BAOAB (magenta), and OBCBO (blue) using a time step of 1 fs, with the solid lines indicating the analytical predictions in Eq. (30) and with the dots indicating the result of numerical simulations. The analytical results are fully reproduced by the simulations. Note that the OBABO per-mode error exhibits dramatic spikes for *ω*_{j,n}Δ*t* = *kπ*, where 1 ≤ *j* ≤ *n*, and for some *k* ≥ 1, which coincide with the loss of ergodicity of that integration scheme. The BAOAB scheme exhibits these resonance instabilities at even values of *k*. However, it is the failure of this per-mode error to sufficiently decay as a function of the mode number for all three of OBABO, BAOAB, and OBCBO that gives rise upon summation to the divergence of the primitive kinetic energy expectation value, as seen for this particular time step value in Fig. 1(d). Since $\omega j,n2sj2\u21921$ as *n* → *∞*, the convergence of $\u2211j=1\u221e|\u27e8KEj\u27e9\u2212\u27e8KEj\u27e9\Delta t|$ reduces to the convergence of the series $\u2211j=1\u221esj2\u2212sj,\Delta t2$, which diverges for both OBABO and OBCBO due to the same reasons as discussed in Sec. III.

## V. DIMENSIONALITY FREEDOM FOR OBCBO VIA FORCE MOLLIFICATION

Sections III and IV have demonstrated that whereas the BCOCB integrator exhibits dimensionality freedom, the OBCBO integrator does not. In the current section, we show that this shortcoming of OBCBO can be addressed by the use of force mollification, in which the external potential energy in Eq. (3) is replaced by

where $\Omega \u0303$ is any positive semi-definite *n* × *n* matrix that has the same eigenvectors as **Ω** [Eq. (5)] while possibly having different eigenvalues. Force mollification has not previously been employed for PIMD although the strategy originates from a variation-of-constants formulation of the solution to Eq. (2);^{36–39} specifically, the protocol in Eq. (31) is a generalization of the mollified impulse method.^{36}

Use of force mollification in the current work can be motivated on physical grounds: In the absence of a physical potential, four of the considered integration schemes (OBABO, BAOAB, OBCBO, and BCOCB) leave invariant the exact free ring-polymer Boltzmann–Gibbs distribution.^{26} Therefore, the loss of any overlap between the exact stationary distribution of the position marginals in the infinite-bead limit for OBABO, BAOAB, and OBCBO must be attributed to the influence of the time evolution from the external potential in the schemes (i.e., the “B” sub-step) as implemented in Eq. (10); the BCOCB scheme does not suffer from this problem. To remove this pathology in the OBCBO scheme, we thus use mollification to taper down the external forces on the high-frequency modes, such that the resulting integration correctly reverts to free ring-polymer motion for those modes, which should become decoupled from the external potential as the frequency increases. The specific appearance of the 1/2 factor in the sinc function argument ensures that the sinc function switches from its high-frequency effect to its low-frequency effect when the period of the Matsubara frequency is commensurate with Δ*t*; the zero-frequency ring-polymer centroid mode is untouched by mollification.

Force mollification requires only a small algorithmic modification of the OBCBO integrator. Specifically, the “B” sub-step in Eq. (10) is replaced with

where the mollified forces are

where $q\u0303=UD\Delta tUTq$ are the mollified bead positions and where *D*_{Δt} is the diagonal matrix of eigenvalues associated with $sinc(\Omega \u0303\Delta t/2)$, i.e.,

where $\omega \u0303j,n$ is the *j*th eigenvalue of $\Omega \u0303$. In practice, the mollified forces are computed in normal mode coordinates as follows:

Starting with the ring-polymer bead position in normal mode coordinates, obtain a copy of the mollified bead positions via

- (b)
Evaluate the external forces at the mollified ring-polymer bead positions, $F(q\u0303)$.

- (c)
Apply the remaining mollification to the forces in Eq. (33) via

We emphasize that in comparison with the standard force update [Eq. (10)], the use of the mollified force update [Eq. (32)] introduces neither additional evaluations of the external forces nor *n* × *n* matrix multiplies associated with the discrete Fourier transform; it, therefore, avoids any significant additional computational cost.

This mollification scheme preserves reversibility and symplecticity as well as local-third order accuracy of the OBCBO scheme with time step. We emphasize that the sinc-function-based mollification scheme in Eq. (32) is not unique, and alternatives can certainly be devised. Even within the functional form of the mollification in Eq. (32), flexibility remains with regard to the choice of the matrix $\Omega \u0303$, which allows for mode-specificity in the way the mollification is applied. A simple choice for this matrix is $\Omega \u0303=\Omega $, such that mollification is applied to all of the non-zero ring-polymer internal modes. With this choice, we arrive at a fully specified integration scheme that replaces the original “B” sub-step in Eq. (10) with the mollified force sub-step in Eq. (32); we shall refer to this force-mollified version of OBCBO integration scheme as “OMCMO.” In Subsection V A, we propose a partially mollified choice for $\Omega \u0303$ that further improves the accuracy.

For the harmonic external potential, all of the previously derived relations for OBCBO [most notably Eqs. (18), (22), (23), (29), and (30)] also hold for OMCMO with Λ suitably replaced by $\Lambda \u0303j=sinc2(\omega j,n\Delta t/2)\Lambda $. Note that $\Lambda \u0303j\u2264\Lambda $ since sinc^{2}(x) ≤ 1 for all x ≥ 0, making clear that the mollification reduces the effect of the external potential on the higher-frequency internal ring-polymer modes.

We now show that mollifying the forces in the B sub-step fixes the pathologies of OBCBO in the infinite-bead limit, by restoring overlap between the sampled and exact stationary distributions. To see this, note that the *j*th summand in Eq. (24) for OMCMO satisfies

where $f(x)=((1\u2212sinc2(x))/x2+sinc2(x))2$, and we have used the infinite-bead limit for the ring-polymer internal-mode frequencies

Since

we obtain^{40}

Again invoking Kakutani’s theorem [Eq. (24)], it follows that the numerical stationary distribution has an overlap with the exact stationary distribution. As a byproduct of this analysis, we can also quantify the amount of overlap between the exact and numerically sampled stationary distributions,^{41} revealing that the total variation distance^{42} between these distributions is given by

In summary, the force mollification strategy introduced here provably removes the pathologies due to the “B” sub-step in the case of a harmonic oscillator potential. Moreover, for any finite number of beads, the total variation distance between the exact and numerically sampled stationary distribution can be bounded by Eq. (39), and thus, OMCMO admits error bounds that are dimension-free.

Before proceeding, we first return to Fig. 1 to compare the accuracy of OMCMO with the un-mollified OBCBO scheme for the internal-mode position marginal of the harmonic oscillator. As seen in Fig. 1(e) for the results with a time step of 1 fs, the per-mode error obtained by the mollified scheme (OMCMO, green line) decays more rapidly with the mode number than does OBCBO. Figure 1(d) further illustrates that upon summation of the per-mode contributions, the OMCMO prediction for the primitive kinetic energy converges to a well-defined asymptote with respect to the number of ring-polymer beads, whereas OBCBO diverges as discussed earlier. Similar behavior is seen for shorter MD time steps [panels (a)–(c)] although the failure of OBCBO becomes less severe with this range of bead numbers as the time step is reduced.

Although it is satisfying that mollification via OMCMO both formally and numerically ameliorates the problems of the OBCBO scheme in the high-bead-number limit, the OMCMO results in Fig. 1 are not ideal since in some cases, the OMCMO error is substantially larger than that of OBCBO when a modest number of beads is used [e.g., 16 beads in panel (d)]. This observation points to a simple and general refinement of the OMCMO scheme, which we discuss in Subsection V A.

### A. Partial mollification

A comparison of the per-mode errors from OBCBO and OMCMO in Fig. 1(e) reveals that lower errors for OMCMO are only enjoyed for internal modes that exceed a particular frequency (indicated by the vertical black line). This observation suggests that if a “crossover frequency” could be appropriately defined, then a refinement to OMCMO could be introduced for which mollification is applied only to the ring-polymer internal modes with frequency that exceed this crossover value.

For the case of a harmonic external potential, this crossover frequency *ω*_{x} can be found by comparing a bound for the per-mode error [Eq. (30)] for OBCBO,

to that for OMCMO,

where *g*(x) = (1 − sinc^{2}(x))/x^{2} + sinc^{2}(x). Since *g*(x) ≥ 1 only when x ≤ 1, we expect better accuracy if mollification is only applied to those ring-polymer internal modes with frequencies *ω*_{j,n} ≥ *ω*_{x}, where *ω*_{x} = 2/Δ*t*. Although this result was derived for the case of a harmonic potential, it does not depend on Λ. We call this resulting partly mollified integration scheme “OmCmO.” This scheme has the nice properties of OMCMO, including strong stability and dimensionality freedom.

Implementation of OmCmO is a trivial modification of OMCMO, requiring only that the diagonal elements of *D*_{Δt} in Eq. (34) are evaluated using

where *j* = 0, …, *n* − 1. In physical terms, the emergence of 2/Δ*t* in the crossover frequency is intuitive since as was previously mentioned, it corresponds to having the ring-polymer mode undergo a full period per time step Δ*t*.

Finally, numerical results for the case of a harmonic potential [Figs. 1(a)–1(d)] reveal that the partially modified OmCmO scheme (cyan) achieves both robust convergence of the primitive kinetic energy with increasing bead number, as well as comparable or better accuracy than the OBCBO and OMCMO integration schemes–as expected. However, it must be emphasized that for all panels of Fig. 1, the BCOCB scheme (which requires no force mollification) is by far the most accurate and stable.

## VI. RESULTS FOR ANHARMONIC OSCILLATORS

Having numerically characterized the performance of the various non-preconditioned PIMD integrators for the case of the harmonic oscillator external potential in Fig. 1, we now turn our attention to anharmonic external potentials. In this section, we consider both a weakly anharmonic potential,

and the more strongly anharmonic quartic potential,

All calculations are performed using *ℏ* = 1, *m* = 1, and *β* = 1. Assuming the system to be at room temperature (300 K), then the thermal time scale corresponds to *βℏ* ≈ 25.5 fs and Λ = *mω*^{2}, where *ω* = 3315 cm^{−1} for Λ = 256. The trajectories are performed with the centroid mode uncoupled from the thermostat (i.e., in the manner of T-RPMD); for the remaining *n* − 1 internal modes, simulations performed with the OBABO and BAOAB schemes use the standard^{18,20} damping schedule of **Γ** = **Ω**, and simulations performed using the Cayley modification (i.e., BCOCB, OBCBO, OMCMO, and OmCmO) use friction $\gamma j=min(\omega j,n,0.9\gamma jmax(\Lambda ),0.9\gamma jmax(0))$ for the *j*th mode, where $\gamma jmax(\Lambda )$ is the friction that saturates the inequality in Eq. (18); for the quartic potential, we set Λ = 1 in this calculation of $\gamma jmax$.

Figures 2(a) and 2(b) present kinetic energy expectation values for the weakly anharmonic potential corresponding to 3315 cm^{−1} at room temperature. For the primitive kinetic energy expectation value, the results obtained using the various integration schemes with time steps of both 0.5 fs [panel (a)] and 1.0 fs [panel (b)] are consistent with the observations for the harmonic potential in Fig. 1; specifically, the integrators without dimensionality freedom (OBABO, BAOAB, and OBCBO) fail to converge with increasing bead number, while the mollified integrators (OMCMO and OmCmO) smoothly converge with increasing bead number, and the partially mollified scheme (OmCmO) is consistently more accurate than OBCBO and OMCMO. However, it is also clear that BCOCB exhibits the best accuracy with increasing bead number, converging to the exact result without a perceivable time step error.

Figures 2(c) and 2(d) present the corresponding results for the virial kinetic energy expectation value,

where $q\xaf$ is the centroid (bead-averaged) position. Whereas the virial kinetic energy for all of the strongly stable integration schemes is well behaved, the OBABO and BAOAB schemes perform erratically at large time steps due to their provable non-ergodicities.^{26} Appealingly, the BCOCB scheme is consistently the most accurate for the virial kinetic energy expectation value as it was for the primitive kinetic energy expectation value.

Figures 3(a)–3(d) show the results of the various numerical integration schemes for the primitive and virial kinetic energy expectation values as a function of the MD time step using 64 ring-polymer beads. Results are shown for both the weakly anharmonic and quartic oscillators. In all cases, the BCOCB scheme is consistently the most accurate across this array of model systems.

Finally, Fig. 3(e) illustrates the use of the BCOCB integrator for the calculation of real-time quantum dynamics via T-RPMD, replacing the often-employed OBABO integration scheme. Using 64 beads, the T-RPMD results are plotted for a range of integration time steps. Strikingly, over the entire range of considered time steps, BCOCB introduces a negligible error in the calculated position time autocorrelation function; it is confirmed that these results are visually indistinguishable from those obtained using the OBABO integrator in the small-time step limit.

## VII. RESULTS FOR LIQUID WATER

Section VI demonstrated the strong performance of the BCOCB integrator for obtaining both PIMD statistics as well as real-time dynamics via the T-RPMD model, in model systems. In this section, we test the accuracy and stability of the various un-mollified integration schemes (OBABO, OBCBO, BAOAB, and BCOCB) in liquid water, a high-dimensional and relatively complex system. Specifically, we consider a periodic 32-molecule water box at a temperature of 298 K and a density of 0.998 g/cm^{3}, as described by the q-TIP4P/F force field, a flexible, 4-point water model introduced in Ref. 43.

In Fig. 4, we compare the accuracy achieved by the different integrators for the average kinetic energy per hydrogen atom as a function of the number of ring-polymer beads. As in Sec. VI, we consider both the primitive [Eq. (26)] and virial [Eq. (45)] estimators for the kinetic energy. For each choice of integrator, time step, and bead number, the primitive and virial estimators for the kinetic energy per hydrogen atom were averaged over a 1-nanosecond trajectory integrated in the manner of T-RPMD, i.e., with the centroid mode uncoupled from the thermostat; for the remaining *n* − 1 internal modes, simulations performed with the OBABO and BAOAB schemes use the standard^{18,20} damping schedule of **Γ** = **Ω**, and simulations performed using the Cayley modification use friction $\gamma j=min{\omega j,n,0.9\gamma jmax(\omega OH2),0.9\gamma jmax(0)}$, where $\gamma jmax(\Lambda /m)$ saturates the inequality in Eq. (18) for the given values of *j* and Λ/*m* at the given time step, and $\omega OH$ is the OH-stretch frequency from the harmonic bending force field term in the q-TIP4P/F force field.^{43} Multi-nanosecond staging PIMD^{8,11} simulations at a time step of 0.1 fs were performed to obtain a bead-converged reference value for the H-atom kinetic energy, plotted as a dashed line in Figs. 4 and 5.

The primitive kinetic energy expectation values in Figs. 4(a) and 4(b) show similar trends to those seen in Figs. 1 and 2 for the harmonic and weakly anharmonic oscillators. For a 0.5-fs time step [Fig. 4(a)], at which all integrators exhibit strong stability for ring polymers with up to 64 beads at the system temperature,^{26} the OBABO, BAOAB, and OBCBO primitive kinetic energy estimates diverge from the converged result as the number of beads increases, in agreement with the result proven in Sec. III that the error in the ring-polymer configurational distribution generated with these schemes grows unboundedly with increasing bead number. At the larger, 0.8-fs time step [Fig. 4(b)], OBABO and BAOAB formally lose strong stability and their respective primitive kinetic energy estimates dramatically diverge for bead numbers greater than 32; the strongly stable OBCBO scheme also yields a divergent result for the same reason as in Fig. 4(a). As seen on the harmonic and anharmonic model systems, the primitive kinetic energy expectation value from the BCOCB integrator monotonically converges to the reference value with increasing bead number, avoiding any perceptible time step error.

Figures 4(c) and 4(d) show the corresponding virial kinetic energy expectation values. For the smaller time step of 0.5 fs, which is a common choice for path-integral simulations of water, all of the integrators perform similarly. However, upon increasing the time step to 0.8 fs, significant differences in the performance of the integrators emerge, with only BCOCB avoiding a perceptible time step error.

To further compare the accuracy and stability of the OBABO, BAOAB, OBCBO, and BCOCB integrators, Fig. 5 considers the average kinetic energy per hydrogen atom obtained using 64 beads over a wide range of time steps. These results show that BCOCB remains remarkably accurate for time steps as large as 1.4 fs for liquid water, which corresponds to the limit of stability for Verlet integration of the centroid mode. In comparison, OBCBO diverges monotonically as the time step increases, reaching unphysical values for the primitive expectation value and yielding sizable error (20%) for the virial expectation value. The erratic performance of both OBABO and BAOAB is due to the emergence of numerical resonance instabilities at time steps greater than 0.6 fs at the employed bead number; indeed, the largest safe time step at which OBABO and BAOAB remain strongly stable for *n* = 64, Δ*t*_{⋆} ≈ 0.63 fs, precedes the range of time steps in Fig. 5 for which these integrators vary erratically.

Extending beyond statistics, we now consider the dynamical properties of liquid water. Given the superiority of the BCOCB scheme for the calculated statistical properties in Figs. 4 and 5, we present results that focus on this scheme in comparison with the most widely used OBABO scheme. In particular, we consider the liquid water infrared absorption spectrum,^{44} which is proportional to *ω*^{2}*Ĩ*(*ω*) where the dipole spectrum $\u0128(\omega )=\u2009\u2009\u222bRdt\u2009e\u2212i\omega tC\u0303\mu \u22c5\mu (t)$ is the Fourier transform of the Kubo-transformed dipole autocorrelation function $C\u0303\mu \u22c5\mu (t)$. The latter is approximated in the RPMD model by^{16} $C\u0303\mu \u22c5\mu (t)=1N\u2211i=1N\mu \xafi(t)\u22c5\mu \xafi(0)$, where *N* is the number of molecules in the liquid, $\mu \xafi(t)$ is the bead-averaged dipole moment of molecule *i* at time *t*, and the angle brackets denote averaging over the ring-polymer thermal distribution. To obtain the time-correlation functions and spectra shown in Fig. 6 for the OBABO and BCOCB integration schemes, 12-nanosecond T-RPMD trajectories were simulated for a ring polymer with 64 beads and time steps ranging from 0.2 fs to 1.4 fs, using the same friction schedule as described for Figs. 4 and 5.

Along each trajectory, the velocities of all degrees of freedom in the system were drawn anew from the Maxwell–Boltzmann distribution every 20 ps; the autocorrelation function was evaluated out to 2 ps by averaging over staggered windows of that time-length within every 20-ps trajectory segment; and exponential-decay extrapolation was used to extend the autocorrelation function before evaluating its numerical Fourier transform to obtain the infrared absorption spectrum.

Figures 6(a) and 6(b) present the dipole autocorrelation functions obtained using the OBABO and BCOCB integrators with a range of time steps. For the OBABO integrator, the calculated correlation function is qualitatively incorrect for time steps as large as 0.8 fs. For the BCOCB integrator, the resulting correlation functions are far more robust with respect to the time step. Although modest differences are seen in the exponential tail of the correlation function, the dynamics on vibrational time scales (see the inset) is largely unchanged as the time step is varied from 0.2 fs to 1.4 fs. Figure 6(c) further emphasizes this point by showing the absorption spectrum that is obtained from the BCOCB time-correlation functions with the various time steps. To minimize bias, we avoided any smoothing of the spectra shown in panel c. It is clearly seen that the librational and bending features (below 2500 cm^{−1}) are visually indistinguishable over the entire range of considered time steps. To clarify the comparison for the stretching region above 3000 cm^{−1}, we smooth the raw spectra in that region by convolution against a Gaussian kernel with a width of 150 cm^{−1} (see the inset). Again, the robustness of the simulated spectrum over this span of time steps is excellent, with the only significant effect due to finite-time step error being a slight blue-shifting of the OH stretching frequency for the results using a 1.4-fs time step, which is nearly three times larger than the typical value employed for the OBABO scheme for simulations with 64 beads. Taken together, these results indicate that the BCOCB integrator provides an excellent description of both PIMD statistics and T-RPMD dynamics in realistic molecular systems, substantially improving the accuracy and stability of previously employed numerical integrators.

## VIII. SUMMARY

In a previous paper,^{26} we showed that essentially all schemes for the non-preconditioned equations of motion of PIMD, including the widely used OBABO scheme, lack strong stability due to the use of exact free ring-polymer time evolution in the “A” sub-step, and we proved that this lack of strong stability gives rise to a lack of ergodicity in the thermostatted trajectories. We further showed that ergodicity can be restored by simply replacing the “A” sub-step with the Cayley transform.

In the current work, we show that a completely distinct—yet equally important—pathology exists in the “B” sub-step of previously developed non-preconditioned PIMD integrators due to the outsized effect of the external potential on the dynamics of the high-frequency ring-polymer modes. Specifically, we show that previous integrators (including OBABO, BAOAB, and OBCBO) yield a numerical stationary distribution for which the overlap with the exact stationary distribution vanishes in the infinite-bead limit. We then show that this pathology is completely avoided in the BCOCB scheme, and we further show that the pathology can be eliminated for the OBCBO scheme by suitably mollifying the “B” sub-step, yielding the dimension-free non-preconditioned PIMD integrators, namely, BCOCB, OMCMO, and OmCmO. Implementation of the dimension-free integration schemes involves no significant additional computational cost, no additional parameters, and no increase in algorithmic complexity in comparison with either OBABO or BAOAB. Furthermore, since the integrators considered here are all non-preconditioned, they can immediately be used for computing the equilibrium statistical properties as well as dynamical properties via the RPMD model. The numerical performance of the BCOCB scheme is particularly striking, yielding results that are markedly better in terms of accuracy and time step stability than any of the other considered integrators. For liquid water, it is shown that BCOCB allows for time steps as large as 1.4 fs while exhibiting minimal time step error in the calculation of both equilibrium expectation values and the dipole absorption spectrum.

## ACKNOWLEDGMENTS

The authors thank Andreas Eberle and Ondrej Marsalek for helpful discussions. N.B.-R. acknowledges support from the National Science Foundation under Award No. DMS-1816378 and the Alexander von Humboldt Foundation. R.K., J.L.R.-R., 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.

### APPENDIX: OTHER SPLITTINGS

There are exactly four locally third-order accurate symmetric splitting schemes that involve one new force evaluation per integration step and that involve splitting the T-RPMD dynamics into Hamiltonian and thermostat parts: OBABO, BAOAB, OABAO, and ABOBA. In Sec. III, we quantified the properties of OBABO, BAOAB, and their Cayley modifications in the case of a harmonic external potential. The corresponding properties of the Cayley modifications of OABAO and ABOBA are given below.

OCBCO is exact in the velocity marginal, but the variance in the position marginal is $(\beta mn)\u22121sj,\Delta t2$, where $sj,\Delta t2=(4m\u2212\Delta t2\Lambda )/(4\Lambda +4m\omega j,n2)$; and

CBOBC is exact in the position marginal, but the variance in the velocity marginal is $(\beta mn)\u22121rj,\Delta t2$, where $rj,\Delta t2=4m/(4m\u2212\Delta t2\Lambda )$.

Numerical experiments confirmed these properties but did not show significant improvement in accuracy compared with BCOCB. Therefore, we did not include numerical results for these schemes.

## REFERENCES

In the special case when Λ = 0, the given condition for OBCBO corrects a sign error in Eq. 37 of Ref. 26.

This inequality comes from using Eq. (37) to write $\u2211j=1\u221ef(\omega j\Delta t/2)=I+II$, where $I=2\u2211j=1\u230a\u210f\beta /(\pi \Delta t)\u230bf(j\pi \Delta t/(\u210f\beta ))$ and $II=2\u2211j=\u2308\u210f\beta /(\pi \Delta t)\u2309\u221ef(j\pi \Delta t/(\u210f\beta ))$. Then, the first term admits the bound I ≤ 2*f*(1)*ℏβ*/(*π*Δ*t*) < 4*ℏ*β/(*π*Δ*t*), and for the second term, we use $II\u2264F(1)+\u210f\beta /(\pi \Delta t)\u222b1\u221eF(x)dx$, where $F(x)=2((1\u2212sinc2(x))/x2+1/x2)2$ is monotone decreasing on [1, ∞) with *F*(1) ≤ 4 and $\u222b1\u221eF(x)dx\u22642$.

This quantification uses (i) *d*_{TV} ≤ 2 *d*_{H}, where *d*_{TV} is the total variation distance and *d*_{H} is the Hellinger distance; and (ii) subadditivity of the squared Hellinger distance, which implies that $dH2(\mu ,\mu \Delta t)\u2264\u2211j=1\u221edH2(N(0,sj2),N(0,sj,\Delta t2))\u2264\u2211j=1\u221e(1\u2212sj2/sj,\Delta t2)2\u22643\u210f\beta /(\pi \Delta t)+2\Delta t4\Lambda 28\u2009m2$.