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.

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 efficient4–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 RPMD12,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 form11,16–25

eΔtLeaΔt2OeΔt2BeΔt2Ae(1a)ΔtOeΔt2AeΔt2BeaΔt2O,
(1)

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 Parrinello22 when a = 1 and the “BAOAB” scheme of Leimkuhler25 when a = 0.

In our previous work,26 we emphasized that earlier PIMD numerical integration schemes had overlooked a fundamental aspect of the exp((Δ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((Δ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.

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

q̇(t)=v(t),v̇(t)=Ω2q(t)+1mnF(q(t))Γv(t)+2βmnΓ1/2Ẇ(t).
(2)

Here, W is an n-dimensional standard Brownian motion; q(t) = (q0(t), …, qn−1(t)) is the vector of positions for the n ring-polymer beads at time t ≥ 0, and v(t) are the corresponding velocities; mn = m/n and β=(kBT)1; and F(q)=Vnext(q), where Vnext is the contribution of the external potential,

Vnext(q)=1nj=0n1V(qj).
(3)

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

Ω2=κn221001121000012110012,
(4)

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

Ω=U  diag(0,ω1,n,,ωn1,n)UT,
(5)

where ωj,n is the jth Matsubara frequency27 given by

ωj,n=2κnsinπj2n,ifjis even2κnsinπ(j+1)2n,else.
(6)

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

Γ=U  diag(0,γ1,,γn1)UT,
(7)

where γj is the friction factor in the jth 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(−βHn(q, v)), where Hn(q, v) is the ring-polymer Hamiltonian defined by

Hn(q,v)=Hn0(q,v)+Vnext(q),
(8)

and Hn0(q,v)=(1/2)mn|v|2+qTΩ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 τ,

qvexp(τA)qv,
(9)

where A=0IΩ20 is the Hamiltonian matrix associated with the free ring polymer, (ii) velocity updates of time step τ due to forces from the external potential,

vv+τ1mnF(q),
(10)

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

vexp(τΓ)v+1βmn(Iexp(2τΓ))1/2ξ,
(11)

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 = /ωj,n for any 1 ≤ jn and k ≥ 1. We also identified a maximum safe time step size Δt = βℏπ/(2n), 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 modification26 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

cay(ΔtA)=(I(1/2)ΔtA)1(I+(1/2)ΔtA).
(12)

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(ΔtA)1/2. While it might be expected that these half-time step updates would, instead, be replaced with cay((Δt/2)A), 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(ΔtA)1/2cay(ΔtA)1/2 = cay(ΔtA). 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 (V ≡ const.).26 

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 jth internal ring-polymer mode with frequency ωj,n, in the presence of a harmonic external potential V(q) = (1/2)Λq2 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

ϱ=UTq  and  φ=UTv,
(13)

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

ϱ̇j(t)φ̇j(t)=Kjϱj(t)φj(t)+02β1mn1γjẆj(t)Kj=Aj+B+Oj,
(14)

where Ẇj is a scalar white-noise, and we have introduced the following 2 × 2 matrices:

Aj=01ωj,n20,B=00Λ/m0, and Oj=000γj.

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,

Σj=1βmnsj2001,  sj2=1Λ/m+ωj,n2.
(15)

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

ϱj(t+Δt)φj(t+Δt)=Mjϱj(t)φj(t)+Rj1/2ξ0η0,
(16)

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

Mj=eaΔt2OjeΔt2BeΔt2Aje(1a)ΔtOjeΔt2AjeΔt2BeaΔt2Oj,Rj=1e2(1a)γjΔtβmnNjPNjT+1eaγjΔtβmn(MjeaΔt2Oj)P(MjeaΔt2Oj)T+P,

where P=0001 and Nj=eaΔt2OjeΔt2BeΔt2Aj. The corresponding step for the Cayley modification is obtained by replacing exp((Δt/2)Aj) in Eq. (16) with cay(ΔtAj)1/2, which is given by

cay(ΔtAj)1/2=14+ωj,n2Δt22Δtωj,n2Δt2.
(17)

A sufficient condition28 for ergodicity of Eq. (16) is

   1>Aj,Δt2cosh2((Δt/2)γj),
(18)

where

Aj,Δt=cos(Δtωj,n)(Λ/m)Δt2ωj,nsin(Δtωj,n).

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

Aj,Δt=1+82(Λ/m)Δt24+ωj,n2Δt2.

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 = /ω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 Σjt that satisfies the linear equation,

   Σj,Δt=MjΣj,ΔtMjT+Rj,

for which the solution is

Σj,Δt=1βmnsj,Δt200rj,Δt2,
(19)

where the variance in the position and velocity marginal are (βmn)1sj,Δt2 and (βmn)1rj,Δt2 with

sj,Δt2=1ωj,n2+ΛΔtωj,nmcot(Δtωj,n)(ΛΔt2m)2,a=11ωj,n2+ΛΔtωj,n2mcot(Δt2ωj,n),a=0,
(20)
rj,Δt2=1,a=12mωj,nΛΔttan(Δt2ωj,n)2mωj,n,a=0.
(21)

For the Cayley modification of Eq. (16),

sj,Δt2=4m4maΔt2Λsj2,
(22)
rj,Δt2=4m(1a)Δt2Λ4m.
(23)

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Σjt. These results have previously been reported for the OBABO [Eqs. (20) and (21), a = 1] and BAOAB [Eqs. (20) and (21), a = 0] schemes8,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 (βmn)1sj2 and (βmn)1sj,Δt2, respectively. By Kakutani’s theorem,30,31 these two distributions have a non-zero overlap if and only if the following series converges:

j=11sjsj,Δt2.
(24)

For OBABO and BAOAB, due to the oscillatory cotangent term appearing in sjt, the limit limj(1sj/sj,Δt)2 does not exist, and therefore, the series does not converge. For OBCBO, the jth summand of this series is

Δt4Λ216m21+4mΔt2Λ4m2,

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.

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 by34,35

KEprim=n2βj=1nmnκn22(qjqj1)2,
(25)
=12β+j=1n112βmnωj,n22ϱj2,
(26)

where the first equality involves a sum over the ring-polymer beads in Cartesian coordinates (with qn = q0), 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.

FIG. 1.

Primitive kinetic energy expectation values for a harmonic potential V(q)=12Λq2, with Λ = 256, = m = 1, and reciprocal temperature β = 1; choosing energies to be in units of kBT at room temperature (300 K), then βℏ ≈ 25.5 fs and Λ = 2, where ω = 3315 cm−1. (a)–(d) For various MD time steps, the primitive kinetic energy expectation value as a function of the number of ring-polymer beads, with the exact kinetic energy indicated as a dashed gray line. The standard error of all visible data points in each plot is smaller than the symbol size. (e) Per-mode error in the variance of position coordinate of the normal modes for simulations run with 128 ring-polymer beads and a time step of 1 fs; solid lines are analytic predictions from Eq. (30) with (20) and (22) defining sj,Δt2 for the different schemes; points indicate the results of numerical PIMD simulations using the various integration schemes. The BCOCB scheme is not shown since it has zero error for all internal modes. The black vertical line indicates the crossover frequency (ωx = 2/Δt) for the error of OBCBO and OMCMO based on the bounds in Eqs. (40) and (41).

FIG. 1.

Primitive kinetic energy expectation values for a harmonic potential V(q)=12Λq2, with Λ = 256, = m = 1, and reciprocal temperature β = 1; choosing energies to be in units of kBT at room temperature (300 K), then βℏ ≈ 25.5 fs and Λ = 2, where ω = 3315 cm−1. (a)–(d) For various MD time steps, the primitive kinetic energy expectation value as a function of the number of ring-polymer beads, with the exact kinetic energy indicated as a dashed gray line. The standard error of all visible data points in each plot is smaller than the symbol size. (e) Per-mode error in the variance of position coordinate of the normal modes for simulations run with 128 ring-polymer beads and a time step of 1 fs; solid lines are analytic predictions from Eq. (30) with (20) and (22) defining sj,Δt2 for the different schemes; points indicate the results of numerical PIMD simulations using the various integration schemes. The BCOCB scheme is not shown since it has zero error for all internal modes. The black vertical line indicates the crossover frequency (ωx = 2/Δt) for the error of OBCBO and OMCMO based on the bounds in Eqs. (40) and (41).

Close modal

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

KEj=12β1ωj,n2sj2,

such that in the infinite-bead limit,

limnj=0n1KEj=4Λm1+2eβΛ/m1.
(27)

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

KEjΔt=12β1ωj,n2sj,Δt2.
(28)

Thus, the per-mode error in kinetic energy is

|KEjKEjΔt|=mnωj,n22ρj,Δt,
(29)

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

ρj,Δt=1βmnsj2sj,Δt2,
(30)

where sjt 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 ρjt = 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, ρjt, 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 = , where 1 ≤ jn, 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 ωj,n2sj21 as n, the convergence of j=1|KEjKEjΔt| reduces to the convergence of the series j=1sj2sj,Δt2, which diverges for both OBABO and OBCBO due to the same reasons as discussed in Sec. III.

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

Ṽnext(q)=Vnext(sinc(Ω̃Δt/2)q),
(31)

where Ω̃ 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

vv+Δt21mnF̃(q),
(32)

where the mollified forces are

F̃(q)=sinc(Ω̃Δt/2)F(q̃)=UDΔtUTF(q̃),
(33)

where q̃=UDΔtUTq are the mollified bead positions and where DΔt is the diagonal matrix of eigenvalues associated with sinc(Ω̃Δt/2), i.e.,

DΔt=  diag(sinc(ω̃0,nΔt/2),,sinc(ω̃n1,nΔt/2)),
(34)

where ω̃j,n is the jth eigenvalue of Ω̃. 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

q̃=UDΔtϱ.
(35)
  • (b)

    Evaluate the external forces at the mollified ring-polymer bead positions, F(q̃).

  • (c)

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

UTF̃(q)=DΔtUTF(q̃).
(36)

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 Ω̃, which allows for mode-specificity in the way the mollification is applied. A simple choice for this matrix is Ω̃=Ω, 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 Ω̃ 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 Λ̃j=sinc2(ωj,nΔt/2)Λ. Note that Λ̃jΛ since sinc2(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 jth summand in Eq. (24) for OMCMO satisfies

1sjsj,Δt21sj2sj,Δt22f(ωjΔt/2)Δt4Λ216m2,

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

ωj=limnωj,n=πjβ,ifjis evenπ(j+1)β,else.
(37)

Since

j=1f(ωjΔt/2)6βπΔt+4,  

we obtain40 

j=11sjsj,Δt26βπΔt+4Δt4Λ216m2.
(38)

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 distance42 between these distributions is given by

dTV(μ,μΔt)6βπΔt+4Δt2Λ2m.
(39)

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 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,

ρj,Δtsj2mnωj,n2βΔt2Λ4mΔt2Λ,
(40)

to that for OMCMO,

ρj,Δtg(ωj,nΔt/2)sj2mnωj,n2βΔt2Λ4mΔt2Λ,
(41)

where g(x) = (1 − sinc2(x))/x2 + sinc2(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

sinc(ω̃j,nΔt/2)=1,for ωj,n<ωxsinc(ωj,nΔt/2),otherwise,
(42)

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.

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,

V(q)=Λ12q2+110q3+1100q4,
(43)

and the more strongly anharmonic quartic potential,

V(q)=14q4.
(44)

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 Λ = 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 standard18,20 damping schedule of Γ = Ω, and simulations performed using the Cayley modification (i.e., BCOCB, OBCBO, OMCMO, and OmCmO) use friction γj=min(ωj,n,0.9γjmax(Λ),0.9γjmax(0)) for the jth mode, where γjmax(Λ) is the friction that saturates the inequality in Eq. (18); for the quartic potential, we set Λ = 1 in this calculation of γ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.

FIG. 2.

Primitive and virial kinetic energy expectation values as a function of bead number for the weakly anharmonic potential corresponding to 3315 cm−1 at room temperature, with results obtained using a time step of 0.5 fs [(a) and (c)] and 1.0 fs [(b) and (d)]. The standard error of all visible data points in each plot is smaller than the symbol size. The exact kinetic energy is indicated with a dashed line.

FIG. 2.

Primitive and virial kinetic energy expectation values as a function of bead number for the weakly anharmonic potential corresponding to 3315 cm−1 at room temperature, with results obtained using a time step of 0.5 fs [(a) and (c)] and 1.0 fs [(b) and (d)]. The standard error of all visible data points in each plot is smaller than the symbol size. The exact kinetic energy is indicated with a dashed line.

Close modal

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

KEvirial=12β12(qq¯)F(q),
(45)

where q¯ 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.

FIG. 3.

Primitive and virial kinetic energy expectation values as a function of the time step for the weakly anharmonic potential corresponding to 3315 cm−1 at room temperature [(a) and (b)] and the quartic potential [(c) and (d)]. The exact kinetic energy is indicated with a dashed line. The standard error of all visible data points in each plot is smaller than the symbol size. The position autocorrelation function (e) for the quartic oscillator at room temperature computed using T-RPMD with the BCOCB integrator. Results are obtained using 64 ring-polymer beads using time steps of Δt = 0.125 fs, 2 fs, 4 fs, and 8 fs.

FIG. 3.

Primitive and virial kinetic energy expectation values as a function of the time step for the weakly anharmonic potential corresponding to 3315 cm−1 at room temperature [(a) and (b)] and the quartic potential [(c) and (d)]. The exact kinetic energy is indicated with a dashed line. The standard error of all visible data points in each plot is smaller than the symbol size. The position autocorrelation function (e) for the quartic oscillator at room temperature computed using T-RPMD with the BCOCB integrator. Results are obtained using 64 ring-polymer beads using time steps of Δt = 0.125 fs, 2 fs, 4 fs, and 8 fs.

Close modal

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.

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/cm3, 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 standard18,20 damping schedule of Γ = Ω, and simulations performed using the Cayley modification use friction γj=min{ωj,n,0.9γjmax(ωOH2),0.9γjmax(0)}, where γjmax(Λ/m) saturates the inequality in Eq. (18) for the given values of j and Λ/m at the given time step, and ωOH is the OH-stretch frequency from the harmonic bending force field term in the q-TIP4P/F force field.43 Multi-nanosecond staging PIMD8,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.

FIG. 4.

Primitive and virial kinetic energy expectation values as a function of the bead number per hydrogen atom in liquid water at 298 K and 0.998 g/cm3 at time step Δt = 0.5 fs [(a) and (c)] and Δt = 0.8 fs [(b) and (d)]. The reference kinetic energy, obtained from a converged staging PIMD simulation at time step Δt = 0.1 fs and bead number n = 256, is indicated with a dashed line. The standard error of all visible data points in each plot is smaller than the symbol size.

FIG. 4.

Primitive and virial kinetic energy expectation values as a function of the bead number per hydrogen atom in liquid water at 298 K and 0.998 g/cm3 at time step Δt = 0.5 fs [(a) and (c)] and Δt = 0.8 fs [(b) and (d)]. The reference kinetic energy, obtained from a converged staging PIMD simulation at time step Δt = 0.1 fs and bead number n = 256, is indicated with a dashed line. The standard error of all visible data points in each plot is smaller than the symbol size.

Close modal
FIG. 5.

Primitive and virial kinetic energy expectation values as a function of the time step per hydrogen atom in liquid water at 298 K and 0.998 g/cm3, as described by a 64-bead ring polymer. The reference kinetic energy, obtained from a converged staging PIMD simulation at time step Δt = 0.1 fs and bead number n = 256, is indicated with a dashed line. The standard error of all visible data points in each plot is smaller than the symbol size.

FIG. 5.

Primitive and virial kinetic energy expectation values as a function of the time step per hydrogen atom in liquid water at 298 K and 0.998 g/cm3, as described by a 64-bead ring polymer. The reference kinetic energy, obtained from a converged staging PIMD simulation at time step Δt = 0.1 fs and bead number n = 256, is indicated with a dashed line. The standard error of all visible data points in each plot is smaller than the symbol size.

Close modal

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 Ĩ(ω)=  RdteiωtC̃μμ(t) is the Fourier transform of the Kubo-transformed dipole autocorrelation function C̃μμ(t). The latter is approximated in the RPMD model by16C̃μμ(t)=1Ni=1Nμ¯i(t)μ¯i(0), where N is the number of molecules in the liquid, μ¯i(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.

FIG. 6.

Dynamical properties of liquid water computed using T-RPMD with the (a) OBABO and [(b) and (c)] BCOCB integration schemes. Panels (a) and (b) present the Kubo-transformed dipole autocorrelation function computed with various time steps, and panel (c) presents the absorption spectrum from the BCOCB correlation function at each time step. The inset to panel (c) presents the OH stretching region with smoothing.

FIG. 6.

Dynamical properties of liquid water computed using T-RPMD with the (a) OBABO and [(b) and (c)] BCOCB integration schemes. Panels (a) and (b) present the Kubo-transformed dipole autocorrelation function computed with various time steps, and panel (c) presents the absorption spectrum from the BCOCB correlation function at each time step. The inset to panel (c) presents the OH stretching region with smoothing.

Close modal

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.

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.

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.

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 (βmn)1sj,Δt2, where sj,Δt2=(4mΔt2Λ)/(4Λ+4mωj,n2); and

  • CBOBC is exact in the position marginal, but the variance in the velocity marginal is (βmn)1rj,Δt2, where rj,Δt2=4m/(4mΔt2Λ).

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.

1.
M.
Parrinello
and
A.
Rahman
,
J. Chem. Phys.
80
,
860
(
1984
).
2.
G. J.
Martyna
,
A.
Hughes
, and
M. E.
Tuckerman
,
J. Chem. Phys.
110
,
3275
(
1999
).
3.
P.
Minary
,
G. J.
Martyna
, and
M. E.
Tuckerman
,
J. Chem. Phys.
118
,
2510
(
2003
).
4.
A.
Beskos
,
G.
Roberts
,
A.
Stuart
, and
J.
Voss
,
Stochast. Dynam.
8
,
319
(
2008
).
5.
A.
Beskos
,
F.
Pinski
,
J.
Sanz-Serna
, and
A.
Stuart
,
Stochastic Process. Appl.
121
,
2201
(
2011
).
6.
J.
Lu
,
Y.
Lu
, and
Z.
Zhou
, arXiv:1811.10995 (
2020
).
7.
Z.
Zhang
,
X.
Liu
,
Z.
Chen
,
H.
Zheng
,
K.
Yan
, and
J.
Liu
,
J. Chem. Phys.
147
,
034109
(
2017
).
8.
J.
Liu
,
D.
Li
, and
X.
Liu
,
J. Chem. Phys.
145
,
024103
(
2016
).
9.
N.
Bou-Rabee
and
J. M.
Sanz-Serna
,
Acta Numer.
27
,
113
(
2018
).
10.
N.
Bou-Rabee
and
A.
Eberle
,
Stoch. PDE (in press)
; arXiv:1909.07962 (
2019
).
11.
M. E.
Tuckerman
,
B. J.
Berne
,
G. J.
Martyna
, and
M. L.
Klein
,
J. Chem. Phys.
99
,
2796
(
1993
).
12.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
121
,
3368
(
2004
).
13.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
 III
,
Annu. Rev. Phys. Chem.
64
,
387
(
2013
).
14.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill
,
1965
).
15.
D.
Chandler
and
P. G.
Wolynes
,
J. Chem. Phys.
74
,
4078
(
1981
).
16.
T. F.
Miller
 III
and
D. E.
Manolopoulos
,
J. Chem. Phys.
123
,
154504
(
2005
).
17.
T. F.
Miller
 III
and
D. E.
Manolopoulos
,
J. Chem. Phys.
122
,
184503
(
2005
).
18.
M.
Ceriotti
,
M.
Parrinello
,
T. E.
Markland
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
133
,
124104
(
2010
).
19.
M.
Ceriotti
,
D. E.
Manolopoulos
, and
M.
Parrinello
,
J. Chem. Phys.
134
,
084104
(
2011
).
20.
M.
Rossi
,
M.
Ceriotti
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
140
,
234116
(
2014
).
21.
M.
Rossi
,
V.
Kapil
, and
M.
Ceriotti
,
J. Chem. Phys.
148
,
102301
(
2018
).
22.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
,
J. Chem. Phys.
126
,
014101
(
2007
).
23.
N.
Bou-Rabee
and
H.
Owhadi
,
SIAM J. Numer. Anal.
48
,
278
(
2010
).
24.
N.
Bou-Rabee
and
E.
Vanden-Eijnden
,
Commun. Pure Appl. Math.
63
,
655
(
2010
).
25.
B.
Leimkuhler
and
C.
Matthews
,
Appl. Math. Res. Express
2013
,
34
.
26.
R.
Korol
,
N.
Bou-Rabee
, and
T. F.
Miller
 III
,
J. Chem. Phys.
151
,
124103
(
2019
).
27.
T.
Matsubara
,
Prog. Theor. Phys.
14
,
351
(
1955
).
28.

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

29.
D.
Li
,
X.
Han
,
Y.
Chai
,
C.
Wang
,
Z.
Zhang
,
Z.
Chen
,
J.
Liu
, and
J.
Shao
,
J. Chem. Phys.
147
,
184104
(
2017
).
30.
S.
Kakutani
,
Ann. Math.
49
,
214
(
1948
).
31.
V. I.
Bogachev
,
Gaussian Measures
(
American Mathematical Society
,
1998
), Vol. 62.
32.
A.
Pérez
and
M. E.
Tuckerman
,
J. Chem. Phys.
135
,
064104
(
2011
).
33.
O.
Marsalek
,
P.-Y.
Chen
,
R.
Dupuis
,
M.
Benoit
,
M.
Méheut
,
Z.
Bačić
, and
M. E.
Tuckerman
,
J. Chem. Theory Comput.
10
,
1440
(
2014
).
34.
J. A.
Barker
,
J. Chem. Phys.
70
,
2914
(
1979
).
35.
M. F.
Herman
,
E. J.
Bruskin
, and
B. J.
Berne
,
J. Chem. Phys.
76
,
5150
(
1982
).
36.
B.
García-Archilla
,
J. M.
Sanz-Serna
, and
R. D.
Skeel
,
SIAM J. Sci. Comput.
20
,
930
(
1998
).
37.
E.
Hairer
and
C.
Lubich
,
SIAM J. Numer. Anal.
38
,
414
(
2000
).
38.
J. M.
Sanz-Serna
,
SIAM J. Sci. Comput.
46
,
1040
(
2008
).
39.
R. I.
McLachlan
and
A.
Stern
,
SIAM J. Sci. Comput.
52
,
1378
(
2014
).
40.

This inequality comes from using Eq. (37) to write j=1f(ωjΔt/2)=I+II, where I=2j=1β/(πΔt)f(jπΔt/(β)) and II=2j=β/(πΔt)f(jπΔt/(β)). Then, the first term admits the bound I ≤ 2f(1)ℏβ/(πΔt) < 4β/(πΔt), and for the second term, we use IIF(1)+β/(πΔt)1F(x)dx, where F(x)=2((1sinc2(x))/x2+1/x2)2 is monotone decreasing on [1, ∞) with F(1) ≤ 4 and 1F(x)dx2.

41.

This quantification uses (i) dTV ≤ 2 dH, where dTV is the total variation distance and dH is the Hellinger distance; and (ii) subadditivity of the squared Hellinger distance, which implies that dH2(μ,μΔt)j=1dH2(N(0,sj2),N(0,sj,Δt2))j=1(1sj2/sj,Δt2)23β/(πΔt)+2Δt4Λ28m2.

42.
A. L.
Gibbs
and
F. E.
Su
,
Int. Stat. Rev.
70
,
419
(
2002
); arXiv:math/0209021.
43.
S.
Habershon
,
T. E.
Markland
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
131
,
024501
(
2009
).
44.
S.
Habershon
,
G. S.
Fanourgakis
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
129
,
074501
(
2008
).