We introduce a scheme for approximating quantum time correlation functions numerically within the Feynman path integral formulation. Starting with the symmetrized version of the correlation function expressed as a discretized path integral, we introduce a change of integration variables often used in the derivation of trajectory-based semiclassical methods. In particular, we transform to sum and difference variables between forward and backward complex-time propagation paths. Once the transformation is performed, the potential energy is expanded in powers of the difference variables, which allows us to perform the integrals over these variables analytically. The manner in which this procedure is carried out results in an open-chain path integral (in the remaining sum variables) with a modified potential that is evaluated using imaginary-time path-integral sampling rather than requiring the generation of a large ensemble of trajectories. Consequently, any number of path integral sampling schemes can be employed to compute the remaining path integral, including Monte Carlo, path-integral molecular dynamics, or enhanced path-integral molecular dynamics. We believe that this approach constitutes a different perspective in semiclassical-type approximations to quantum time correlation functions. Importantly, we argue that our approximation can be systematically improved within a cumulant expansion formalism. We test this approximation on a set of one-dimensional problems that are commonly used to benchmark approximate quantum dynamical schemes. We show that the method is at least as accurate as the popular ring-polymer molecular dynamics technique and linearized semiclassical initial value representation for correlation functions of linear operators in most of these examples and improves the accuracy of correlation functions of nonlinear operators.

The numerical computation of equilibrium quantum time correlation functions remains one of the outstanding challenges in the theoretical treatment of nuclear quantum effects. This challenge constitutes a serious stumbling block to the accurate determination of both transport properties, which are determined from time correlation functions via the Green-Kubo theory, and vibrational spectra, which are fundamentally linked to integral transforms of time correlation functions. In principle, if the eigenvalues and eigenvectors of a quantum system could be determined, then the calculation of time correlation functions, which requires the time propagation of one of the operators in the correlation function, would be rendered straightforward. However, for condensed-phase systems, this is generally intractable. As an alternative to the direct calculation of eigenstates, the Feynman path integral formulations of quantum mechanics1–4 and quantum statistical mechanics are most often the approach of choice, and for static quantum properties, path integral techniques are both exact within the choice of the Hamiltonian and computationally efficient. Because static properties are derived from the trace of the Boltzmann operator exp(−βĤ), where Ĥ is the Hamiltonian and β = 1/kBT, the path integral formulation is tantamount to sampling the configuration space of a set of ring-polymers composed of pseudoparticles or “beads” indexed by the discretization of an imaginary time parameter τ = βℏ determined by the temperature.5 Unfortunately, when it comes to quantum time correlation functions, the path integral formulation does not yield a computationally tractable problem, as the time propagators needed to evolve an operator in the Heisenberg picture, exp (±iĤt/), give rise to highly oscillatory phase factors in the path integral that are extremely difficult to converge.

Among the various approximations that are made in the calculation of quantum time correlation functions, the simplest method is to replace the quantum correlation function with its classical counterpart, which clearly constitutes a complete neglect of quantum effects. For vibrational spectra, quantum corrections can be incorporated via frequency-dependent prefactors based on an exact relation between classical and quantum time correlation functions for a harmonic oscillator. Not surprisingly, this approximation is satisfactory if the correlation function in question does not deviate strongly from either the harmonic or classical limits. Another class of highly popular techniques that exploit the connection between quantum and classical correlation functions are the centroid molecular dynamics6–9 (CMD) and ring-polymer molecular dynamics10 (RPMD) methods. These methods exploit the path integral molecular dynamics (PIMD) methodology,2,5,11–13 which is a widely used approach to evaluate imaginary time path integrals for static properties, ascribing a quantum dynamical character to the molecular dynamics evolution of the ring-polymer centroid or ring-polymer bead positions and using this dynamics to compute classical-like dynamical correlations of functions of these centroid or bead positions. These produce approximations to the Kubo-transformed quantum time correlation function. The primary difference between CMD and RPMD is the choice of fictitious dynamical masses of the beads and the manner in which thermostats are employed to produce Boltzmann sampling. However, because these are essentially imaginary time methods, they do not explicitly include any phase information. Consequently, systems with strong quantum coherence are treated rather poorly within the CMD and RPMD approximations. Furthermore, while various methods have been proposed,7,14–18 CMD and RPMD methods continue to exhibit difficulties in evaluating nonlinear correlation functions. Other approaches, including the thermal Gaussian molecular dynamics (TGMD) methods,19,20 which are conceptually similar to CMD, provide an alternative mapping of a quantum system to a classical system using an effective potential derived from the variational Gaussian wave packet approximation.21 While computationally efficient, TGMD provides slightly improved results when compared to RPMD and CMD for some correlation functions.

Beyond CMD and RPMD, it is also possible to employ semiclassical approaches that do treat phase factors albeit approximately via linearization techniques.22–36 These methods exist in a hierarchy, with some requiring generation of classical trajectories in very large numbers in order to converge the semi-classical expressions with an associated high computational overhead but possessing the ability to handle quantum coherence,22 and others requiring fewer trajectories but being more likely to break down when strong quantum coherence is present. The basic idea of most semiclassical trajectory-based approaches is the exploitation of a transformation to sum and difference variables between forward and backward time propagation paths, and the linearization is performed with respect to the difference variables. This type of transformation can be applied to the various formulations of the time correlation function including the standard24–26,28–30,33–35 and symmetrized versions31–33 and can even be carried out to second order31 or be used to create trajectory-based Monte Carlo algorithms.32 Finally, in certain cases, it has been shown that the phase factors can be computed via quadrature techniques using points sampled from the corresponding positive-definite quantum Boltzmann distribution.37–40 The putative scaling of such an approach, however, would be expected to grow exponentially with system size and is, therefore, difficult to apply to more than a small number of quantum degrees of freedom, although Makri and co-workers37–40 have shown that optimal Monte Carlo sampling of the grid points improves the scaling considerably.

In this paper, we consider a somewhat different perspective on the calculation of quantum time correlation functions, which, though inspired by the aforementioned semiclassical ideas, nevertheless, aims at a purely imaginary-time sampling scheme capable of capturing approximate phase information explicitly. The starting point is the symmetrized time correlation function alluded to above, which uses complex-time forward and backward propagators to evolve the quantum operators and which can be written as a path integral that possesses a positive-definite density and an oscillatory phase factor, both of which depend on the same path integral bead variables.41,42 The bead variables in this formulation of the symmetrized time correlation function are then transformed to sum and difference variables along the discretized complex-time paths, which, as previously noted, is not new.24–26,28–35 However, the manner and context in which it is applied here are. Once the transformation is made, the potential is then expanded as a power series in the difference variables, and the resulting integrals over these variables are then performed analytically. Following this procedure in the present context gives rise to a purely imaginary time path integral in the remaining sum variables that are subject to a modified potential. In this way, we calculate the approximate quantum time correlation function as a standard path-integral sampling problem rather than a costly trajectory-sampling problem. Nevertheless, phase information is included up to second order in the path difference variables. The resulting imaginary time path integral can be evaluated by any numerical scheme commonly used for path integrals, including Monte Carlo,43 molecular dynamics,13 hybrid Monte Carlo,13,44 or any path-integral enhanced sampling scheme.45,46 Moreover, the approximation can be systematically improved via a cumulant expansion carried out to any desired order. Within such a cumulant expansion, the integrals over the path difference variables can be performed analytically at each order, and the remaining path integral can be shown to be real and positive definite at any order and can, in principle, be evaluated using the same imaginary-time path integral sampling schemes mentioned above. Even at second order, the resulting correlation functions are found, in various cases, to be of comparable accuracy to those produced by the methods discussed above for linear operators and are generally of improved accuracy for nonlinear operators. Thus, it is possible to compute approximate quantum time correlation functions in which both the approximation made and a path to systematic improvement of the approximation are manifestly clear.

This paper is organized as follows: In Sec. II, we lay out the theory behind our new approach. Following this, a number of simple numerical examples are presented. Computational details of these examples are presented in Sec. III, and the results are shown in Sec. IV. Conclusions, further discussion, and prospects for future developments are given in Sec. V. Finally, in  Appendixes A and  B, we derive the first cumulant correction at third order in the path difference variables and the Suzuki-Chin47,48 formulation of our scheme.

The theoretical development begins with the symmetrized quantum time correlation function GAB(t), which is defined as

GAB(t)=1ZTrÂeiĤτc*/B^eiĤτc/.
(1)

In this representation, Ĥ is the Hamiltonian, Z is the canonical partition function, Â and B^ are the operators of interest, and τc is a complex time variable such that τc = tiβℏ/2. Thus, in the symmetrized time correlation function framework, the operator B^ is evolved in the complex time variable τc that mixes temperature and time. The symmetrized correlation function GAB(t) and the true time correlation function, CAB(t), are related through their Fourier transforms via the relation

C̃AB(ω)=eβω/2G̃AB(ω).
(2)

In the analysis to be presented here, we will consider a single particle in one dimension with Hamiltonian Ĥ=p^2/2 m+U(x^). Assuming  and B^ are functions only of the position operator x^, we proceed by evaluating the trace in the coordinate basis, which gives

GAB(t)=1Zdxdxx|eiHτc*/|xb(x)x|eiHτc/|xa(x).
(3)

Unlike CAB(t), GAB(t) is a real object since the two matrix elements are complex conjugates of each other, and a(x) and b(x′), which are the eigenvalues of  and B^ at x and x, respectively, are real. The integrand in Eq. (3) can be interpreted as follows: start at the point x and measure the operator  at that point; the point x is then propagated forward in complex time τc to x′; at x′, the operator B^ is measured; finally, the point x′ is propagated backward in complex time back to x. Using the Trotter factorization method, Krilov and Berne41 demonstrated that GAB(t) can be written in a path integral formulation by expressing each matrix element as a discrete path integral with P beads. Thus, the path-integral representation of GAB(t) contains 2P beads and takes the form

GAB,P(t)=1Zdx1dx2Pa(x1)b(xP+1)×ρ(x1x2P)eiϕ(x1x2P),
(4)

where ρ(x1x2P) is a positive-definite distribution given by

ρ(x1x2P)=mP2π|τc|Pexpk=12PmPβ4|τc|2xk+1xk2+β2PU(xk)x2P+1=x1,
(5)

and the phase factor ϕ(x1x2P) is defined to be

ϕ(x1x2P)=mPt2|τc|2k=1Pxk+1xk2k=P+12Pxk+1xk2tPk=2PU(xk)k=P+22PU(xk).
(6)

The condition x2P+1 = x1 arises from the trace condition. Consequently, the path integral we obtain involves a ring polymer consisting of 2P beads. As Eqs. (4)–(6) show, both the positive definite distribution and the phase factor depend on β and t; they also depend on the same set of path-integral variables x1, …, x2P, which has some advantages for numerical convergence over the standard correlation function CAB(t), as it allows the phase factor to be treated as a weight in the configurational average. However, the rapid fluctuations of the phase factor caused by small changes in the complex-time paths still cause severe convergence issues.

Since the positive-definite distribution and phase factor in the symmetrized correlation function involve the same set of path-integral bead variables, a useful way to tame the phase factor is to introduce a change of the integration variables x1, …, x2P to new variables that represent an average of the forward and backward complex-time paths and the difference between these two paths. One possibility for such a transformation is the variable set ri = (xi + xP+i)/2, si = xixP+i, where i = 1, …, P. In this transformation, each point along the forward path is combined with the point along the reverse path having the same imaginary time index. This transformation has a few appealing properties, in particular, since x1 = x2P+1, it follows that rP+1 = r1 and sP+1 = −s1. Thus, if we expand the potential appearing in ρ and ϕ to second order in s1, …, sP and then integrate analytically over s1, …, sP, the remaining path integral over r1, …, rP will involve a ring polymer of P beads in an effective potential with contributions arising from the analytical integration over the path difference variables. However, as the impetus for expanding the potential in powers of si is that these variables are small, the unfortunate fact remains that the points xi and xP+i are separated by P intervening beads, which gives them a maximal separation in the original 2P-bead ring polymer so that the assumption under which we expand the potential in powers of si is not valid. Consider, however, the following transformation:24–26,28–35

r1=x1,rP+1=xP+1,ri=xi+x2(P+1)i2,  i=2,,P,si=xix2(P+1)i,  i=2,,P.
(7)

The variables ri and si are also path average and difference variables, respectively; however, in this transformation, early points along the forward path are matched with corresponding late points along the reverse path and vice versa, as illustrated in Fig. 1. This transformation has the following inverse:

x1=r1,xP+1=rP+1,xi=ri+12si,  i=2,,P,x2(P+1)i=ri12si,  i=2,,P.
(8)

In many condensed-phase problems, there is significant decoherence, suggesting that the path difference variables si can be expected to be small. Thus, the transformation in Eq. (7), which makes these variables explicit, allows the resulting path integral to be expanded in powers of si, an expansion that would be expected to converge quickly as these powers increase. For practical implementation, if the expansion is truncated at second order, then the resulting integrals over si can be performed analytically. When this is done, the resulting path integral over ri will involve an open or worm-like chain. Under this transformation, the harmonic bead-bead coupling terms in the distribution function ρ can be easily worked out, and we find the following:

k=12PxkxK+12  =i=1Pxixi+12+i=1Px2(P+1)i1x2(P+1)i2  =2i=1Pri+1ri2+12i=2P1si+1si2+12s22+sP2.
(9)

Similarly, in the phase factor, we obtain

k=1Pxk+1xk2k=P+12Pxk+1xk2  =i=1Pxi+1xi2i=1Px2(P+1)i1x2(P+1)i2  =2i=2P1ri+1risi+1si+2s2r2r1   2sPrP+1rP.
(10)

For the potential energy, we perform a power series expansion about si = 0 on the two types of terms that result from the transformation,

U(xi)=Uri+12siU(ri)+12U(ri)si+18U(ri)si2,U(x2(P+1)i)=Uri12siU(ri)12U(ri)si+18U(ri)si2.
(11)

A similar second-order expansion of the potential around the difference variables has been implemented in the local harmonic approximation/local Gaussian approximation (LHA/LGA) for the linearized semiclassical initial value representation (LSC-IVR) methods.49,50 The justification for the second-order expansion of the potential is to evaluate appropriate integrals analytically as shown below. Thus, in the positive-definite distribution, the potential energy becomes

k=12PU(xk)=i=1PU(xi)+i=1PU(x2(P+1)i)U(r1)+U(rP+1)+2i=2PU(ri)+14i=2PU(ri)si2,
(12)

and in the phase factor,

k=2PU(xk)k=P+22PU(xk)=i=2PU(xi)i=2PU(x2(P+1)i)=i=2PU(ri)si.
(13)

Putting the potential and harmonic coupling expressions together, ρ transforms to

ρ(r1rP+1,s2sP)=mP2π|τc|PemPβ4|τc|22i=1Pri+1ri2+12i=2P1si+1si2+12(s22+sP2)×eβ2Pi=1Pi=2P2U(ri)+14U(ri)si2+U(r1)+U(rP+1),
(14)

and the phase factor becomes

ϕ(r1rP+1,s2sP)=mPt2|τc|22i=2P1(ri+1ri)(si+1si)+2s2(r2r1)2sP(rP+1rP)tPi=2PU(ri)si.
(15)
FIG. 1.

Schematic representation of the transformation of Eq. (7). Specifically, a closed ring polymer of size 2P where P = 4 is transformed to an open chain of P + 1 path integral beads. The forward complex-time path is represented by blue and the backward complex-time path is represented by red. The transformation matches (dotted lines) each early point along the forward path with a corresponding late point along the backward path.

FIG. 1.

Schematic representation of the transformation of Eq. (7). Specifically, a closed ring polymer of size 2P where P = 4 is transformed to an open chain of P + 1 path integral beads. The forward complex-time path is represented by blue and the backward complex-time path is represented by red. The transformation matches (dotted lines) each early point along the forward path with a corresponding late point along the backward path.

Close modal

Note that the distribution can be separated into a purely ri-dependent term and a term that depends on both ri and si according to

ρ(r1rP+1,s2sP)=ρ1(r1rP+1)ρ2(r2rP,s2sP),
(16)

where

ρ1(r1rP+1)=mP2π|τc|PemPβ4|τc|2i=1P2(ri+1ri)2×eβ2Pi=2P2U(ri)+U(r1)+U(rP+1)
(17)

and

ρ2(r2rP,s2sP)=emPβ4|τc|212i=2P1si+1si2+12(s22+sP2)×eβ2Pi=2P14U(ri)si2.
(18)

We now define a square tridiagonal (P − 1) × (P − 1) matrix M as

Mij(r)=2A+β4PU(ri)δijAδi,j+1Aδi+1,j,
(19)

where rr2, , rP, and

A=mPβ4|τc|2,
(20)

and a (P − 1)-dimensional vector S = (s2, …, sP). In terms of these, the distribution ρ2 can be rewritten as

ρ2(r2rP,s2sP)=e12STM(r)S.
(21)

The phase factor can be rewritten in a similar manner by introducing a (P − 1)-dimensional vector K defined as

Ki(r)=γ2riri1ri+1tPU(ri),  i=2,,P,
(22)

where r′ ≡ r1, …, rP+1, and

γ=mPt|τc|2.
(23)

In terms of K and S, the phase becomes

ϕ(r1rP+1,s2sP)=K(r)S.
(24)

Thus, the overall expression for GAB(t) within our approximation becomes

GAB(t)=1Zdr1drP+1ρ1(r1rP+1)ar1brP+1×ds2dsPe12STM(r)S+iK(r)S.
(25)

Within this approximation, the integration over S is simply a Gaussian integral which can be performed analytically, and while the second-order expansion and resulting integration are not new, the motivation of the expansion in this derivation allows for the elimination of the underlying imaginary terms in the integrals under the assumption that the eigenvalues of M are all positive definite. This last condition should be obeyed if P is sufficiently large, which contrasts with the conditions obtained in trajectory-based methods.31 The result of the integration is

ds2dsPe12STM(r)S+iK(r)S=2πP1det[M(r)]12e12K(r)TM1(r)K(r),
(26)

and the symmetrized correlation function that results involves an imaginary-time path integral over open chains in the variables r1, …, rP+1 is

GAB(t)=2π(P1)/2Zdr1drP+1ρ1(r1rP+1)ar1brP+1×e12KT(r)M1(r)K(r)+lndet[M(r)].
(27)

Eq. (27) indicates that the open chain is subject to an effective or modified potential determined not only by the original potential U but also by the matrix M(r) and the vector K(r′). However, as Eq. (27) is also just an imaginary-time path integral, any of the usual sampling schemes typically applied for path integrals can be employed, including Monte Carlo,43 molecular dynamics,13 hybrid Monte Carlo,13,44 or any path-integral enhanced sampling scheme.45,46 Here, although it might not be optimal, for simplicity, we have employed a standard PIMD scheme. We thus introduce fictitious momenta and corresponding fictitious masses and write Eq. (27) as

GAB(t)=1Zdp1dpP+1dr1drP+1a(r1)b(rP+1)eβH̃(r,p),
(28)

where pp1, …, pP+1, and where the bead Hamiltonian is given by

H̃=i=1P+1pi22m+i=1P12mωP2(ri+1ri)2+1Pi=2PU(ri)+12PU(r1)+U(rP+1)+Ṽ(r1,,rP+1).
(29)

Every term in this Hamiltonian, with the exception of the last term and ωp=P/|τc|, is identical to standard open PIMD.51–53 The last term is an additional potential resulting from the Gaussian integration and is defined as

Ṽ(r1,,rP+1)=12βKT(r)M1(r)K(r)+ln(det[M(r)]).
(30)

Note that, unlike U, the potential (r1, …, rP+1) couples all of the beads, i.e., it cannot be written as a sum over terms that are functions of only one bead variable. It also depends on both β and t.

At this point, a few comments are in order. First, since the approach is based on an open-chain path integral formulation of the symmetrized correlation function, we assign the acronym OPSCF to the method, which stands for Open Path Symmetrized Correlation Function. Second, although we have truncated the expansion in the path difference variables at second order, the scheme can be systematically improved by going to higher orders in si and employing a cumulant expansion. In  Appendix A, we derive an example, computing to the first cumulant coming from an expansion up to third order in si. Importantly, the expansion and subsequent analytical integration that lead to Eq. (30) involve an expansion of the potential only about si = 0 rather than a global expansion of the potential as in Ref. 25. Thus, the full potential dependence on ri is retained along with the extra terms in the effective potential appearing in Eq. (30). The quality of this effective potential determines the accuracy of the approximate correlation functions produced by OPSCF. Third, and following on the previous point, it is possible to improve the accuracy with the choice of P by, for example, designing colored-noise thermostatting approaches54,55 or by using a higher-order Trotter expansion, for example, the Suzuki-Chin approximation.47,48 Note that these schemes cannot be used to accelerate RPMD or CMD calculations, indicating a possible advantage of sampling schemes over trajectory-based methods. In  Appendix B, we derive the corresponding Suzuki-Chin version of Eq. (27). The result, which is manifestly real, is the appearance of new, additional terms in the effective potential in Eq. (30). Fourth, the operators in Eq. (27) need not be linear functions of position. Since the transformation does not affect the a(x1) = a(r1) and b(xP+1) = b(rP+1) terms in the integral, the approximation should be of similar quality for both linear and nonlinear operators. Fifth, as previously mentioned, the path integral in Eq. (27) can be evaluated by any type of sampling algorithms, including those designed for Monte Carlo calculations.56 Sixth, the value of t is treated as an initial parameter in the simulation, and thus each value of t requires a separate simulation regardless of the sampling scheme. Finally, for general potential energy functions, the canonical ensemble is not preserved within the OPSCF scheme, and we will discuss this issue shortly.

The OPSCF method approximates a 2P ring polymer by a P + 1 open chain in an effective potential. This should not be surprising as the symmetrized correlation function, in contrast to the standard correlation function CAB(t), is constructed such that the same bead variables define both the positive-definite distribution and the phase function. Consequently, any approximation made to tame the phase factor must be made to the positive-definite (canonical) distribution as well. Thus, the validity of the method relies on the accuracy of the second-order expansion of the potential made in both contributions so the correlation function. Note that the canonical ensemble distribution is generated only at t = 0 in the symmetrized correlation function, and since the phase factor vanishes at t = 0, it is always possible to treat t = 0 as a special case for which an exact 2P-bead cyclic path integral calculation is performed without making the OPSCF transformation and subsequent power-series approximation. If this is done, then the canonical ensemble will be preserved at t = 0. In the applications to be presented below, however, we will not employ this trick, as it is important to assess the degree to which OPSCF approximates the true canonical distribution at t = 0. The OPSCF method does preserve detailed balance, however, which is also true for CMD and RPMD but is not true for some LSC-IVR implementations. The only changes to the derivation needed to prove preservation of detailed balance are the sign of the imaginary term in Eq. (25) and an exchange of the bead indices at which operators are evaluated, i.e., [b(r1) and a(rP+1)]. Clearly, the change of sign does not affect the result of the analytical integration of Eq. (25) over s2, …, sP. Moreover, the symmetry of the open chain allows the evaluation of the operators at either endpoint as long as both endpoints sample the entire phase space. It must be noted, however, that the argument made above assumes the validity of the positive-definite distribution ρ(r1, …, rP+1) of the OPSCF approximation. However, as this distribution is only exact for the harmonic oscillator, preservation of detailed balance in anharmonic systems will only be valid with respect to this (approximate) distribution.

All of these characteristics need to be contrasted to the CMD and RPMD approaches, in which it is not possible to perform the calculations by any means other than molecular dynamics, higher-order factorizations cannot be used, and which degrade when applied to nonlinear operators. The fact that the OPSCF scheme presented here can be improved in a systematic and practically computable way also stands in contrast to CMD and RPMD for which it is not straightforward to see how to obtain systematically better results than what the basic schemes produce. Furthermore, the OPSCF transformation and subsequent power-series expansion are not explicitly limited to the symmetrized correlation function but could also be applied to the Kubo-transformed correlation function. The OPSCF method should also be contrasted with trajectory-based schemes, which also require the generation of dynamical trajectories and carry a computational overhead when attempting to accurately capture quantum coherence. Finally, we point out that the OPSCF scheme is exact for harmonic potentials, for which the expansion in si is exact, and converges to the correct classical limit, where the difference between the forward and backward path vanishes. The obvious downside is the need to compute the second derivative of the potential.

In Secs. III and IV, we will examine a number of simple one-dimensional examples that are commonly used to study the performance of the CMD, RPMD, and semiclassical approaches. We first describe the computational details of our simulations followed by a presentation of the results for each of the examples.

Two correlation functions on which we will focus in this work are the position autocorrelation function (Â(x^)=B^(x^)=x^) and the position-squared autocorrelation function (Â(x^)=B^(x^)=x^2). We explored four different potentials using the Open Path Symmetrized Correlation Function (OPSCF) method. These are a harmonic oscillator

U(x)=12x2,
(31)

a mildly anharmonic potential

U(x)=1100x4+110x3+12x2,
(32)

a quartic potential

U(x)=14x4,
(33)

and a symmetric double well potential

U(x)=110x2522.
(34)

These various potentials have been extensively used to benchmark approximation schemes for computing quantum time correlation functions.7,10,17,18,32,33,57–61 In all examples, we use natural units in which = 1. Two different temperatures were explored for each potential. Specifically, β = 1 and β = 8, except for the symmetric double well, in which the low temperature was β = 4. The position-squared correlation function, Gx2x2(t), was excluded from the main body of the text for the mildly anharmonic potential and quartic potential. We found that Gx2x2(t) for these potentials at the lowest temperatures considered here were so small as to be within the sampling error that no useful information concerning the performance of the method could be discerned.62,63 At the higher temperatures, however, we were able to obtain quality Gx2x2(t) curves that we could benchmark against exact results. Section IV will discuss the particular issues encountered in generating Gx2x2(t) at low temperatures, specifically for the harmonic oscillator and double well potentials.

We implemented the stochastic isokinetic thermostat64 scheme solely to ensure proper canonical sampling. The OPSCF simulations were performed using staging variables13 for open chains.53 The time step of the simulations ranged from 0.005 to 0.05 natural units of time depending on the temperature, the value of t, and the bead number P. We chose the value of P such that |τc|/P remains approximately constant as t is increased, and as P is increased, the time step must typically be decreased, or multiple time-step methods can be employed.16,81 Figure S1 of the supplementary material demonstrates convergence of the correlation functions with respect to the bead number. Simulations were run for a minimum of 107 steps for β = 1 and 2 × 107 steps for β = 4, 8 until proper convergence is reached. Standard errors of the data points of each simulation are depicted as error bars in various figures. Finally, the Suzuki-Chin factorization for open path-integral molecular dynamics (OPIMD)47,48,65,66 (see  Appendix B) was implemented for the mildly anharmonic and quartic potentials (Figs. S2–S5 of the supplementary material) in order to determine if a higher-order factorization scheme improved the efficiency by allowing smaller values of P.

We benchmarked the results of the OPSCF method against RPMD correlation functions as well as correlation functions obtained from the linearized semiclassical initial value representation (LSC-IVR). For the RPMD simulations, 100 000 configurations were accumulated from a long thermostatted PIMD simulation of 108 steps for the various potentials and temperatures. From these configurations, we sampled the initial velocities from the Maxwell-Boltzmann distribution. Trotter numbers of 8 and 32 were used for the high and low temperatures, respectively. From each of these initial phase-space points, we performed the RPMD simulations in primitive coordinates with a time step of 0.0005 natural units of time to converge the dynamics fully, following the protocol of Ref. 57. Each RPMD simulation was run for 80 000 steps, which is the equivalent to 40 natural units of time. Since RPMD produces the Kubo transformed correlation functions, KAB(t), the symmetrized correlation function for RPMD was computed by performing a Fourier transform of KAB(t) and generating the Fourier transform of GAB(t) via

G̃AB(ω)=βω2sinh(βω2)K̃AB(ω),
(35)

followed, finally, by an inverse Fourier transform to give the RPMD approximation to GAB(t). Figures S6–S9 in the supplementary material provide the Kubo transformed correlation functions directly from the RPMD simulations prior to their transformation to GAB(t).

For the LSC-IVR calculations, the thermal Gaussian approximation (TGA-LSC-IVR) was used following the scheme devised by Liu and Miller.27 Specifically, 41 imaginary time trajectories were generated for the potentials at the desired temperatures. For these 41 trajectories, the differential equations were solved using a fourth-order Runge-Kutta algorithm. Subsequently, approximately 107 real time trajectories were sampled from each imaginary time trajectory and the dynamics were evolved using the velocity Verlet integrator with a time step of 0.02. TGA-LSC-IVR generates the symmetrized autocorrelation function, GAA(t), directly. While TGA-LSC-IVR might not be the most robust semiclassical method, it has been shown to perform well in condensed-phase systems28 and is an appropriate semiclassical method to compare to OPSCF. Exact quantum time correlation functions were calculated by solving the eigenvalues and eigenvectors of each system using a harmonic basis set, which was performed in Mathematica.67 In these calculations, 100 basis functions were sufficient to converge all exact quantum time correlation functions.

The OPSCF method is exact for harmonic systems, as is also the case for RPMD, CMD, and TGA-LSC-IVR. Therefore, simulations of the harmonic oscillator provide a benchmark for understanding the strengths and weaknesses associated with numerical evaluation of Eq. (27). Figure 2(a) depicts the position autocorrelation function Gxx(t) at β = 1. RPMD, TGA-LSC-IVR, and OPSCF reproduce the exact result as expected. Figure S1 of the supplementary material shows the convergence of the OPSCF correlation functions as a function of the number of beads. In Sec. III, we discussed the convergence criterion for OPSCF calculations, specifically, that the value |τc|/P should be kept approximately constant. In this way, the proposed criterion resembles the convergence criteria often used in imaginary-time path integral simulations (i.e., β/P) remain constant as β is varied. In OPSCF, however, |τc| depends on both time and temperature. Thus, each of the data points in the correlation function plots to be presented in this section corresponds to the minimum value of P required to reach convergence.

FIG. 2.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for the harmonic oscillator, Eq. (31), for β = 1. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (dark to orange: P = 8, P = 16, P = 32).

FIG. 2.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for the harmonic oscillator, Eq. (31), for β = 1. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (dark to orange: P = 8, P = 16, P = 32).

Close modal

Figure 2(b) demonstrates the increased difficulty of obtaining correlation functions for non-linear operators over linear operators. This problem has challenged schemes such as CMD and RPMD, and various approaches have been introduced to address it.7,14–18 Specifically, Fig. 2(b) shows the position-squared autocorrelation function Gx2x2(t) at β = 1. RPMD misses some of the extrema, particularly the minima, of this correlation function. The OPSCF method and TGA-LSC-IVR, by contrast, are seen to reproduce the exact result within statistical error bars. Equation (27) is agnostic to the functional form of operators  and B^ as long as the operators are solely functions of positions. The accuracy of Gx2x2(t) suggests that the OPSCF method should reproduce any linear and nonlinear correlation function for the harmonic oscillator as long as the configuration space is sufficiently sampled. This highlights the strength of the OPSCF approximation in representing correlation functions for nonlinear operators.

Reducing the temperature of the harmonic oscillator ensures that higher energy states are less accessible, thus creating a ground-state dominance, thereby increasing the quantum character of the system. In principle, since the OPSCF method is exact for the harmonic oscillator, we should be able to reproduce the exact correlation functions for both linear and nonlinear operators. Gxx(t) [Fig. 3(a)] is reproduced very well with OPSCF, RPMD, and TGA-LSC-IVR although a larger value of P is needed for the OPSCF and RPMD compared to β = 1. Some of the data points lie just outside the exact curve but are within statistical accuracy denoted by the error bars.

FIG. 3.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for the harmonic oscillator, Eq. (31), oscillator for β = 8. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (P = 32, orange; P = 64, dark yellow).

FIG. 3.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for the harmonic oscillator, Eq. (31), oscillator for β = 8. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (P = 32, orange; P = 64, dark yellow).

Close modal

Finally, we discuss the efficiency of the OPSCF method in reproducing Gx2x2(t) for the harmonic oscillator at β = 8 [Fig. 3(b)]. As stated, the OPSCF method is exact in the harmonic oscillator and should reproduce any position-dependent correlation function. However, Fig. 3(b) demonstrates that while many points are within statistical error of the exact correlation function, the data points do not produce the correct periodicity compared to the exact curve. Given that the oscillations for the Gx2x2exact(t) range within an order of 10−4, the level of sampling is not able to resolve such small differences. Note that while the periodicity of RPMD is quite close to the exact result, the curve is shifted considerably from the exact curve. Interestingly, however, this numerical resolution issue is less severe for the Kubo-transformed correlation function, Kx2x2(t) (Fig. S6 of the supplementary material), which RPMD approximates directly. As Fig. S6(d) of the supplementary material shows, the amplitude of the oscillations is larger than that of Gx2x2(t). Furthermore, the results from TGA-LSC-IVR are excluded, as the evaluation of the symmetrized correlation function is numerically unfavorable for this method. However, it has been shown that LSC-IVR is able to recover the correlation function Cx2x2(t) exactly for the harmonic oscillator at low temperatures.62,63 This issue regarding Gx2x2(t) at low temperatures is present in both the mildly anharmonic potential and the quartic potentials. Therefore, we have excluded the data for Gx2x2(t) at β = 8 from this present work owing to the difficulty of converging the curve to a sufficiently small numerical error that a meaningful comparison could be made. We will, however, continue to present Gx2x2(t) for β = 1.

The mildly anharmonic potential, Eq. (32), and its associated time correlation functions have been extensively studied as a nontrivial example to test approximation schemes.10,17,18,32,57–59 The slight anharmonicity in the potential produces small deviations from the simple harmonic oscillator. We can, therefore, test the validity of the second-order truncation of the potential energy within the OPSCF method using this anharmonic potential. Once again, at β = 1, Fig. 4(a), all three methods reproduce the correlation function, Gxx(t). All methods exhibit a very slight shift in frequency at longer times. When comparing to the harmonic oscillator, the OPSCF requires a larger value of P to accurately capture the autocorrelation function.

FIG. 4.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for mildly anharmonic potential, Eq. (32) for β = 1. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (dark red to dark yellow: P = 8, P = 16, P = 32, P64).

FIG. 4.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for mildly anharmonic potential, Eq. (32) for β = 1. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (dark red to dark yellow: P = 8, P = 16, P = 32, P64).

Close modal

In Fig. 4(b), we see that the OPSCF does remarkably well for the nonlinear case of Gx2x2(t) compared to RPMD, and OPSCF compares well against TGA-LSC-IVR. In particular, RPMD accurately captures some of the extrema of the curve but fails at the minima, while OPSCF outperforms RPMD at many of the extrema. So far, the OPSCF method has proven to be an effective tool at the higher temperature β = 1 even for Gx2x2(t). The Suzuki-Chin factorization implemented for this potential at β = 1 [see Fig. S2(b) of the supplementary material] improves the performance of Gx2x2(t) compared to the traditional Trotter factorization. Specifically, smaller values of P are needed to converge some of the individual data points.

At β = 8, we begin to see the limitations of all three methods. Interestingly, Fig. 5 demonstrates that OPSCF produces the correct frequency in the curve but with a shift in the amplitude, while for RPMD and TGA-LSC-IVR, there is a noticeable shift in the frequency, but the amplitude is more accurate. Furthermore, one may notice that the RPMD curve does not begin at the correct value of Gxx(0). This is a result of the Fourier transform described in Eq. (35). Figure S7(b) of the supplementary material demonstrates the accuracy of the RPMD simulations compared to the exact Kubo-transformed correlation functions, Kxx(t), for the mildly anharmonic potential. Even TGA-LSC-IVR does not have the correct t = 0 limit in this case. However, the semiclassical method is able to recover the correct t = 0 limit for anharmonic potential using different correlation functions,62,63 thus demonstrating the difficulty in converging the symmetrized correlation function at low temperatures. The OPSCF method is able to capture the correct t = 0 limit at this temperature, but it begins to degrade soon after. Even a slight anharmonicity in the potential will be affected by the truncation of the potential to second order in Eq. (11). In  Appendix A, we propose an expansion of the potential up to third order in si allowing us to calculate the first cumulant. In theory, this would improve the results, specifically for t > 0. As stated in  Appendix A, this cumulant will be incorporated into the time-dependent phase factor, and thus t = 0 would not be affected by the first cumulant. This is ideal for the mildly anharmonic potential since the short time limit is accurately obtained for the OPSCF method. The Suzuki-Chin implementation of the OPSCF method for β = 8, Fig. S3 of the supplementary material, provides moderate improvement to Gxx(t), especially within the first few data points.

FIG. 5.

Symmetrized position autocorrelation function for mildly anharmonic potential, Eq. (32), for β = 8. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (P = 32, red; P = 64, dark yellow). Black dotted line used to clarify the frequency of the correlation function from the OPSCF method.

FIG. 5.

Symmetrized position autocorrelation function for mildly anharmonic potential, Eq. (32), for β = 8. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (P = 32, red; P = 64, dark yellow). Black dotted line used to clarify the frequency of the correlation function from the OPSCF method.

Close modal

Generating accurate time correlation functions for the quartic potential, Eq. (33), has been difficult for most approximation schemes.6,10,32,33,58–60 None of the methods presented replicate the long time behavior of the correlation function. In fact, RPMD and OPSCF methods are intriguingly similar, while TGA-LSC-IVR falls short. Gxx(t) at β = 1, Fig. 6(a), for RPMD and OPSCF decay soon after the first oscillation. However, neither extrema of the first oscillation are perfectly captured. Most current methods fail to reproduce multiple oscillations. At a slightly lower temperature, β = 2, Matsubara dynamics68,69 are able to reproduce more oscillations compared to RPMD and CMD for the quartic potential. While Matsubara dynamics are computationally expensive, its basis in semiclassical theory may provide integral information in improving the OPSCF method to produce multiple oscillations for the quartic potential.

FIG. 6.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for purely quartic potential, Eq. (33), for β = 1. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (dark red to dark yellow: P = 8, P = 16, P = 32, P = 64).

FIG. 6.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for purely quartic potential, Eq. (33), for β = 1. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (dark red to dark yellow: P = 8, P = 16, P = 32, P = 64).

Close modal

The OPSCF method outperforms both RPMD and TGA-LSC-IVR for Gx2x2(t) (Fig. 6) by locating the correct extrema in the first oscillation. This has not been achieved either with RPMD or with CMD nor have improvements to CMD, such as the use of CMD in conjunction with maximum entropy analytical continuation16 (CMD-MEAC) succeeded in this regard. Thus far, one of the strengths of OPSCF is the accuracy of generating nonlinear quantum correlation functions for β = 1. Capturing the early extrema in Gx2x2(t) is a positive outcome for the OPSCF method even if the OPSCF method cannot extend past the first oscillation; the short time behavior of the correlation function is well reproduced.

The low-temperature regime still proves to be a challenge for the OPSCF method as currently implemented (Fig. 7). The low temperature of β = 8 demonstrates the first instance in this paper when GxxOPSCF(0)Gxxexact(0). t = 0 is a special case within the OPSCF scheme when the phase factor in Eq. (4) disappears resulting in a 2P imaginary-time path integral. Thus, by performing the transformation of Eq. (7), we are effectively approximating a closed-chain imaginary time path integral with an open-chain imaginary time path integral containing a modified form of the additional potential of Eq. (30), Ṽ=12βln(det[M(r)]), which includes the approximate second-order expansion of the potential energy around si. In other words, the t = 0 statistics are approximated as they are at other times, and this stands in contrast to other semiclassical approaches that correctly reproduce the equilibrium statistics. This is a consequence of the use of the symmetrized correlation function, as the same bead variables appear in both the positive-definite density and the phase factor, and hence, approximation made to the latter also affects the former. The data at t = 0 confirm that the truncation is insufficient for the quartic potential at low temperatures, while at β = 1, quantum effects are less pronounced and we find the truncation to be sufficient in the short time limit. Note, however, that the t = 0 value is still a slight improvement over RPMD, which also fails to produce the correct t = 0 value of the correlation function. TGA-LSC-IVR interestingly correctly captures the t = 0 limit even though it failed to do so for the anharmonic potential. In a similar manner to the mildly anharmonic potential, a cumulant expansion provides a systematic improvement to the method. However, it is clear that a higher-order cumulant would be necessary (not shown in  Appendix A) in order to bring in a correction term at t = 0. Moreover, as previously noted, t=0 can be treated as a special case, and the exact 2P-path integral could be performed without approximation in order to obtain the correct value at this point.

FIG. 7.

Symmetrized correlation function for the position autocorrelation function for purely quartic potential, Eq. (33), for β = 8. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (P = 32, orange; P = 64, dark yellow). Black dotted line used to clarify the frequency of the correlation function from the OPSCF method.

FIG. 7.

Symmetrized correlation function for the position autocorrelation function for purely quartic potential, Eq. (33), for β = 8. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (P = 32, orange; P = 64, dark yellow). Black dotted line used to clarify the frequency of the correlation function from the OPSCF method.

Close modal

For t > 0, the OPSCF method also compares poorly to the exact correlation function, Gxx(t), in the low temperature regime. The data points fail to describe both the amplitude and frequency of the correlation function. This system would need to be evaluated more extensively to determine if and how a systematic improvement could be developed. One such approach could be the incorporation of the first- and second-order cumulant expansions. As stated above, at t = 0, the second-order expansion of the potential energy is insufficient in describing the quartic potential at low temperatures, and this would hold true for t > 0. Even though, TGA-LSC-IVR correctly obtains the short time limit, this method decays too rapidly compared to RPMD. However, a more robust semiclassical method would presumably perform better at longer times. RPMD does a much better job in the low temperature regime over a long time compared to β = 1. It has been argued by Voth and co-workers6,58,60 that the low temperature regime of the quartic oscillator creates a two-state system that effectively evolves in a quasi-harmonic environment. Both RPMD and CMD are exact in the harmonic limit and thus will perform reasonably well in this regime. However, the dynamics of the higher frequency modes in the RPMD scheme cause the dampening effect of the correlation function further out in time. Although the OPSCF method is also exact in the harmonic limit, it appears to be more sensitive to the deviations produced by the quartic potential. Unlike the mildly anharmonic potential, the Suzuki-Chin factorization (Figs. S4 and S5 of the supplementary material) does not appear to affect the convergence of the correlation function with respect to P.

The symmetric double well potential, Eq. (34), is the final potential used to test the OPSCF method and has been explored for a variety of the methods currently discussed.58,61–63 Not only does this potential deviate greatly from the harmonic oscillator but also it has been studied in both enhanced sampling methods45 and quantum transition state theory,70 specifically in understanding quantum tunneling.

Figure 8(a) shows the position correlation function in the short time regime for β = 1. All methods perform similarly for Gxx(t). However, OPSCF improves when examining Gx2x2(t) [Fig. 8(b)] confirming the OPSCF method’s ability to determine accurately the nonlinear correlation function for the higher temperature (β = 1) compared to RPMD and TGA-LSC-IVR. While the higher temperature regime favors OPSCF, the method continues to struggle with the lower temperatures. Figure 9 illustrates the limitations of the method at β = 4. An intermediate temperature was used instead of β = 8 to ensure that sufficient sampling over the barrier could be achieved. The OPSCF method reproduces the correlation function Gxx(t) at short time only [Fig. 9(a)], while Gx2x2(t) [Fig. 9(b)] is rather poor at all times. By contrast, RPMD captures the short time limit of both correlation functions rather well. TGA-LSC-IVR is excluded for this temperature since convergence of the symmetrized correlation function was poor, meaning that we could not run enough imaginary-time or real-time trajectories to yield converged results. However, TGA-LSC-IVR has performed well in the short time limit for the double well potential at low temperatures for the Kubo and standard correlation functions.62,63 Systematic improvements to the OPSCF method, specifically first-and second-order cumulants, could lead to more accurate results.

FIG. 8.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for symmetric double well potential, Eq. (34), for β = 1. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (P = 8, dark red).

FIG. 8.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for symmetric double well potential, Eq. (34), for β = 1. Exact, black; RPMD, blue; TGA-LSC-IVR, green; OPSCF, dots (P = 8, dark red).

Close modal
FIG. 9.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for symmetric double well potential, Eq. (34), for β = 4. Exact, black; RPMD, blue; OPSCF, dots (P = 32, orange).

FIG. 9.

Symmetrized correlation function for the position (a) and position-squared (b) autocorrelation functions for symmetric double well potential, Eq. (34), for β = 4. Exact, black; RPMD, blue; OPSCF, dots (P = 32, orange).

Close modal

Another issue arises with the OPSCF method in regards to the symmetric double well potential. Specifically, the OPSCF simulations for the double well potential begin to fail for times t > 3. A preliminary analysis of the forces produced by the simulations revealed that the forces spike to large values (many orders of magnitude) throughout the simulation. This alters the dynamics and the simulations no longer follow the equations of motions defined by the Hamiltonian of Eq. (29). We are currently investigating the nature of this behavior and if either an alternative sampling scheme or other improvements can be made to alleviate this problem. Monte Carlo would provide an alternative sampling mechanism if the evaluation of forces proves to be the cause of the failed simulations.

At this point, a comment is in order on how one might use OPSCF to compute a quantum rate since rates require the evaluation of a different nonlinear correlation function. The alternative methods of RPMD,71,72 CMD,73 and LSC-IVR50,74 have the ability to calculate rate constants by utilizing the flux-side correlation function. The symmetrized version of this correlation function is

Gfs(t)=1ZTrF^eiĤτc*/ĥeiĤτc/,
(36)

where F^=12 mp^δ(x^(0)x)+δ(x^(0)x)p^ is the flux operator and ĥ=h(x^(t)x) is the Heaviside function, and x is the location of the transition state (dividing surface) such that the reaction rate is

k(β)=limt1Qr(β)Gfs(t),
(37)

where Qr(T) is the reactant partition function. In a similar fashion to Eq. (3), the symmetrized version of the flux-side correlation function can be written as

Gfs(t)=1Zdxdxdxdpx|eiHτc*/|xx|eiHτc/|x×h(x)peip(xx)12mδ(xx)+δ(xx).
(38)

The incorporation of a momentum-dependent operator (flux operator) yields a few more terms in Eq. (38) compared to Eq. (3). Path integral representations of symmetrized correlation functions of momentum operators do not result in a closed ring polymer [see Eq. (3) and Fig. 1 for position operators], but instead an open chain is formed. Therefore, the transformation of Eq. (7) from an open chain to a different open chain and the second-order expansion of the potential become much less straightforward for momentum operators. However, calculating momentum operators can be accomplished by taking an appropriate time-derivative of position operators. In the case of the flux-side correlation function, a reaction rate can be expressed by a time derivative of the side-side correlation function

k(β)=limt1Qr(β)ddtGss(t).
(39)

The symmetrized side-side correlation function,

Gss(t)=1Z TrĥeiĤτc*/1ĥeiĤτc/,
(40)

is solely dependent on position operators and can be used within the OPSCF method resulting in an estimator for the side-side correlation function,

Gss(t)=1Zdp1dpP+1dr1drP+1h(r1x)×1h(rP+1x)eβH̃(r,p).
(41)

The side-side correlation function has not been explored extensively in the literature outside its relation to the flux-side correlation function. However, it has been shown that the symmetrized side-side correlation function can be analytically solved for a simple parabolic barrier.75 Since the transformation and second-order expansion of the potential described in Sec. II are exact for quadratic potentials, one can recover analytically the side-side symmetrized correlation function for the parabolic barrier within the OPSCF method. For more complex systems, Eq. (41) is a proposed sampling scheme to determine side-side correlation functions necessary for reaction rates. Written in this form, the estimator yields an interesting picture such that the estimator is non-zero when one end of the open chain is on the “product” side of the transition state and the other end is on the “reactant” side. While reaction rates can be calculated in principle, further exploration of the OPSCF method in relation to quantum rates is necessary in order to understand any limitations as well as improvements to the current derivation of the method.

We make one final point on the OPSCF method pertaining to the double well potential. An interesting advantage of the OPSCF method is the potential use of enhanced sampling methods to obtain quantum time correlation functions. As stated previously, Eq. (27) can use any preferred sampling technique to produce quantum time correlation functions. Since the observables of interest, a(x1) and b(xP+1), are functions of the chain endpoints, enhanced sampling methods can be used to ensure that the end points of the open chain are sampled sufficiently. Furthermore, we can propose an enhanced sampling technique over the entire configuration space in systems where the open chain cannot thermally access regions of phase space, not just the end to end distances. There is an emerging line of research focused on the development of enhanced sampling of “rare” events in quantum mechanics.45,46 While this has not been explored in the present work, quantum time correlation functions of “rare” events would be an interesting area of research to explore, which neither CMD or RPMD can achieve on their own.

Given the challenge of computing quantum time correlation functions numerically, new perspectives on this problem are needed. In this paper, we have proposed an approach that draws inspiration from semiclassical methodology but that aims to work within an imaginary-time bead sampling paradigm rather than a trajectory-based one. In so doing, we have generated a scheme that can leverage all of the available sampling and enhanced sampling schemes applicable for imaginary-time path integrals, thereby offering the potential for considerable savings over trajectory-based formulations. The approach exploits a transformation that is often used in semiclassical schemes, namely, to sums and differences of the forward and backward complex-time paths of the symmetrized correlation function GAB(t), together with an expansion of the potential in powers of the difference variables. The transformation is applied, in the present context, in such a way that, after analytical integration over the difference variables, the result is an open-chain imaginary-time path integral subject to a modified potential that couples all of the path-integral beads. Therefore, we have given the new method the acronym of the OPSCF (Open Path Symmetrized Correlation Function) method. Despite it being an imaginary time formulation, it, nevertheless, explicitly accounts for phase information up to second-order in the path difference variables. The new scheme is applied to a set of one-dimensional problems that are often used to benchmark new approximate quantum dynamical methods, as exact results are available for them. It is important to note that the new formulation can be derived within a higher-order Trotter factorization, such as the Suzuki-Chin47,48 technique (see  Appendix B). In addition, higher-order terms in the path difference variables can be employed within a cumulant expansion approach, thus suggesting one possible path for systematic improvement of the approximation. The integrals over the path difference variables can still be done analytically, and the result is manifestly real and formulated as an imaginary-time sampling problem. An interesting consequence of this formulation as an imaginary-time sampling problem, as can be seen in the mathematical structure of Eq. (27), is that the behavior of a correlation function in time t can be rather easily anticipated qualitatively: One only needs to predict the rough behavior of the distribution of the endpoints of an open chain as the length of the chain grows, which necessarily occurs when t increases.

In the present paper, we employed open path-integral molecular dynamics (OPIMD) to sample the resulting expression for the approximate correlation function, however, any number of sampling methods could be applied, including Monte Carlo,43 molecular dynamics,13 hybrid Monte Carlo,13,44 or path-integral metadynamics.45,46 At the highest temperature studied (T = 1), the results are at least as accurate as the popular ring-polymer molecular dynamics (RPMD) approach and the TGA-LSC-IVR method for linear operators and show improved accuracy for nonlinear operators. At the lower temperatures studied, we start to see the breakdown of the second-order expansion. Interestingly, while RPMD and TGA-LSC-IVR also breaks down at these lower temperatures, the behavior of these methods is rather different, something that merits further investigation. At low temperatures, all methods reproduce the short-time dynamics reasonably well for linear operators. For nonlinear operators, the case of the harmonic oscillator shows that OPSCF is considerably better in the short-time limit at the low temperature studied.

We feel that the OPSCF scheme offers an interesting perspective on the formulation of quantum time correlation functions as an imaginary-time sampling problem and that it shows considerable promise. Future work on the OPSCF method will need to address several issues. First, different sampling techniques need to be explored to determine which can best treat the modified potential (r1, …, rP+1) appearing in Eq. (27). Given that this potential involves the Hessian, it is expected that either Monte Carlo or reweighted molecular dynamics will be the most efficient. However, this raises a second issue that must be addressed, namely, how best to incorporate the Hessian. Is the exact, analytical Hessian truly needed or could it be approximated in some way? Third, we will need to see how much the higher-order terms in the path difference variables and the cumulant expansion approach improve the accuracy of the results. Fourth, it could be useful to use the OPSCF results as input to a Bayesian76 or other types of machine learning schemes as part of an iterative improvement approach. Fifth, instead of performing a second-order expansion and integrating over the difference variables, it would be interesting to consider a complex Langevin equation approach77–80 applied to both r and s variables, allowing the latter to be expanded up to a few higher orders. Complex Langevin dynamics may improve the OPSCF method for position correlation functions and also allow for direct calculation of momentum correlation functions. Sixth, the OPSCF transformation could be applied to the Kubo-transformed correlation function to see if this function is more amenable to the second-order expansion in the path-difference variables. Towards the end of Sec. IV D, we discussed the modification of the OPSCF method needed for the calculation of momentum correlation functions. As noted, the equivalent of Eq. (4) for momentum operators does not result in a cyclic ring polymer but rather, an open chain, which makes the transformation of Eq. (7) and the second-order expansion less straightforward. While time derivatives of position operators can determine some momentum-dependent operators, it is necessary to devise alternative schemes, like complex Langevin dynamics, that may be beneficial in the computation of momentum correlation functions and general correlation functions for systems when the second-order expansion fails. Finally, it will be important to test the approach on a true condensed-phase system, such as liquid para-H2 or water in order to compare to other methods that have been used to obtain spectral or transport properties from correlation functions for these systems.

See supplementary material for additional figures reporting (1) convergence of the OPSCF correlation functions of a harmonic oscillator at β = 1 and β = 8 as a function of the number of path-integral beads; (2) OPSCF correlation functions of the weakly anharmonic and quartic potentials computed using the Suzuki-Chin factorization at β = 1 and β = 8; and (3) Kubo-transformed correlation functions from RPMD simulations for all four potentials studied.

The authors gratefully acknowledge Nancy Makri for informative and stimulating discussions. This work was supported by the National Science Foundation Award No. CHE-1566085 and the Margaret Strauss Kramer Fellowship.

As discussed throughout the text, cumulant expansions provide a route to the systematic improvement to the OPSCF method by rendering the integrals generated when including higher-order terms in si amenable to analytical evaluation. Depending on the order of expansion, different contributions to the method would be improved. Specifically, the positive definite distribution, Eq. (14), involves only even orders of si while the phase factor, Eq. (15), involves only odd orders of si. As an illustration of this idea, this section of the appendix will examine the first cumulant corresponding to an expansion of the potential energy up to third order in si, which occurs through the phase factor and, therefore, can only affect t ≠ 0 results. We begin with the expansion of the potential energy,

i=2PU(xi)i=P+22PU(xi)=i=2PU(ri)si+124i=2PU(ri)si3.
(A1)

This results in a modified phase factor different from Eq. (15) while the distribution ρ remains the same,

ϕ(r1rP+1,s2sP)=mPt2|τc|22i=2P1(ri+1ri)(si+1si)+2s2(r2r1)2sP(rP+1rP)tPi=2PU(ri)sit24Pi=2PU(ri)si3.
(A2)

If we assume that the third-order term is merely a correction, we can write the phase factor as

ϕ(r1rP+1,s2sP)=mPt2|τc|22i=2P1(ri+1ri)(si+1si)+2s2(r2r1)2sP(rP+1rP)tPi=2PU(ri)si+δϕ,
(A3)

where

δϕ=t24Pi=2PU(ri)si3.
(A4)

The integration of S = (s2, …, sP) can now be written as

ds2dsPe12STMS+iKSeiδϕ=eiδϕ,
(A5)

where the average is taken with respect to the second-order expansion. We approximate the average of the phase correction,

eiδϕei<δϕ>=eit24Pi=2PU(ri)si3.
(A6)

From the above equation, we can derive a correction term to the OPSCF method. Specifically, we can define vectors V3,i = U‴(ri) and S3,i=si3, where i = 2, , P such that

i=2PU(ri)si3=V3TS3,
(A7)
V3TS3=dse12STMS+iKTSV3TS3.
(A8)

In order to perform this integral, we must first redefine M = OTDO where D is a diagonal matrix, O is the orthogonal matrix which diagonalizes M, s = OTσ where, since detO = 1, ds = , and KT = λTO,

V3TS3=dσe12σTDσ+iλTσV3TOTσ3.
(A9)

Expanding the cubic function accordingly results in a sum of Gaussian integrals

V3TS3=i,j,k,lV3iOjiOkiOli×dσ2dσPe12i=2PDiiσi22iλiσiσjσkσl.
(A10)

By completing the square,

V3TS3=i,j,k,lV3iOjiOkiOlidσ2dσP×ei=2Pλi22Diie12i=2PDiiσiiλiDii2σjσkσl,
(A11)

we can define an appropriate substitution, yi=σiiλiDii, whose inverse is σi=yi+iλiDii to solve the set of integrals,

i,j,k,lV3iOjiOkiOlidy2dyPei=2Pλi22Diie12i=2PDiiyi2  ×yj+iλjDjjyk+iλkDkkyl+iλlDll.
(A12)

Expanding the terms yields

yj+iλjDjjyk+iλkDkkyl+iλlDll  =yjykylλkλlDkkDllyjλjλlDjjDllykλkλjDkkDjjyl  +iλlDllyjyk+iλkDkkyjyl+iλjDjjylykiλjλkλlDjjDkkDll.
(A13)

We recognize that the first four terms have odd exponents in at least one of the ym terms, which, when multiplied by the Gaussian, will integrate to 0. Furthermore, for the next three terms, the integral is non-zero when the indices of ym are the same, and the final term is a constant. Thus, the only remaining terms, which interestingly are all imaginary, will be

yj+iλjDjjyk+iλkDkkyl+iλlDll  =iλlDllyj2δjk+iλkDkkyl2δjl+iλjDjjyk2δkliλjλkλlDjjDkkDll.
(A14)

When the Gaussian integrals are performed, the result we obtain is

V3TS3=m=2P2πDmm1/2eλm22Dmmi,j,k,lV3iOjiOkiOli×iλlDllDjjδjk+iλkDkkDllδjl+iλjDjjDkkδkliλjλkλlDjjDkkDll  .
(A15)

If necessary, Eq. (A15) can be transformed back to depend explicitly on the original vector K(r′) and matrix M(r). Specifically,

m=2P2πDmm1/2eλm22Dmm=2πP1det(M(r))1/2e12KT(r)M1(r)K(r)
(A16)

and

k=2P3iV3kMkk1Yk̃iV3kYk̃3,
(A17)

where

Yk̃=j=2PKjMjk1.
(A18)

With this short hand, we can write the average of the correction term as

δϕ=it24P2πP1det(M(r))1/2e12KT(r)M1(r)K(r)×k=2P3Mkk1Yk̃2V3kYk̃
(A19)

from which it can be seen that exp(iδϕ⟩) is purely real. While this functional form is rather complex and would be cumbersome to implement directly into a molecular dynamics scheme, this clearly shows a way to improve the OPSCF method. In principle, higher-order cumulants can be derived as well to better improve the methodology, and the resulting path integral will remain positive definite.

The derivation of the OPSCF method using the Suzuki-Chin47,48,65,66 factorization follows many of the similar steps as in Sec. II. However, we will point out some of the differences in this appendix. Unlike the Trotter factorization, the Suzuki-Chin factorization decomposes the complex-time propagator as follows:

e2ϵĤ=eϵUe^/3eϵT^e4ϵUm^/3eϵT^eϵUe^/3,
(B1)

where ϵ = c/2 and ϵ=iτc*/2 for the forward and backward complex-time propagators, respectively, and Ue^ and Um^, in one-dimensions for a single particle, are defined by

Ue^=Û+αϵ226mÛx2,Um^=Û+(1α)ϵ2212mÛx2,
(B2)

where α is a constant in the range of [0,1] and specifically chosen to be zero in our implementation. Following closely the derivation of the path integral representation of forward and backward complex-time propagators of Krilov,76 we obtain new expressions for the phase factor and positive-definite distribution that contain first-order derivatives of the potential,

ρ(x1x2P)=mP2π|τc|Pek=12PmPβ4|τc|2xk+1xk2×eβ3Pk=1PU(x2k1)+2U(x2k)+β2212t224mP2αU(x2k1)2+(1α)U(x2k)2x2P+1=x1,
(B3)
ϕ(x1x2P)=mPt2|τc|2k=1Pxk+1xk2k=P+12Pxk+1xk2t3P2s=2P/2U(x2k1)s=P/2+2PU(x2k1)+4s=1P/2U(x2k)s=P/2+1PU(x2k)+3β224t212mP2αs=2P/2U(x2k1)2s=P/2+2PU(x2k1)2+(1α)s=1P/2U(x2k)2s=P/2+1PU(x2k)2.
(B4)

After performing the transformation described in Eq. (7), we obtain a separable distribution function, ρ(r1rP+1, s2sP) = ρ1(r1rP+1)ρ2(r2rP, s2sP),

ρ1(r1rP+1)=mP2π|τc|PemPβ4|τc|22i=1Pri+1ri2eβ3PU(r1)+U(rP+1)+k=2P/22U(r2k1)+k=1P/24U(r2k)×eβ3Pβ2212t224mP2αU(r1)2+αU(rP+1)2+k=2P/2α2U(r2k1)2+k=1P/2(1α)2U(r2k)2,
(B5)
ρ2(r2rP,s2sP)=emPβ4|τc|212i=2P1si+1si2+12(s22+sP2)eβ3Pk=2P/214U(r2k1)s2k12+k=1P/212U(r2k)s2k2×eβ3Pβ2212t224mP2k=2P/2α12s2k12U(r2k1)U(r2k1)+U(r2k1)2×eβ3Pβ2212t224mP2k=1P/2(1α)12s2k2U(r2k)U(r2k)+U(r2k)2.
(B6)

The phase factor, under the transformation, is now defined as

ϕ(r1rP+1,s2sP)=mPt2|τc|22i=2P1(ri+1ri)(si+1si)+2s2(r2r1)2sP(rP+1rP)2t3Pk=2P/2U(r2k1)s2k1+k=1P/22U(r2k)s2k+3β224t212mP2×k=2P/2αU(r2k1)U(r2k1)s2k1+k=1P/2(1α)U(r2k)U(r2k)s2k.
(B7)

With the new definition of ρ2, the matrix M with the same (P − 1) × (P − 1) dimensions is

Mij(r)=2A+β3PUM(ri)δijAδi,j+1Aδi+1,j,
(B8)

where we introduce a potential energy UM that is dependent on the parity of the bead number,

UM(r2k)=U(r2k)+β2212t224mP2*(1α)U(r2k)U(r2k)+U(r2k)2,UM(r2k1)=12U(r2k1)+β2212t224mP2*αU(r2k1)U(r2k1)+U(r2k1)2.
(B9)

A similar procedure is done for the phase factor and the vector K in which we introduce the potential UK,

Ki(r)=γ2riri1ri+1t3PUK(ri),
(B10)
UK(r2k)=4U(r2k)+3β224t26mP2×(1α)U(r2k)U(r2k),UK(r2k1)=2U(r2k1)+3β224t26mP2×αU(r2k1)U(r2k1).
(B11)

With these definitions of M(r) and K(r′), the rest of the derivation follows exactly as in Sec. II. The resulting Hamiltonian for a path integral scheme is as follows:

H̃=k=1P+1pk22m+k=1P12mωp2(rk+1rk)2+13PU(r1)+U(rP+1)+ασU(r1)2+U(rP+1)2+23Pk=2P/2U(r2k1)+σαU(r2k1)2+k=1P/22U(r2k)+σ(1α)U(r2k)2+12βṼ(r1,,rP+1).
(B12)

In theory, one could introduce a weighting function to mimic a more traditional OPSCF scheme,66 but we implemented the above Hamiltonian for this current project. As noted in Sec. IV, this higher-order factorization improves the correlation functions for the anharmonic potential but is insufficient for the quartic oscillator.

1.
R. P.
Feynman
, “
Space-time approach to non-relativistic quantum mechanics
,”
Rev. Mod. Phys.
20
,
367
387
(
1948
).
2.
R. P.
Feynamn
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill
,
New York
,
1965
).
3.
R. P.
Feynman
,
Statistical Mechanics: A Set of Lectures
(
Addison-Wesley
,
Reading, MA
,
1998
).
4.
R. P.
Feynman
,
Feynman’s Thesis: A New Approach to Quantum Theory
, edited by
L. M.
Brown
(
World Scientific
,
Singapore
,
2005
).
5.
D.
Chandler
and
P. G.
Wolynes
, “
Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids
,”
J. Chem. Phys.
74
,
4078
4095
(
1981
).
6.
J.
Cao
and
G. A.
Voth
, “
A new perspective on quantum time correlation functions
,”
J. Chem. Phys.
99
,
10070
10073
(
1993
).
7.
J.
Cao
and
G. A.
Voth
, “
The formulation of quantum statistical mechanics based on the Feynman path centroid density. II. Dynamical properties
,”
J. Chem. Phys.
100
,
5106
5117
(
1994
).
8.
J.
Cao
and
G. A.
Voth
, “
The formulation of quantum statistical mechanics based on the Feynman path centroid density. III. Phase space formalism and analysis of centroid molecular dynamics
,”
J. Chem. Phys.
101
,
6157
6167
(
1994
).
9.
J.
Cao
and
G. A.
Voth
, “
The formulation of quantum statistical mechanics based on the Feynman path centroid density. IV. Algorithms for centroid molecular dynamics
,”
J. Chem. Phys.
101
,
6168
(
1994
).
10.
I. R.
Craig
and
D. E.
Manolopoulos
, “
Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics
,”
J. Chem. Phys.
121
,
3368
3373
(
2004
).
11.
M.
Parrinello
and
A.
Rahman
, “
Study of an F center in molten KCl
,”
J. Chem. Phys.
80
,
860
867
(
1984
).
12.
R. W.
Hall
and
B. J.
Berne
, “
Nonergodicity in path integral molecular-dynamics
,”
J. Chem. Phys.
81
,
3641
(
1984
).
13.
M. E.
Tuckerman
,
B. J.
Berne
,
G. J.
Martyna
, and
M. L.
Klein
, “
Efficient molecular dynamics and hybrid Monte Carlo algorithms for path integrals
,”
J. Chem. Phys.
99
,
2796
2808
(
1993
).
14.
S.
Habershon
,
B. J.
Braams
, and
D. E.
Manolopoulos
, “
Quantum mechanical correlation functions, maximum entropy analytic continuation, and ring polymer molecular dynamics
,”
J. Chem. Phys.
127
,
174108
(
2007
).
15.
D. R.
Reichman
,
P.-N.
Roy
,
S.
Jang
, and
G. A.
Voth
, “
A Feynman path centroid dynamics approach for the computation of time correlation functions involving nonlinear operators
,”
J. Chem. Phys.
113
,
919
929
(
2000
).
16.
F.
Paesani
and
G. A.
Voth
, “
Nonlinear quantum time correlation functions from centroid molecular dynamics and the maximum entropy method
,”
J. Chem. Phys.
129
,
194113
(
2008
).
17.
A.
Horikoshi
and
K.
Kinugawa
, “
Effective potential analytic continuation approach for real time quantum correlation functions involving nonlinear operators
,”
J. Chem. Phys.
122
,
174104
(
2005
).
18.
V.
Krishna
and
G. A.
Voth
, “
Evaluation of nonlinear quantum time correlation functions within the centroid dynamics formulation
,”
J. Phys. Chem. B
110
,
18953
18957
(
2006
).
19.
I.
Georgescu
and
V. A.
Mandelshtam
, “
Molecular dynamics with quantum fluctuations
,”
Phys. Rev. B
82
,
094305
(
2010
).
20.
I.
Georgescu
,
J.
Deckman
,
L. J.
Fredrickson
, and
V. A.
Mandelshtam
, “
Thermal Gaussian molecular dynamics for quantum dynamics simulations of many-body systems: Application to liquid para-hydrogen
,”
J. Chem. Phys.
134
,
174109
(
2011
).
21.
P. A.
Frantsuzov
and
V. A.
Mandelshtam
, “
Quantum statistical mechanics with Gaussians: Equilibrium properties of van der Waals clusters
,”
J. Chem. Phys.
121
,
9247
9256
(
2004
).
22.
M. F.
Herman
and
E.
Kluk
, “
A semiclassical justification for the use of non-spreading wavepackets in dynamics calculations
,”
Chem. Phys.
91
,
27
34
(
1984
).
23.
H.
Wang
,
X.
Sun
, and
W. H.
Miller
, “
Semiclassical approximations for the calculation of thermal rate constants for chemical reactions in complex molecular systems
,”
J. Chem. Phys.
108
,
9726
9736
(
1998
).
24.
X.
Sun
,
H. B.
Wang
, and
W. H.
Miller
, “
Semiclassical theory of electronically nonadiabatic dynamics: Results of a linearized approximation to the initial value representation
,”
J. Chem. Phys.
109
,
7064
7074
(
1998
).
25.
Q.
Shi
and
E.
Geva
, “
Nonradiative electronic relaxation rate constants from approximations based on linearizing the path-integral forward-backward action
,”
J. Phys. Chem. A
108
,
6109
6116
(
2004
).
26.
S.
Bonella
,
D.
Montemayor
, and
D. F.
Coker
, “
Linearized path integral approach for calculating nonadiabatic time correlation functions
,”
Proc. Natl. Acad. Sci. U. S. A.
102
,
6715
6719
(
2005
).
27.
J.
Liu
and
W. H.
Miller
, “
Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions
,”
J. Chem. Phys.
125
,
224104
(
2006
).
28.
J.
Liu
and
W. H.
Miller
, “
Linearized semiclassical initial value time correlation functions using the thermal Gaussian approximation: Applications to condensed phase systems
,”
J. Chem. Phys.
127
,
114506
(
2007
).
29.
J.
Liu
and
W. H.
Miller
, “
Test of the consistency of various linearized semiclassical initial value time correlation functions in application to inelastic neutron scattering from liquid para-hydrogen
,”
J. Chem. Phys.
128
,
144511
(
2008
).
30.
E. R.
Dunkel
,
S.
Bonella
, and
D. F.
Coker
, “
Iterative linearized approach to nonadiabatic dynamics
,”
J. Chem. Phys.
129
,
114106
(
2008
).
31.
S.
Bonella
,
M.
Monteferrante
,
C.
Pierleoni
, and
G.
Ciccotti
, “
Path integral based calculations of symmetrized time correlation functions. I
,”
J. Chem. Phys.
133
,
164104
(
2010
).
32.
S.
Bonella
,
M.
Monteferrante
,
C.
Pierleoni
, and
G.
Ciccotti
, “
Path integral based calculations of symmetrized time correlation functions. II
,”
J. Chem. Phys.
133
,
164105
(
2010
).
33.
M.
Monteferrante
,
S.
Bonella
, and
G.
Ciccotti
, “
Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging
,”
Mol. Phys.
109
,
3015
3027
(
2011
).
34.
P. F.
Huo
and
D. F.
Coker
, “
Communication: Partial linearized density matrix dynamics for dissipative, non-adiabatic quantum evolution
,”
J. Chem. Phys.
135
,
201101
(
2011
).
35.
P. F.
Huo
and
D. F.
Coker
, “
Consistent schemes for non-adiabatic dynamics derived from partial linearized density matrix propagation
,”
J. Chem. Phys.
137
,
22A535
(
2012
).
36.
J. A.
Poulsen
,
G.
Nyman
, and
P. J.
Rossky
, “
Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method
,”
J. Chem. Phys.
119
,
12179
12193
(
2003
).
37.
V.
Jadhao
and
N.
Makri
, “
Iterative Monte Carlo for quantum dynamics
,”
J. Chem. Phys.
129
,
161102
(
2008
).
38.
V.
Jadhao
and
N.
Makri
, “
Iterative Monte Carlo with bead-adapted sampling for complex-time correlation functions
,”
J. Chem. Phys.
132
,
104110
(
2010
).
39.
V.
Jadhao
and
N.
Makri
, “
Iterative Monte Carlo path integral with optimal grids from whole-necklace sampling
,”
J. Chem. Phys.
133
,
114105
(
2010
).
40.
C. O.
Baltaretu
and
N.
Makri
, “
Iterative Monte Carlo formulation of real-time correlation functions
,”
J. Chem. Phys.
133
,
164103
(
2010
).
41.
G.
Krilov
and
B. J.
Berne
, “
Real time quantum correlation functions. I. Centroid molecular dynamics of anharmonic systems
,”
J. Chem. Phys.
111
,
9140
9146
(
1999
).
42.
M. E.
Tuckerman
,
Statistical Mechanics: Theory and Molecular Simulation
(
Oxford University Press
,
Oxford, UK
,
2010
).
43.
E. L.
Pollock
and
D. M.
Ceperley
, “
Simulation of quantum many-body systems by path-integral methods
,”
Phys. Rev. B
30
,
2555
2568
(
1984
).
44.
S.
Duane
,
A. D.
Kennedy
,
B. J.
Pendleton
, and
D.
Roweth
, “
Hybrid Monte Carlo
,”
Phys. Lett. B
195
,
216
(
1987
).
45.
R. G.
Quhe
,
M.
Nava
,
P.
Tiwary
, and
M.
Parrinello
, “
Path integral metadynamics
,”
J. Chem. Theory Comput.
11
,
1383
1388
(
2015
).
46.
M.
Nava
,
R. G.
Quhe
,
F.
Palazzesi
,
P.
Tiwary
, and
M.
Parrinello
, “
de Broglie swapping metadynamics for quantum and classical sampling
,”
J. Chem. Theory Comput.
11
,
5114
5119
(
2015
).
47.
M.
Suzuki
, “
Hybrid exponential product formulas for unbounded operators with possible applications to Monte Carlo simulations
,”
Phys. Lett. A
201
,
425
428
(
1995
).
48.
S. A.
Chin
, “
Symplectic integrators from composite operator factorizations
,”
Phys. Lett. A
226
,
344
348
(
1997
).
49.
Q.
Shi
and
E.
Geva
, “
Semiclassical theory of vibrational energy relaxation in the condensed phase
,”
J. Phys. Chem. A
107
,
9059
9069
(
2003
).
50.
J.
Liu
and
W. H.
Miller
, “
A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids
,”
J. Chem. Phys.
131
,
074113
(
2009
).
51.
D. M.
Ceperley
and
E. L.
Pollock
, “
Path-integral computation of the low-temperature properties of liquid 4He
,”
Phys. Rev. Lett.
56
,
351
354
(
1986
).
52.
D. M.
Ceperley
, “
Path integrals in the theory of condensed helium
,”
Rev. Mod. Phys.
67
,
279
355
(
1995
).
53.
J. A.
Morrone
,
V.
Srinivasan
,
D.
Sebastiani
, and
R.
Car
, “
Proton momentum distribution in water: An open path integral molecular dynamics study
,”
J. Chem. Phys.
126
,
234504
(
2007
).
54.
M.
Ceriotti
,
M.
Parrinello
,
T. E.
Markland
, and
D. E.
Manolopoulos
, “
Efficient stochastic thermostatting of path integral molecular dynamics
,”
J. Chem. Phys.
133
,
124104
(
2010
).
55.
M.
Ceriotti
,
D. E.
Manolopoulos
, and
M.
Parrinello
, “
Accelerating the convergence of path integral dynamics with a generalized Langevin equation
,”
J. Chem. Phys.
134
,
084104
(
2011
).
56.
G.
Ciccotti
and
S.
Meloni
, “
Temperature accelerated Monte Carlo (TAMC): A method for sampling the free energy surface of non-analytical collective variables
,”
Phys. Chem. Chem. Phys.
13
,
5952
5959
(
2011
).
57.
A.
Pérez
,
M. E.
Tuckerman
, and
M. H.
Müser
, “
A comparative study of the centroid and ring-polymer molecular dynamics methods for approximating quantum time correlation functions from path integrals
,”
J. Chem. Phys.
130
,
184105
(
2009
).
58.
S.
Jang
and
G. A.
Voth
, “
A derivation of centroid molecular dynamics and other approximate time evolution methods for path integral centroid variables
,”
J. Chem. Phys.
111
,
2371
2384
(
1999
).
59.
J.
Liu
and
W. H.
Miller
, “
Real time correlation function in a single phase space integral beyond the linearized semiclassical initial value representation
,”
J. Chem. Phys.
126
,
234110
(
2007
).
60.
T. D.
Hone
,
P. J.
Rossky
, and
G. A.
Voth
, “
A comparative study of imaginary time path integral based methods for quantum dynamics
,”
J. Chem. Phys.
124
,
154103
(
2006
).
61.
B. J.
Ka
and
G. A.
Voth
, “
An efficient and accurate implementation of centroid molecular dynamics using a Gaussian approximation
,”
J. Phys. Chem. A
109
,
11609
11617
(
2005
).
62.
J.
Liu
and
W. H.
Miller
, “
An approach for generating trajectory-based dynamics which conserves the canonical distribution in the phase space formulation of quantum mechanics. I. Theories
,”
J. Chem. Phys.
134
,
104101
(
2011
).
63.
J.
Liu
, “
Two more approaches for generating trajectory-based dynamics which conserves the canonical distribution in the phase space formulation of quantum mechanics
,”
J. Chem. Phys.
134
,
194110
(
2011
).
64.
B.
Leimkuhler
,
D. T.
Margul
, and
M. E.
Tuckerman
, “
Stochastic, resonance-free multiple time-step algorithm for molecular dynamics with very large time steps
,”
Mol. Phys.
111
,
3579
3594
(
2013
).
65.
S.
Jang
,
S.
Jang
, and
G. A.
Voth
, “
Applications of higher order composite factorization schemes in imaginary time path integral simulations
,”
J. Chem. Phys.
115
,
7832
7842
(
2001
).
66.
A.
Pérez
and
M. E.
Tuckerman
, “
Improving the convergence of closed and open path integral molecular dynamics via higher order Trotter factorization schemes
,”
J. Chem. Phys.
135
,
064104
(
2011
).
67.
Mathematica, Version 11.2,
Wolfram Research, Inc.
,
2017
.
68.
T. J. H.
Hele
,
M. J.
Willatt
,
A.
Muolo
, and
S. C.
Althorpe
, “
Boltzmann-conserving classical dynamics in quantum time-correlation functions: ‘Matsubara dynamics
,’”
J. Chem. Phys.
142
,
134103
(
2015
).
69.
T. J. H.
Hele
,
M. J.
Willatt
,
A.
Muolo
, and
S. C.
Althorpe
, “
Communication: Relation of centroid molecular dynamics and ring-polymer molecular dynamics to exact quantum dynamics
,”
J. Chem. Phys.
142
,
191101
(
2015
).
70.
G. A.
Voth
, “
Feynman path integral formulation of quantum mechanical transition-state theory
,”
J. Phys. Chem.
97
,
8365
8377
(
1993
).
71.
I. R.
Craig
and
D. E.
Manolopoulos
, “
Chemical reaction rates from ring polymer molecular dynamics
,”
J. Chem. Phys.
122
,
084106
(
2005
).
72.
I. R.
Craig
and
D. E.
Manolopoulos
, “
A refined ring polymer molecular dynamics theory of chemical reaction rates
,”
J. Chem. Phys.
123
,
034102
(
2005
).
73.
E.
Geva
,
Q.
Shi
, and
G. A.
Voth
, “
Quantum-mechanical reaction rate constants from centroid molecular dynamics simulations
,”
J. Chem. Phys.
115
,
9209
9222
(
2001
).
74.
C.
Venkataraman
and
W. H.
Miller
, “
Chemical reaction rates using the semiclassical Van Vleck initial value representation
,”
J. Chem. Phys.
126
,
094104
(
2007
).
75.
W. H.
Miller
,
S. D.
Schwartz
, and
J. W.
Tromp
, “
Quantum mechanical rate constants for bimolecular reactions
,”
J. Chem. Phys.
79
,
4889
4898
(
1983
).
76.
G.
Krilov
,
E.
Sim
, and
B. J.
Berne
, “
On the Bayesian approach to calculating time correlation functions in quantum systems; reaction dynamics and spectroscopy
,”
Chem. Phys.
268
,
21
34
(
2001
).
77.
G.
Parisi
, “
On complex probabilities
,”
Phys. Lett. B
131
,
393
395
(
1983
).
78.
J. R.
Klauder
, “
A Langevin approach to fermion and quantum spin correlation functions
,”
J. Phys. A: Math. Gen.
16
,
L317
(
1983
).
79.
L. L.
Salcedo
, “
Representation of complex probabilities
,”
J. Math. Phys.
38
,
1710
1722
(
1997
).
80.
D.
Weingarten
, “
Complex probabilities on RN as real probabilities on CN and an application to path integrals
,”
Phys. Rev. Lett.
89
,
240201
(
2002
).
81.
M. E.
Tuckerman
,
G. J.
Martyna
, and
B. J.
Berne
, “
Reversible multiple time-scale molecular dynamics
,”
J. Chem. Phys.
97
,
1990
(
1992
).

Supplementary Material