We introduce a new semiclassical (SC) framework, the Mixed Quantum-Classical Initial Value Representation (MQC-IVR), that can be tuned to reproduce existing quantum-limit and classical-limit SC approximations to quantum real-time correlation functions. Applying a modified Filinov transformation to a quantum-limit SC formulation leads to the association of a Filinov parameter with each degree of freedom in the system; varying this parameter from zero to infinity controls the extent of quantization of the corresponding mode. The resulting MQC-IVR expression provides a consistent dynamic framework for mixed quantum-classical simulations and we demonstrate its numerical accuracy in the calculation of real-time correlation functions for a model 1D system and a model 2D system over the full range of quantum- to classical-limit behaviors.

Understanding and characterizing the role of quantum effects in complex chemical and biological systems is an important and ongoing challenge in several areas of research including the rational design of materials for renewable energy harvesting. Theoretical simulation methods are uniquely suited to uncover reaction mechanisms at an atomistic level; however, the exponential costs of exact quantum mechanical simulations for large systems make the development of approximate methods essential.

Mixed quantum-classical (MQC) methods attempt to incorporate quantum effects by treating a subset of the system degrees of freedom (DOFs) quantum mechanically while using classical molecular dynamics (MD) to simulate the remaining degrees of freedom.1–6 However, it has been shown that such methods introduce uncontrolled approximations in the feedback forces between modes described at two different levels of theory.7 More rigorous multi-physics approaches include mixed quantum-classical Liouville (MQCL) methods8–13 and mixed quantum-classical Bohmian dynamics14–16 that can, in principle, incorporate zero-point energy, tunneling, and quantum coherence effects. In particular, MQCL methods have been employed in high dimensional simulations by partitioning the system according to relative atomic masses. An alternate approach is to employ a consistent dynamic framework for all system DOFs. In this class of methods, path-integral based methods such as Ring Polymer Molecular Dynamics (RPMD)17–20 and centroid molecular dynamics (CMD)21,22 have had a great deal of success in describing condensed-phase charge transfer reactions where quantum coherence effects do not play a significant role. Yet other path-integral based methods such as the mapping-variable RPMD23,24 and partially linearized density matrix methods25–27 are able to capture some quantum coherences in simulations of nonadiabatic dynamics.

In this paper, we derive a new method based on the semiclassical Initial Value Representation (SC-IVR) that employs a consistent dynamic framework to implement the MQC idea: a few system DOFs are treated in the quantum-limit of the SC theory while the remaining modes are described using classical-limit SC methods. Classical-limit SC methods for real-time correlation functions like the Linearized SC (LSC)-IVR28–30 can describe quantum tunneling and zero-point energy, while the more accurate quantum-limit SC methods like the coherent state Double Herman-Kluk (DHK)-IVR30–32 also capture quantum coherence effects. Here, we employ a modified Filinov filtration technique33–35 to smooth the oscillatory integrand in the DHK-IVR formulation and obtain a classical-limit SC formulation much like the LSC-IVR. This work closely follows an earlier derivation of the Generalized Forward-Backward IVR (GFB-IVR),30 where modified Filinov filtration was used to smooth the oscillatory integrand in the DHK-IVR formulation and obtain a more computationally tractable quantum-limit formulation, the Forward-Backward IVR expression. In our MQC-IVR formulation, the extent to which a given mode contributes to the oscillatory integrand is controlled by the value of its associated Filinov parameter. This feature allows us to implement a quantum-classical partitioning of the system without introducing boundary problems associated with existing multi-physics approaches.

This paper is organized as follows: in Sec. II, we present a derivation of the MQC-IVR formulation and show that for the two limiting values of the Filinov parameter, 0 and ∞, MQC-IVR reduces to existing quantum- and classical-limit SC-IVR methods, respectively. In Sec. III, we numerically demonstrate the accuracy of MQC-IVR in calculations of real-time correlation functions for a 1D model system and a 2D model system over a full range of behaviors from the quantum to the classical limit. We summarize our results in Sec. IV.

Dynamic phenomena in complex systems can be described in terms of real-time correlation functions. In quantum mechanics, a general time correlation function of two operators A ˆ and B ˆ is written as

C AB ( t ) = Tr A ˆ e i ħ H ˆ t B ˆ e i ħ H ˆ t ,
(1)

where H ˆ is the Hamiltonian of the system under consideration. In this section, we provide a derivation of the mixed quantum-classical IVR approximation to Eq. (1) and discuss its behavior in limiting cases.

We start with a brief review of the modified Filinov transformation34,35 that has been used before in a variety of SC-IVR calculations to smooth oscillatory integrands and improve Monte Carlo statistics.30,36,37 Consider the integral

I = d z g ( z ) e i f ( z ) ,
(2)

where z is an N-dimensional vector and g and f are, in general, complex functions. The phase factor eif(z) makes the integrand oscillatory and difficult to evaluate using standard Monte Carlo methods. The modified Filinov approximation to the integral in Eq. (2) is given by

I ( c ) = d z g ( z ) e i f ( z ) F ( z ; c ) ,
(3)

where F(z; c) is the smoothing factor,

F ( z ; c ) = det I + i c 2 f z 2 1 2 exp 1 2 f z c f z ,
(4)

and c is a diagonal matrix of the Filinov parameters with elements that take values in the interval [0, ∞). When c = 0, the integral in Eq. (3) reduces to the original expression, Eq. (2). In the limit c → ∞, the modified integral I(c) gives the stationary phase approximation to Eq. (2),

lim c I ( c ) = j g ( z j ) e i f ( z j ) det 1 2 π i 2 f z 2 z = z j 1 2 ,
(5)

where the sum runs over all points of stationary phase. Intermediate values of the Filinov parameters in c lead to a smoothing of the original integrand and, therefore, better Monte Carlo statistics.

We apply the modified Filinov transformation34,35 to the DHK-IVR expression for a real-time correlation function. Our derivation follows closely that of the GFB-IVR method.30 Consider the following position-space matrix element:

q f | U ˆ | q i = q f | e i ħ H ˆ t B ˆ e i ħ H ˆ t | q i .
(6)

The SC version of Eq. (6) can be obtained using the generalized Herman–Kluk IVR approximation31 for the time-evolution operator of an N-dimensional system, where the initial and final coherent states have widths γ0 and γt, respectively,

e i ħ H ˆ t = 1 ( 2 π ħ ) N d q 0 d p 0 C t ( p 0 , q 0 , γ 0 , γ t ) × e i ħ S t ( p 0 , q 0 ) | p t q t γ t p 0 q 0 γ 0 | .
(7)

Here, (p0, q0) and (pt, qt) represent the position and momentum corresponding to the initial and final coherent states of width γ0 and γt, respectively, St(p0, q0) is the classical action along the trajectory, and the Herman–Kluk prefactor is given by

C t ( p 0 , q 0 , γ 0 , γ t ) = det 1 2 γ t 1 2 M q q γ 0 1 2 + γ t 1 2 M p p γ 0 1 2 i ħ γ t 1 2 M q p γ 0 1 2 + i ħ γ t 1 2 M p q γ 0 1 2 1 2 ,
(8)

with M α β = α t β 0 being elements of the monodromy matrix. The coherent state wavefunctions in Eq. (7) are

x | p q γ = det ( γ ) π N 1 4 × exp 1 2 ( x q ) γ ( x q ) + i ħ p ( x q ) .
(9)

In the rest of the paper, we do not include the γ symbol in the coherent state labels to simplify notation. Substituting Eq. (7) for the two propagators in Eq. (6), we obtain the DHK-IVR expression for the position matrix element,

q f | U ˆ | q i = 1 ( 2 π ħ ) 2 N d q 0 d p 0 d q t d p t C t ( p 0 , q 0 ) C t ( p t , q t ) e i ħ [ S t ( p 0 , q 0 ) + S t ( p t , q t ) ] × p 0 q 0 | q i q f | p 0 q 0 p t q t | B ˆ | p t q t .
(10)

The complex phase of the integrand in Eq. (10) is written as

f ( z ) = S t ( p 0 , q 0 ) + S t ( p t , q t ) + i 2 ( q 0 q i ) γ 0 ( q 0 q i ) + p 0 ( q 0 q i ) + i 2 ( q 0 q f ) γ 0 ( q 0 q f ) p 0 ( q 0 q f ) + i 4 Δ q γ t Δ q + i 4 Δ p γ t 1 Δ p + 1 2 ( p t + p t ) Δ q ,
(11)

where the collective vector z = ( q t , p t , q 0 , p 0 ) , Δ x = x t x t , x (p,q) are changes in position and momentum or “jumps” in the classical trajectory at time t due to the action of operator B ˆ , and we set ħ = 1. We note that Eq. (11) excludes any contribution from the operator B ˆ to the complex phase of the integrand in Eq. (10).

Applying the modified Filinov transform to Eq. (10), as shown in Eqs. (3) and (4), requires evaluation of the first and second derivatives of the phase term in Eq. (11). Following the GFB-IVR derivation,30 we factorize the first derivative of Eq. (11) with respect to z and express it in terms of a matrix K and vector y,

K y = M p q b T + i M q q b T γ 0 0 i 2 γ t 1 2 I M p p b T + i M q p b T γ 0 0 1 2 I i 2 γ t 1 0 i γ 0 1 2 M p q f T i 2 M q q f T γ t 1 2 M q q f T i 2 M p q f T γ t 1 0 I 1 2 M p p f T i 2 M q p f T γ t 1 2 M q p f T i 2 M p p f T γ t 1 q 0 q f q 0 q i Δ q Δ p ,
(12)

where the superscript f/b indicates components of the monodromy matrix along the forward/backward trajectory and I is the N-dimensional identity matrix. The second derivative of f(z) is written as

2 f ( z ) z 2 = z Ky = K y z = KJ = K M q q b M q p b 0 0 0 0 I 0 I 0 M q q f M q p f 0 I M p q f M p p f ,
(13)

where, as is typical in modified Filinov filtration of semiclassical integrands,30 dependence of the monodromy matrix on the integration variable is ignored. Substituting Eqs. (12) and (13) into Eq. (4) and defining the new matrix of Filinov parameters as c ̃ = K T cK yields

F ( z ; c ̃ ) = det ( K T + i c ̃ T J ) 1 2 det ( K T ) 1 2 exp 1 2 y c ̃ y .
(14)

The vector y defined in Eq. (12) has components dependent on positions of the initial (qi) and final (qf) states. In order to remove this dependence, the Filinov parameter matrix c is chosen to have the form

c ̃ = 0 0 0 0 0 0 0 0 0 0 c q 0 0 0 0 c p ,
(15)

where cq and cp are N × N diagonal matrices, and the exponential factor in Eq. (14) becomes

exp 1 2 y c ̃ y = exp 1 2 Δ q c q Δ q exp 1 2 Δ p c p Δ p .
(16)

Combining the results of Eqs. (10), (14), and (16) gives the mixed quantum-classical IVR expression for the matrix element of operator U ˆ ,

q f | U ˆ | q i = 1 ( 2 π ) 2 N d q 0 d p 0 d Δ q d Δ p q f | p 0 q 0 p t q t | B ˆ | p t q t e i [ S t ( p 0 , q 0 ) + S t ( p t , q t ) ] × C t ( p 0 , q 0 ) C t ( p t , q t ) det ( K T + i c T J ) 1 2 det ( K T ) 1 2 e 1 2 Δ q c q Δ q e 1 2 Δ p c p Δ p p 0 q 0 | q i .
(17)

Substituting the MQC-IVR approximation to U ˆ from Eq. (17) in Eq. (1) yields the MQC-IVR expression for the real-time correlation function,

C AB ( t ) = 1 ( 2 π ) 2 N d q 0 d p 0 d Δ q d Δ p p 0 q 0 | A ˆ | p 0 q 0 p t q t | B ˆ | p t q t × e i [ S t ( p 0 , q 0 ) + S t ( p t , q t ) ] D ( p 0 , q 0 , Δ q , Δ p , c q , c p ) e 1 2 Δ q c q Δ q e 1 2 Δ p c p Δ p ,
(18)

and the MQC-IVR prefactor, D(p0, q0, Δq, Δp, cq, cp), for a general operator B ˆ is provided in Appendix  A along with a detailed derivation.

For simplicity, here we consider a position-space operator B ˆ B ( q ˆ ) . In this case, the matrix elements of operator B ( q ˆ ) can be evaluated by setting the width of coherent states at time t, γt → ∞, such that

lim γ t det ( γ t ) 1 2 p t q t | B ( q ˆ ) | p t q t = ( 4 π ) N 2 δ ( Δ q ) B q t + q t 2 .
(19)

This idea was previously used in deriving the Exact Forward-Backward (EFB-IVR) formulation from the DHK-IVR and introduces no additional approximations.38,39 Inserting Eq. (19) into Eq. (18) and evaluating the integral over Δq yields the following expression for the real-time correlation function:

C AB ( t ) = 2 2 N π 3 N 2 d q 0 d p 0 d Δ p p 0 q 0 | A ˆ | p 0 q 0 B ( q t ) e i [ S t ( p 0 , q 0 ) + S t ( p t , q t ) ] D q ( p 0 , q 0 , Δ p , c p ) e 1 2 Δ p c p Δ p ,
(20)

where the prefactor is given by

D q = det γ 0 1 c p M p p b i γ 0 M q p b M p p f γ 0 + i M p q f + γ 0 M q q b + i M p q b M q q f i M q p f γ 0 + M p p b i γ 0 M q p b c p 1 M q q f i M q p f γ 0 1 2 .
(21)

The MQC-IVR expression for a momentum-space operator B ˆ B ( p ˆ ) is also provided in Appendix  B.

In Sec. II B, we derived the MQC-IVR approximation to the time-correlation function that includes a matrix of Filinov parameters, c, with values in the interval [0, ∞). Specifically, the Filinov parameter appears in gaussian exponentials and the prefactor of the MQC-IVR formulation (Eq. (18)). In the limit c → 0, the gaussian exponential factors in Eq. (18) vanish, the prefactor becomes

D 0 = C t ( p 0 , q 0 ) C t ( p t , q t ) ,
(22)

and the resulting expression for the correlation function is identical to the DHK-IVR formulation. The MQC-IVR method, thus, approaches the full semiclassical or quantum-limit IVR description of dynamics as c → 0.

For the other limiting case, c → ∞, we have

lim c e 1 2 Δ q c q Δ q e 1 2 Δ p c p Δ p D ( p 0 , q 0 , Δ q , Δ p , c q , c p ) = ( 2 π ) N δ ( Δ q ) δ ( Δ p ) D ,
(23)

where

D = lim c det ( c q c p ) 1 2 D q ( p 0 , q 0 , Δ q , Δ p , c q , c p ) = det ( 2 γ 0 ) 1 γ 0 M q q full + M p p full γ 0 i γ 0 M q p full γ 0 + i M p q full 1 2 .
(24)

Here, M α β full = α 0 β 0 are components of the monodromy matrix calculated over the full (forward and backward) trajectory. Inserting Eq. (23) in Eq. (18), and evaluating integrals over Δq, Δp analytically leads to the condition q t = q t , p t = p t indicating that the trajectory does not undergo any change in phase space coordinates at time t. Consequently, the forward and backward trajectories coincide and elements of the full monodromy matrix become M α β full = δ α β I in Eq. (24), reducing the prefactor to D = 1. The final expression for the MQC-IVR correlation function in the limit c → ∞ is

C AB ( t ) = 1 ( 2 π ) N d q 0 d p 0 × p 0 q 0 | A ˆ | p 0 q 0 p t q t | B ˆ | p t q t .
(25)

The correlation function in Eq. (25) is an alternate formulation to the LSC-IVR where Husimi functions are employed in the place of Wigner functions.40,41 Therefore, in the c → ∞ limit, the MQC-IVR method approaches a classical-limit SC-IVR formulation. This feature of the MQC-IVR method arises from treating the complex phase of the integrand in Eq. (10) independent of operator B ˆ , distinguishing it from the GFB-IVR formulation previously derived by Thoss, Wang, and Miller.30 

Characterizing the limiting behavior serves to demonstrate that the MQC-IVR method tunes between classical-limit and quantum-limit SC-IVR formulations continuously. Further, in a system with many DOFs, different values of the Filinov parameter can be associated with different modes effectively using different levels of theory to describe each degree of freedom. MQC-IVR thus provides a consistent semiclassical framework for mixed quantum-classical simulations of the dynamics in complex systems.

We demonstrate the numerical accuracy of MQC-IVR in calculating the average position as a function of time for a 1D anharmonic potential,

V ( x ) = 1 2 m x ω x 2 x 2 0 . 1 x 3 + 0 . 1 x 4 ,
(26)

where mx = 1 a.u and ω x = 2 a.u . This potential has been used before to quantify the extent to which various SC-IVR based methods incorporate quantum effects.30 Substituting A ˆ = | ψ i ψ i | and B ˆ = x ˆ in Eq. (20), we obtain

x ( t ) = 2 2 π 3 2 d q 0 d p 0 d Δ p p 0 q 0 | ψ i ψ i | p 0 q 0 × e i [ S t + S t ] x t D q ( p 0 , q 0 , Δ p , c p ) e 1 2 c p Δ p 2 ,
(27)

where the initial state wavefunction is chosen to be a coherent state, |ψi〉 ≡ |piqi〉, with parameters qi = 1, pi = 0, and γ = mxωx. Numerical integration of Eq. (27) is performed using a standard Monte Carlo simulation, where the initial phase space points are sampled according to

ω ( q 0 , p 0 , Δ p ) = | p 0 q 0 | ψ i | 2 e 1 2 c p Δ p 2 .
(28)

Semiclassical trajectories are time evolved using a fourth-order symplectic integrator,42 with a time step of Δt = 0.05 a.u. Conservation of classical energy along the trajectory is monitored using 1 E ( t ) / E ( 0 ) < ε , where the tolerance parameter is set to ε = 10−4, and the correct branch of the complex square root in the MQC-IVR prefactor is chosen by tracking the Maslov index.

Figure 1 compares MQC-IVR correlation functions for two limiting values of the Filinov parameter, cp, against the results of exact quantum mechanical and classical MD simulations. Exact quantum results are obtained by constructing a matrix representation of Eq. (1) in the eigenbasis of the Hamiltonian operator. Eigenfunctions and eigenvalues of the Hamiltonian are determined using the Discrete Variable Representation (DVR) method with a uniform grid.43 The quantum results show a recurrence in oscillations of 〈x(t)〉 at t > 40 a.u., unlike the classical simulation, indicating the importance of quantum effects in this system. The MQC-IVR correlation function with cp = 0.05 captures long-time oscillations with the correct frequency and a slightly reduced amplitude. It has been shown previously that the full recurrence amplitude is reproduced in the cp = 0 limit.30 For a large value of the Filinov parameter, cp = 100, the MQC-IVR result is indistinguishable from a classical simulation, exhibiting damping in the oscillatory structure and establishing the importance of quantum interference effects. In Fig. 2, we gradually change cp values from 0.05 a.u. to 100 a.u. and demonstrate that the MQC-IVR correlation function continuously tunes between quantum- and classical-limit semiclassical dynamics. Gradually increasing the Filinov parameter leads to smaller amplitude recurrences in the average position with time, and for cp > 10, the MQC-IVR results reproduce classical-limit dynamics. The number of semiclassical trajectories required to converge MQC simulations decreases with increasing cp values as reported in Table I.

FIG. 1.

The average position of the 1D anharmonic oscillator in Eq. (26) as a function of time calculated with MQC-IVR for two values of the Filinov parameter is shown in red (cp = 0.05) and green (cp = 100). Results of quantum mechanical (black) and classical (blue) simulations are shown for comparison.

FIG. 1.

The average position of the 1D anharmonic oscillator in Eq. (26) as a function of time calculated with MQC-IVR for two values of the Filinov parameter is shown in red (cp = 0.05) and green (cp = 100). Results of quantum mechanical (black) and classical (blue) simulations are shown for comparison.

Close modal
FIG. 2.

Dynamics of the average position of the anharmonic oscillator in Eq. (26) calculated for a range of Filinov parameters. Shown are the results obtained with cp = 0.05 (red), cp = 1 (green), cp = 5 (blue), cp = 10 (magenta), cp = 100 (cyan).

FIG. 2.

Dynamics of the average position of the anharmonic oscillator in Eq. (26) calculated for a range of Filinov parameters. Shown are the results obtained with cp = 0.05 (red), cp = 1 (green), cp = 5 (blue), cp = 10 (magenta), cp = 100 (cyan).

Close modal
TABLE I.

Number of semiclassical trajectories, Ntraj, and maximum absolute statistical error, ϵmax, in the MQC-IVR calculations of 〈x(t)〉 for different values of the Filinov parameter.

cp Ntraj ϵmax (a.u./102)
0.05  4.0 × 106  4.4 
1.0  2.5 × 105  4.0 
5.0  3.0 × 105  3.5 
10.0  3.5 × 104  4.8 
100.0  1.0 × 103  4.6 
cp Ntraj ϵmax (a.u./102)
0.05  4.0 × 106  4.4 
1.0  2.5 × 105  4.0 
5.0  3.0 × 105  3.5 
10.0  3.5 × 104  4.8 
100.0  1.0 × 103  4.6 

In Sec. III A, we demonstrated the ability of the MQC-IVR method to continuously tune between quantum- and classical-limit semiclassical dynamics. In this section, we employ a simple 2D model system to numerically illustrate the features that make MQC-IVR a promising method for large-scale atomistic simulations. In complex systems, coupling between many DOFs typically leads to a damping of quantum coherences. This allows us to use larger cp(>0) values and hence significantly fewer SC trajectories to calculate correlation functions in the quantum-limit leading to a reduction in computational cost.

We construct a 2D model where the 1D anharmonic oscillator used previously is coupled to a relatively high-mass harmonic mode,

V ( x , y ) = 1 2 m x ω x 2 x 2 0 . 1 x 3 + 0 . 1 x 4 + 1 2 m y ω y 2 y 2 + k x y ,
(29)

where mass my = 25 a.u., frequency ωy = 0.333 a.u., and the bilinear coupling parameter k = 2 a.u. The initial state wavefunction, |ψi〉, is the product of two coherent states,

x y | ψ i = γ x γ y π 2 1 4 e γ x 2 ( x q x ) 2 + i p x ( x q x ) × e γ y 2 ( y q y ) 2 + i p y ( y q y ) ,

where, in atomic units, qx = qy = 1, px = py = 0, γ x = 2 , and γy = 8.333.

We calculate the average position as a function of time for the x-mode. In Fig. 3, we present a set of MQC-IVR results for limiting values of the Filinov parameter for the x-mode, cx = 0.05 and cx = 108, coupled to two different values for the y-mode, cy = 1 and cy = 108. The decrease in oscillatory behavior of the exact quantum correlation function at long times compared to the isolated anharmonic mode demonstrates the reduction in quantum interference effects due to the coupling with a higher mass mode. As before, the classical result exhibits no long time recurrences and the MQC-IVR simulation with cx = cy = 108 is in excellent agreement with the classical result in Fig. 3(a).

FIG. 3.

The average position as a function of time for the x-mode of 2D potential (29). Panel (a) shows MQC-IVR results calculated with cx = cy = 108 (red) compared to a classical simulation (black). Panel (b) shows comparison between MQC-IVR results obtained with cx = 0.05, cy = 1 (red) and cx = 0.05, cy = 108 (green) and the exact quantum calculation (black).

FIG. 3.

The average position as a function of time for the x-mode of 2D potential (29). Panel (a) shows MQC-IVR results calculated with cx = cy = 108 (red) compared to a classical simulation (black). Panel (b) shows comparison between MQC-IVR results obtained with cx = 0.05, cy = 1 (red) and cx = 0.05, cy = 108 (green) and the exact quantum calculation (black).

Close modal

From Fig. 3(b), it is evident that both MQC-IVR simulations with cx = 0.05, cy = 1 and cx = 0.05, cy = 108 correctly describe the long time oscillatory behavior of the anharmonic mode. Therefore, it is sufficient to set the Filinov parameter only for the light mode to a small value in order to numerically obtain the quantum-limit correlation function. Using large values of the Filinov parameter for the y-mode makes the integrand less oscillatory and therefore less computationally expensive.

In Fig. 4, we show the reduction in statistical error for the two MQC-IVR calculations presented in Fig. 3(b). In both cases, the error grows with propagation time, but the growth rate as well as the absolute error is substantially smaller for cy = 108 than for cy = 1. Moreover, since the coupling between DOFs reduces the role of quantum interference effects, we can employ larger cx values than in simulations of the isolated anharmonic mode without loss of accuracy. The results of MQC-IVR simulations with intermediate values for cx are shown in Fig. 5; it is evident that all values of cx ≤ 1 capture quantum recurrence to a similar extent, while cx = 5 begins to resemble the classical result.

FIG. 4.

The absolute statistical error in 〈x(t)〉 from MQC-IVR calculations performed with cx = 0.05, cy = 1 (red) and cx = 0.05, cy = 108 (green) using the same number of semiclassical trajectories Ntraj = 2 × 106. The error is calculated by averaging over 5 blocks of data.

FIG. 4.

The absolute statistical error in 〈x(t)〉 from MQC-IVR calculations performed with cx = 0.05, cy = 1 (red) and cx = 0.05, cy = 108 (green) using the same number of semiclassical trajectories Ntraj = 2 × 106. The error is calculated by averaging over 5 blocks of data.

Close modal
FIG. 5.

Dynamics of the average position for the x-mode of the 2D system, Eq. (29), calculated with MQC-IVR for a range of Filinov parameters for the x-mode. Shown are the results obtained with cx = 0.05 (red), cx = 0.5 (green), cx = 1 (blue), and cx = 5 (gold). In all calculations, cy = 108. The inset focuses on oscillations at long times. For visual clarity, statistical error bars are not shown in the inset.

FIG. 5.

Dynamics of the average position for the x-mode of the 2D system, Eq. (29), calculated with MQC-IVR for a range of Filinov parameters for the x-mode. Shown are the results obtained with cx = 0.05 (red), cx = 0.5 (green), cx = 1 (blue), and cx = 5 (gold). In all calculations, cy = 108. The inset focuses on oscillations at long times. For visual clarity, statistical error bars are not shown in the inset.

Close modal

Table II reports the number of SC trajectories required to converge the MQC-IVR results presented in Fig. 5. It is clear that increasing the Filinov parameter substantially improves Monte Carlo statistics, effectively decreasing the computational cost of the simulation. We note that a time step of Δt = 0.025 a.u. is required for the cx = 0.05 simulation whereas Δt = 0.05 a.u. was found to be sufficient for other cx values. In more complex systems where the quantum mode is coupled to several more classical DOFs, we expect the role of quantum interference effects to be further reduced allowing for use of larger values of Filinov parameters and, consequently, fewer semiclassical trajectories without loss of accuracy. The MQC-IVR formulation is, thus, a promising approach for simulating the dynamics of complex systems.

TABLE II.

Number of semiclassical trajectories Ntraj required to converge the MQC-IVR calculations of 〈x(t)〉 for different values of the Filinov parameter, cx. In all calculations, cy = 108 and the maximum absolute statistical error ϵmax < 5 × 10−2 a.u.

cx Ntraj
0.05  2.0 × 106 
0.5  5.0 × 105 
1.0  4.5 × 105 
5.0  3.0 × 105 
108  4.0 × 103 
cx Ntraj
0.05  2.0 × 106 
0.5  5.0 × 105 
1.0  4.5 × 105 
5.0  3.0 × 105 
108  4.0 × 103 

We also calculate the average position as a function of time for the heavier harmonic mode. Similar to the previous case, the MQC-IVR results obtained with cy = 0.05, cx = 108 (the y-mode is in the quantum limit and the x-mode is classical) and with cx = cy = 108 are in good agreement with the exact quantum and classical results, respectively (Fig. 6). Interestingly, setting the Filinov parameter for the heavier y-mode to a large value, cy = 108, while treating the x-mode in the quantum limit, cx = 1, provides a good approximation to the quantum-limit correlation function with only minor differences in the amplitude of oscillations. We present convergence numbers for these simulations in Table III, and note that, as expected, far fewer trajectories are needed in general to converge our calculations for the heavier y-mode correlation function. In addition, we are able to treat the lighter x-mode in the classical limit and still obtain numerically accurate results for the quantum-limit y-mode correlation function.

FIG. 6.

The average position as a function of time for the y-mode of the 2D system in Eq. (29). Shown are the MQC-IVR results obtained with cy = 0.05, cx = 108 (red); cy = 108, cx = 1 (green); and cy = cx = 108 (blue). Quantum mechanical (black) and classical (magenta) results are presented for comparison.

FIG. 6.

The average position as a function of time for the y-mode of the 2D system in Eq. (29). Shown are the MQC-IVR results obtained with cy = 0.05, cx = 108 (red); cy = 108, cx = 1 (green); and cy = cx = 108 (blue). Quantum mechanical (black) and classical (magenta) results are presented for comparison.

Close modal
TABLE III.

Number of semiclassical trajectories Ntraj required to converge the MQC-IVR calculations of 〈y(t)〉 for different values of Filinov parameters. The maximum absolute statistical error ϵmax < 5 × 10−2 a.u.

cy cx Ntraj
0.05  108  4.0 × 104 
108  1.0  8.5 × 104 
108  108  1.5 × 103 
cy cx Ntraj
0.05  108  4.0 × 104 
108  1.0  8.5 × 104 
108  108  1.5 × 103 

Qualitatively, we explain the results shown in Fig. 6 by examining the role of operator B. If B acts only on a subset of modes at time t to change the corresponding coherent state positions and momenta, then, it is reasonable to expect that a good approximation to the semiclassical correlation function will be obtained by describing this same subset of modes in the quantum limit and the rest of the system in the classical limit. In fact, this treatment bears a strong resemblance to the Forward-Backward IVR approximation previously derived as one limit of the GFB-IVR formulation.30 We expect this feature of MQC-IVR to prove particularly useful in high-dimensional systems where the lighter (more quantum) modes are not the modes being observed.

In this paper, we derive a novel semiclassical approach, MQC-IVR, that provides a uniform dynamic framework in which individual modes of a complex system can be quantized to different extents. We numerically demonstrate that MQC-IVR calculations of correlation functions can tune continuously between computationally expensive quantum-limit SC methods like DHK-IVR and inexpensive, classical-limit SC methods like Husimi IVR. Using a simple 2D model system, we further show that this approach allows for the use of different levels of semiclassical theory for different DOFs, leading to a considerable reduction in computational effort. These features make MQC-IVR a promising approach to incorporate subtle quantum effects due to a few DOFs in large-scale simulations of complex biological and chemical systems.

The authors sincerely thank Bill Miller for several valuable discussions and comments. This work was partially supported by a Young Investigator grant from the Army Research Office (Grant No. W911NFD-13-1-0102). N.A. additionally acknowledges support from a NSF CU-ADVANCE grant and Cornell start-up funds.

Within the MQC-IVR formalism, the general expression for the pre-factor is given by

D = C t ( p 0 , q 0 ) C t ( p t , q t ) det ( K T + i c T J ) 1 2 det ( K T ) 1 2 ,
(A1)

where matrices K and J are defined in Eqs. (12) and (13), respectively. Here, we briefly show how Eq. (A1) can be simplified. The underlying algebra is based on the general expression for the determinant of a block matrix: if A, B, C, and D are components of a block matrix and either A or D is invertible, then

det A B C D = det ( A ) det ( D CA 1 B ) = det ( D ) det ( A BD 1 C ) .
(A2)

We start by considering the denominator in Eq. (A1). We can write

det ( K T ) = det M p q b + i γ 0 M q q b M p p b + i γ 0 M q p b 0 0 0 0 i γ 0 I i 2 γ t 1 2 I 1 2 M p q f i 2 γ t M q q f 1 2 M p p f i 2 γ t M q p f 1 2 I i 2 γ t 1 1 2 M q q f i 2 γ t 1 M p q f 1 2 M q p f i 2 γ t 1 M p p f .
(A3)

Multiplying the last row in Eq. (A3) by t from the left and adding to the third row, followed by moving the last row of the resulting determinant to the second row yields

det ( K T ) = det M p q b + i γ 0 M q q b M p p b + i γ 0 M q p b 0 0 1 2 I i 2 γ t 1 1 2 M q q f i 2 γ t 1 M p q f 1 2 M q p f i 2 γ t 1 M p p f 0 0 i γ 0 I 0 0 M p q f i γ t M q q f M p p f i γ t M q p f = det M p q b + i γ 0 M q q b M p p b + i γ 0 M q p b 1 2 I i 2 γ t 1 det i γ 0 I M p q f i γ t M q q f M p p f i γ t M q p f .
(A4)

The last line of Eq. (A4) is obtained using Eq. (A2), and the product of determinants can be written as

det ( K T ) 1 2 = det ( 2 i γ 0 ) 1 2 C t ( p 0 , q 0 ) C t ( p t , q t ) ,
(A5)

where Ct(p0, q0) and C t ( p t , q t ) are the Herman–Kluk prefactors (Eq. (8)) for the forward and backward trajectories, respectively.

In order to simplify the determinant in the numerator of Eq. (A1), we write the matrix KT + icTJ in terms of 2N × 2N blocks,

K T + i c T J = XM b Y Z Z T M f ,
(A6)

where X = i γ 0 I 0 0 , Y = 0 0 i γ 0 I , Z = i 2 γ t + i c q 1 2 I 1 2 I i 2 γ t 1 + i c p , and Mf/b are the full forward/backward monodromy matrices. Using Eq. (A2) for the determinant yields

det ( K T + i c T J ) = det ( Z ) × det ( Y + XM b Z 1 Z T M f ) .
(A7)

Since components of Z are diagonal matrices, its determinant and inverse can be easily evaluated,

det ( Z ) = det 1 2 ( c q + γ t ) c p + c q ( γ t 1 + c p ) = det 1 2 G ,
(A8)
Z 1 = 2 i 2 γ t 1 + i c p G 1 G 1 G 1 2 i 2 γ t + i c q G 1 .
(A9)

Now, let us consider the second multiplier on the right hand side of Eq. (A7). Due to the structure of matrices X and Y, the resulting determinant can be written as

det ( Y + XM b Z 1 Z T M f ) = det T 11 T 12 i γ 0 I = det ( T 11 i T 12 γ 0 ) ,
(A10)

where T = XMbZ−1ZTMf and subscripts refer to corresponding matrix elements. Substituting Eqs. (A5)–(A10) into Eq. (A1) gives

D = det ( 4 i γ 0 ) 1 2 × det ( G ( T 11 i T 12 γ 0 ) ) 1 2 .
(A11)

Evaluating matrix elements T11 and T12 gives the final expression for the MQC-IVR prefactor,

D = 2 N 2 det ( γ 0 1 G ) 1 2 det 1 2 M p p b i γ 0 M q p b G 1 + I M p p f γ 0 + i M p q f + γ 0 M q q b + i M p q b 1 2 γ t 1 + c p G 1 M p p f γ 0 + i M p q f + 1 2 γ 0 M q q b + i M p q b G 1 + I M q q f i M q p f γ 0 + M p p b i γ 0 M q p b 1 2 γ t + c q G 1 M q q f i M q p f γ 0 1 2 .

The MQC-IVR approximation for the time-dependent correlation function for a momentum-space operator B ˆ B ( p ˆ ) is obtained by evaluating the γt → 0 limit in Eq. (18) to reduce the coherent states at time t to momentum eigenstates. The resulting expression is

C AB ( t ) = 2 2 N π 3 N 2 d q 0 d p 0 d Δ q p 0 q 0 | A ˆ | p 0 q 0 B ( p t ) e i [ S t ( p 0 , q 0 ) + S t ( p t , q t ) ] × e i p t Δ q D p ( p 0 , q 0 , Δ q , c q ) e 1 2 Δ q c q Δ q ,
(B1)

and the corresponding prefactor is given by

D p = det | γ 0 1 c q [ M p p b i γ 0 M q p b M p p f γ 0 + i M p q f + γ 0 M q q b + i M p q b M q q f i M q p f γ 0 + γ 0 M q q b + i M p q b c q 1 M p p f γ 0 + i M p q f ] | 1 2 .
(B2)
1.
J. C.
Tully
,
J. Chem. Phys.
93
,
1061
(
1990
).
2.
A. W.
Jasper
,
S. N.
Stechmann
, and
D. G.
Truhlar
,
J. Chem. Phys.
116
,
5424
(
2002
).
3.
A. V.
Akimov
and
O. V.
Prezhdo
,
Phys. Rev. Lett.
113
,
153003
(
2014
).
4.
P.
Ehrenfest
,
Z. Phys.
45
,
455
(
1927
).
5.
C.
Zhu
,
A. W.
Jasper
, and
D. G.
Truhlar
,
J. Chem. Phys.
120
,
5543
(
2003
).
6.
X.
Li
,
J. C.
Tully
,
H. B.
Schlegel
, and
M. J.
Frisch
,
J. Chem. Phys.
123
,
084106
(
2005
).
7.
J.
Caro
and
L. L.
Salcedo
,
Phys. Rev. A
60
,
842
(
1999
).
8.
C. C.
Martens
and
J.-Y.
Fang
,
J. Chem. Phys.
106
,
4918
(
1997
).
9.
A.
Donoso
and
C. C.
Martens
,
J. Phys. Chem. A
102
,
4291
(
1998
).
10.
R.
Kapral
and
G.
Ciccotti
,
J. Chem. Phys.
110
,
8919
(
1999
).
11.
Q.
Shi
and
E.
Geva
,
J. Chem. Phys.
121
,
3393
(
2004
).
13.
A.
Kelly
,
R.
van Zon
,
J.
Schofield
, and
R.
Kapral
,
J. Chem. Phys.
136
,
084101
(
2012
).
14.
E.
Gindensperger
,
C.
Meier
, and
J. A.
Beswick
,
J. Chem. Phys.
113
,
9369
(
2000
).
15.
E.
Gindensperger
,
C.
Meier
, and
J. A.
Beswick
,
Adv. Quantum Chem.
47
,
331
(
2004
).
16.
S.
Garaschuk
and
M. V.
Volkov
,
J. Chem. Phys.
137
,
074115
(
2012
).
17.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
121
,
3368
(
2004
).
18.
A. R.
Menzeleev
,
N.
Ananth
, and
T. F.
Miller
III
,
J. Chem. Phys.
135
,
074106
(
2011
).
19.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
III
,
Annu. Rev. Phys. Chem.
64
,
387
(
2013
).
20.
A. R.
Menzeleev
,
F.
Bell
, and
T. F.
Miller
III
,
J. Chem. Phys.
140
,
064103
(
2013
).
21.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
100
,
5106
(
1994
).
22.
S.
Jang
and
G. A.
Voth
,
J. Chem. Phys.
111
,
2371
(
1999
).
23.
N.
Ananth
and
T. F.
Miller
III
,
J. Chem. Phys.
133
,
234103
(
2010
).
24.
N.
Ananth
,
J. Chem. Phys.
139
,
124102
(
2013
).
25.
Q.
Shi
and
E.
Geva
,
J. Phys. Chem. A
108
,
6109
(
2004
).
26.
P.
Huo
and
D. F.
Coker
,
Mol. Phys.
110
,
1035
(
2012
).
27.
P.
Huo
,
T. F.
Miller
III
, and
D. F.
Coker
,
J. Chem. Phys.
139
,
151103
(
2013
).
28.
H.
Wang
,
X.
Sun
, and
W. H.
Miller
,
J. Chem. Phys.
108
,
9726
(
1998
).
29.
X.
Sun
,
H.
Wang
, and
W. H.
Miller
,
J. Chem. Phys.
109
,
4190
(
1998
).
30.
M.
Thoss
,
H.
Wang
, and
W. H.
Miller
,
J. Chem. Phys.
114
,
9220
(
2001
).
31.
K. G.
Kay
,
J. Chem. Phys.
100
,
4377
(
1994
).
33.
V. S.
Filinov
,
Nucl. Phys. B
271
,
717
(
1986
).
34.
N.
Makri
and
W. H.
Miller
,
Chem. Phys. Lett.
139
,
10
(
1987
).
35.
N.
Makri
and
W. H.
Miller
,
J. Chem. Phys.
89
,
2170
(
1988
).
36.
B. W.
Spath
and
W. H.
Miller
,
J. Chem. Phys.
104
,
95
(
2004
).
37.
H.
Wang
,
D. E.
Manolopoulos
, and
W. H.
Miller
,
J. Chem. Phys.
115
,
6317
(
2001
).
38.
W. H.
Miller
,
J. Chem. Phys.
125
,
132305
(
2006
).
39.
N.
Ananth
,
C.
Venkataraman
, and
W. H.
Miller
,
J. Chem. Phys.
127
,
084114
(
2007
).
40.
Y.
Zhao
and
N.
Makri
,
Chem. Phys.
280
,
135
(
2002
).
41.
N. J.
Wright
and
N.
Makri
,
J. Phys. Chem. B
108
,
6816
(
2004
).
42.
M. L.
Brewer
,
J. S.
Hulme
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
106
,
4832
(
1997
).
43.
D. T.
Colbert
and
W. H.
Miller
,
J. Chem. Phys.
86
,
1982
(
1992
).