Quantum-Classical Path Integral with Self ... - ACS Publications

Subscriber access provided by Syracuse University Library

Article

Quantum-Classical Path Integral with SelfConsistent Solvent-Driven Reference Propagators Tuseeta Banerjee, and Nancy Makri J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/jp4043123 • Publication Date (Web): 26 Jul 2013 Downloaded from http://pubs.acs.org on July 29, 2013

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Quantum-Classical Path Integral with Self-Consistent Solvent-Driven Reference Propagators

Tuseeta Banerjee1 and Nancy Makri1,2,*

1

Department of Chemistry, University of Illinois,

600 S. Goodwin Avenue, Urbana, Illinois 61801

2

Department of Physics, University of Illinois,

1110 W. Green Street, Urbana, Illinois 61801

* Corresponding author : [email protected], (217) 333-6589

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 38

Abstract Efficient procedures for evaluating the quantum-classical path integral (QCPI) [J. Chem. Phys. 2013, 137, 22A552] are described. The main idea is to identify a trajectory-specific reference Hamiltonian that captures the dominant effects of the classical “solvent” degrees of freedom on the dynamics of the quantum “system”. This time-dependent reference is used to construct a system propagator that is valid for large time increments. Residual “quantum memory” interactions are included via the path integral representation of the density matrix, which converges with large time steps. Two physically motivated reference schemes are considered. The first involves the dynamics of the solvent unperturbed by the system, which forms the basis for the “classical path” approximation. The second is based on solvent trajectories determined self-consistently with the evolution of the system, according to the time-dependent self-consistent field or Ehrenfest model. Application to dissipative two-level systems indicates that both reference schemes allow a substantial increase of the path integral time step, leading to rapid convergence of the path sum. In addition, the time-dependent reference propagators automatically weigh state-to-state coupling against solvent reorganization in the determination of transition probabilities, further enhancing the convergence of the path integral.

Keywords: quantum dynamics, density matrix, interference, decoherence, transitions, nonadiabatic


2

Page 3 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


I. Introduction As is well known, the difficulty of solving the full Schrödinger equation increases very rapidly with the number of coupled degrees of freedom. To date, quantum mechanical effects in the dynamics of condensed phase processes can be accounted for accurately only by making simplified assumptions for the system’s environment.

Another very appealing approach is to restrict the

quantum mechanical treatment to just a few degrees of freedom, approximating the dynamics of the environment by classical trajectories.

Unfortunately, the fundamental differences between the

classical mechanical formalism and the Schrödinger equation make this task extremely difficult. A number of quantum-classical approximations are available. The oldest and crudest of these is generally known as the classical path model. This propagates the wavefunction of the quantum mechanical system along the trajectory of the solvent, which is assumed unperturbed by changes in the state of the system. However, the correct description of chemical processes usually requires accounting for the “back reaction”, i.e., the effects of the system on the classical trajectories of the solvent. The simplest way of achieving this is known as the Ehrenfest model.1 In the context of electronic transitions, the Ehrenfest model implies that the trajectory of the classical nuclei evolves under a force specified by the average of the Born-Oppenheimer potentials describing the various electronic states. As is well known, this assumption can lead to unphysical results.2 The Ehrenfest model is intimately connected to the time-dependent self-consistent field approximation.3,4 Pechukas has described a semiclassical approximation in which the classical trajectory is determined along with the evolution of the quantum wavefunction between initial and final states in a self-consistent manner.5

The procedure is rigorous and accurate, but the self-consistent determination of the

trajectory is very demanding, making the method impractical for large-scale simulation. A different class of quantum-classical approximations aims at practical schemes that avoid the Ehrenfest average. Originally proposed by Tully, surface hopping methods2,6 allow trajectory


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 38

jumps among electronic states with probabilities determined by the magnitude of the nonadiabatic coupling strength.

The surface hopping method corrects the most serious shortcomings of the

Ehrenfest model, allowing for population branching,7 and has been applied to proton transfer reactions,8 multiple proton transfer and proton-coupled electron transfer.9 Other surface hopping schemes utilize the Pechukas force or a mean-field force for short time intervals.10-13 Similar in spirit to surface hopping is the multiple spawning method, which calculates quantum mechanical transition amplitudes using local basis sets that surround classical trajectories.14 Bohm’s quantum trajectory formulation has also been exploited in the development of approximate simulation methods for nonadiabatic dynamics.15-17 Yet another treatment involves solving the quantum-classical Liouville equation subject to a low-order expansion in the ratio of masses corresponding to the quantum and classical particles.18,19 A different idea for extending the capabilities of classical trajectory methods beyond a single Born-Oppenheimer state was originally proposed by Meyer and Miller20 and further developed by Stock and Thoss.21 This approach involves mapping each electronic state on a pair of action-angle variables, thus replacing discrete states by continuous degrees of freedom which are directly amenable to standard classical or semiclassical trajectory treatments on a single potential surface, avoiding the need for a mixed quantum-classical treatment.22 The path integral formulation of quantum mechanics23,24 offers a very attractive starting point for the development of quantum-classical propagation methods, because the local, trajectory-like character of the quantum paths circumvents the averaging involved in many of the above approximations.

However, simulation of a dynamical process typically requires hundreds or

thousands of path integral steps.

Since the number of paths for the quantum system grows

exponentially with the number of time steps, evaluation of the full path sum appears prohibitive. Our group has been pursuing rigorous and practical quantum-classical path integral (QCPI) methods


4

Page 5 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


through the development of forward-backward semiclassical influence functional techniques25,26 and algorithms for dealing with the time-nonlocal nature of the quantum dynamics.27-29 In particular, we have shown30 that the quantum mechanical contribution to nonlocality is associated exclusively with the “back reaction”, i.e., the effects of the quantum system on the environment. For a near-classical solvent, the classical dynamics of the latter provide a reasonable starting point; the interaction with the quantum particle of interest, which perturbs the solvent dynamics, is captured through the “quantum influence function”, which is amenable to more efficient path integral treatment in this framework.31,32 In the present paper we describe procedures for further accelerating QCPI calculations. The basic idea is that simple approximations to the dynamics, such as the classical path approximation (properly averaged with respect to the density matrix of the solvent) or the Ehrenfest model, are excellent on intermediate timescales. We thus use such approximations within each time step of the path integral. By using the solvent trajectory specified by these models as a time-dependent reference, we obtain propagators that capture important dynamics (including the classical decoherence induced by the unperturbed solvent) in a single time step which can exceed significantly the time step allowed by the propagator of the bare system. Quantum-memory corrections arising from the “back reaction” of the system on the bath are captured through a full summation over system paths. Calculations on models show that using these improved reference propagators it is often possible to capture the dynamics of interest within a few path integral time steps, obviating iterative decompositions of the path sum. More generally, we introduce system-independent reference trajectories that take into account the system-solvent interaction according to a chosen ansatz, and use these trajectories to construct improved propagators that allow large time steps. We argue that time-dependent self-consistent field (or Ehrenfest) approximations are often ideally suited to this task, allowing time steps comparable to


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 38

the solvent-induced decoherence time. For our model systems, use of these propagators leads to convergence of the full path integral with relatively minor effort. Thus, we expect this approach to lead to substantial savings when applied to QCPI simulations of condensed phase or biological processes. Section II reviews the QCPI formulation and describes the construction of propagators from a particular ansatz for the solvent dynamics. In section III we discuss two excellent candidates for such an ansatz, the unperturbed solvent and the TDSCF or Ehrenfest models. In section IV, these ideas are illustrated with applications on a two-level system coupled to a harmonic bath, for which accurate results for comparison are available.33,34 Section V presents some concluding remarks.

II. Quantum-Classical Path Integral with Time-Dependent Reference Throughout this paper the quantum mechanical system is described by the coordinate s and momentum ps , while x, p describe collectively the coordinates of the bath or solvent particle(s) with mass mb , which are to be treated via classical trajectories. The total Hamiltonian has the form

Hˆ = H 0 ( sˆ, pˆ s ) + Tb ( pˆ ) + Vb ( xˆ ) + Vint ( sˆ, xˆ )

(2.1)

where H 0 is the Hamiltonian of the quantum system, Tb is the kinetic energy of the bath or solvent, Vb is the potential function for the bare solvent, and Vint is the potential interaction between the solvent and the quantum system. We are interested in the time evolution of the system’s reduced density matrix,


6

Page 7 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ˆ

ˆ

ρ red ( sn± ; nδ t ) = Trb sn+ e− iHt / h ρˆ (0)eiHt / h sn− ,

(2.2)

where δ t is an elementary time step to be specified later. Adopting a classical or semiclassical description of the solvent dynamics, the reduced density matrix can be expressed in the form25,26,31,35

ρ red ( sn± ; nδ t ) = ∫ dx0 ∫ dp0 P ( x0 , p0 )Q( x0 , p0 ; sn± )

(

where P ( x0 , p0 ) is the solvent phase space density and Q x0 , p0 ; sn±

(2.3)

) is the “quantum influence

function”,31,32 which contains all dynamical effects on the solvent due to its interaction with the quantum system. The quantum influence function is given in terms of a sum with respect to all (forward and backward) system paths, according to the expression

(

) ∫

∫

Q x0 , p0 ; sn± = ds0± L dsn±−1 sn+ e −iH 0 δ t / h sn+−1 L s1+ e −iH 0 δ t / h s0+ ρ red ( s0± ;0) ˆ

ˆ

±

±

±

× s0− eiH 0 δ t / h s1− L sn−−1 eiH 0 δ t / h sn− eiΦ ( x0 , p0 ; s0 , s1 ,K, sn )/ h ˆ

ˆ

(2.4)

Here Φ is the trapezoid-discretized action integral along a solvent classical trajectory launched from x0 , p0 and integrated subject to a force dictated by the instantaneous coordinate of the system along the particular path, and δ t is the path integral time step. The precise form of the solvent trajectory and action integral depends on the particular treatment of the solvent dynamics. For example, a classical or linearized path integral (LPI) treatment35,36 leads to trajectories that experience the average of the forces exerted by the system at its instantaneous forward and backward configurations,


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

f ( x, kδ t ) = −

Page 8 of 38

∂  1  Vb ( x) + Vint ( sk+ , x) + Vint ( sk− , x)  ,  ∂x  2 

(

)

(2.5)

while a forward-backward semiclassical treatment of the solvent dynamics (FBSD)25,26 leads to separate forward and backward trajectories. For concreteness we focus on the purely classical treatment of the bath, for which the phase space density is given by the Wigner transform37 of the initial density operator,

P ( x0 , p0 ) = (2π h)−1 ∫ d ∆x0 x0 + 12 ∆x0 ρˆ (0) x0 − 12 ∆x0 e− ip0 ∆x0 / h ,

(2.6)

while the action is given by the expression

 Φ ( x0 , p0 ; s0± , s1± ,K, sn± ) =  12 Vint ( s0+ , x0 ) + 

n −1

∑V

+ + 1 int ( sk , xk ) + 2 Vint ( sn , xn )

k =1

n −1

 − 12 Vint ( s0− , x0 ) − Vint ( sk− , xk ) − 12 Vint ( sn− , xn )  δ t k =1 

(2.7)

∑

where xk = x(kδ t ) . Even though the action integral involves the full potential, only the system-bath interaction terms appear in Eq. (2.7), as the solvent-solvent potential terms enter with opposite signs and thus cancel. It is easily seen from Eq. (2.7) that the action is small in the limit of short time ( n δ t → 0 ) or weak system-bath coupling. Our strategy is to evaluate the system propagators by full quantum mechanics. Note that it is possible to also treat the dynamics of the quantum system at the semiclassical level, arriving at expressions that involve classical trajectories for the coordinates of all (system and solvent) degrees


8

Page 9 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of freedom.38,39 While similar in spirit to QCPI, the main challenges of semiclassical-classical treatments are associated with convergence of the oscillatory phase space integrals and are different from those encountered in our QCPI treatment, as discussed below. The main difficulty in the evaluation of the QCPI expression is the number of terms entering the path sum. This becomes clear by representing the quantum system in terms of M states or sites. The use of discrete variable representations40 of the system coordinate is advantageous in this context, as it allows replacement of each of the 2n integrals in Eq. (2.4) by a discrete sum with respect to a small number of states that have a local, position-like character, leading to a simple representation of the system-bath interaction potential.41 Because each forward-backward system path combination generally leads to a distinct trajectory, Eq. (2.4) requires the evaluation of M 2n classical trajectories. The dependence of the classical trajectories on system paths is the QCPI manifestation of “memory”,31 familiar from influence functional treatments. This memory prevents decomposition of the path sum into a series of single-step operations. On the other hand, direct evaluation of the full path sum is practical only if the number of path integral time steps is small. Consider first the ensemble-averaged classical path (EACP) approximation, in which the classical trajectory from each initial condition is replaced by that of the pure solvent, neglecting effects induced by system-bath interactions. In this approximation all trajectories emanating from the same phase space point are identical, amounting to an external time-dependent potential that drives the system dynamics. Thus, the path sum can be replaced by an iterative procedure that involves n sequential matrix-vector multiplications. Recent work in our group has shown31,32 that the classical path approximation does an excellent job of capturing the decoherence of the quantum system induced by the bath in the limit of high temperature and/or weak system-bath coupling. In fact, the EACP approximation reproduces precisely the real part of the influence functional in the case of a harmonic bath,31 which is dominant in high-temperature/weak-coupling regimes, as the influence of


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 38

the quantum system on the majority of the bath trajectories is very minor in these limits. We also note that the effects of the system-bath interaction accumulate over time, so the EACP approximation can be very good for sufficiently short times even under unfavorable low-temperature or strong coupling conditions.31 This observation can be exploited to increase the path integral time step. To this end, we note that the action Φ depends on the coordinates of the system paths in two ways: (i) through the explicit dependence of the “local” Lagrangian and (ii) through the implicit dependence of the solvent trajectory on the instantaneous force exerted by the system (which makes the action nonlocal in time). We thus partition the action integral into two parts,

(

)

(

)

(

Φ x0 , p0 ; s0± , s1± ,K, sn± = Φ ref s0± , s1± ,K, sn± ; x ref ( x0 , p0 ) + ∆Φ x0 , p0 ; s0± , s1± ,K , sn±

)

(2.8)

of which the first (the “reference” part) depends on the system paths only through the explicit terms in the Lagrangian and the “reference trajectory” x ref (t ′) is the same for all system paths. Using the discrete time notation xkref ≡ x ref (k δ t ), k = 0,K , n , the reference action is

 Φ ref s0± , s1± ,K, sn± ; x ref ( x0 , p0 ) =  12 Vint s0+ , x0ref + 

(

)

(

−

1 V 2 int

(

n −1

) ∑V ( s

s0− , x0ref

int

+ ref k , xk

k =1

n −1

) + ∑V ( int

)+

sk− , xkref (k

k =1

1 V 2 int

δ t)) +

(s

+ ref n , xn

1 V 2 int

(

)

sn− , xnref

)

  

(2.9)

An obvious choice for this reference part is the action along the unperturbed solvent trajectory, but as we argue later, other choices may prove even better. Suppose that interaction with the system


10

Page 11 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


perturbs the solvent trajectory only slightly for times up to ∆t = nref δ t , such that ∆Φ along a path { s0± , s1± ,K , sn±ref } is small and thus be replaced by its trapezoid rule approximation,

∆Φ ( x0 , p0 ; s0± , s1± ,K , sn±ref ) (2.10)

1  ∆Vint ( s0+ , x0 ) + ∆Vint ( sn+ , xn ) − ∆Vint ( s0− , x0 ) − 12 ∆Vint ( sn− , xn )  ∆t ref ref  2

where ∆Vint ( sk± , xk ) = Vint ( sk± , xk ) − Vint ( sk± , xkref ) . With this assumption, Eq. (2.4) becomes

(

) ∫

∫

Q x0 , p0 ; sn±ref = ds0± L dsn±ref −1 e +

ref

× e −iVint ( s1 , x1 −

ref

× eiVint ( s0 , x0 ×e

)δ t / h

)δ t /2 h

− iVint ( sn+ref , xnrefref )δ t /2 h

sn+ref e −iH 0δ t / h sn+ref −1 e ˆ

+

s1+ e− iH 0δ t / h s0+ e −iVint ( s0 , x0 ˆ

−

s0− eiH 0δ t / h s1− eiVint ( s1 , x1 ˆ

iVint ( sn−ref −1 , xnrefref −1 )δ t / h

ref

sn−ref −1 eiH 0δ t / h sn−ref e ˆ

ref

)δ t /2 h

)δ t / h

− iVint ( sn+ref −1 , xnrefref −1 )δ t / h

L

ρred ( s0± ;0) (2.11)

L

iVint ( sn−ref , xnrefref )δ t /2 h

i  ∆Vint ( s0+ , x0 ) +∆Vint ( sn+ref , xn ref ) −∆Vint ( s0− , x0 ) − 12 ∆Vint ( sn−ref , xn ref )  ∆t /2 h

×e 

Inserting the reference trajectory phase factors in the short time propagators and noting that the reference trajectory values xkref are independent of the path integral variables sk± , one sees that Eq. (2.11) is the system density matrix propagated by the time-dependent Hamiltonian

(

Hˆ ref (t ) = Hˆ 0 + Vˆint x ref (t )

)

(2.12)

and the quantum influence function (thus the desired reduced density matrix) is given by the expression


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(

) ∫

† Q x0 , p0 ; sn±ref = ds0± sn+ref U ref (∆t ,0; x0 , p0 ) s0+ ρ red ( s0± ;0) s0− U ref (0, ∆t; x0 , p0 ) sn−ref

×e

−

i ∆t  ∆Vint ( s0+ , x0 ) +∆Vint ( sn+ ref , xn ref ) −∆Vint ( s0− , x0 ) − 12 ∆Vint ( sn− ref , xn ref )  2h 

Page 12 of 38

(2.13)

where U ref is the time evolution operator corresponding to the time-dependent Hamiltonian, Eq. (2.12). The latter, and its matrix elements appearing in Eq. (2.13), can be obtained by solving the time-dependent Schrödinger equation in the basis of system states,

ih

∂ s′ Uˆ ref ( t ; t0 ; x0 , p0 ) s′′ = s′ Hˆ ref (t )Uˆ ref ( t; t0 ; x0 , p0 ) s′′ , ∂t

(2.14)

or, equivalently, by decomposing Eq. (2.11) into nref sequential multiplications that involve the M 2 × M 2 short-time propagator matrix.

We emphasize that the time-dependent reference

Hamiltonian depends on the trajectory initial condition, thus the reference propagators need to be reevaluated for each ( x0 , p0 ) pair. For small ∆t (i.e., nref ~ 1 ), Eq. (2.13) is equivalent to the original expression. As ∆t increases, the effects of the system-bath interaction become more pronounced and may eventually cause the approximation (2.10) to break down. In spite of its simplicity, the procedure described by Equations (2.12)-(2.14) can often capture rather faithfully important features of the dynamics even if ∆t is comparable to the entire time length of interest. For example, if the reference action is

evaluated from trajectories of the pure solvent, as in the EACP approximation, recent work on spinboson models has shown that the procedure captures all bath-induced decoherence of a classical origin.31 The EACP approximation fails to account for “back reaction” effects, which lead to a


12

Page 13 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


distinct bath trajectory for each different forward-backward system path and thus are associated with nonlocality of a quantum mechanical nature. Such effects become increasingly important as the temperature is lowered and/or the system-bath coupling is increased, restricting the accuracy of the classical approximation to short times. If a better reference action could be identified (based on trajectories independent of system paths), the propagation scheme described by Equations (2.13) and (2.14) would be valid over a longer time period. The variables s1± ,K, sn±ref −1 are now redundant, so we relabel the system variable that corresponds to the end of the ∆t interval as s1± . Thus Eq. (2.13) is written as

(

) ∫

† Q x0 , p0 ; s1± = ds0± s1+ U ref (∆t ,0; x0 , p0 ) s0+ ρ red ( s0± ;0) s0− U ref (0, ∆t; x0 , p0 ) s1−

×e

−

i ∆t  ∆Vint ( s0+ , x0 ) +∆Vint ( s1+ , x1 ) −∆Vint ( s0− , x0 ) − 12 ∆Vint ( s1− , x1 )  2h 

(2.15)

Our strategy is to use a physically motivated approximation to obtain solvent trajectories which, while being independent of system paths, capture the dynamics of the composite system-bath interaction as well as possible within a time interval ∆t , and to use these trajectories to construct a time-dependent reference Hamiltonian that will be used to construct a system propagator for this time interval. This propagator is to be used in the full path integral representation of the reduced density matrix to generate accurate results over longer times. As argued above, a reference propagator obtained from a physically motivated approximation should allow much larger time steps compared to the bare system propagator, requiring fewer time slices in the path integral. Furthermore, if the allowed path integral time step is comparable to the span of the quantum nonlocality, evaluation of the path integral via an iterative decomposition32 should be relatively inexpensive.


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 38

Propagators for subsequent time steps are obtained with the same procedure, using the reference trajectory segments in each time interval to construct the driven system Hamiltonian. Using the improved reference propagators, the quantum influence function is written in the form

(

) ∫

Q x0 , p0 ; sN± = ds0± L dsN± −1 sN+ Uˆ ref ( N ∆t ,( N − 1)∆t ; x0 , p0 ) s N+ −1 L

∫

×

s1+

† Uˆ ref ( ∆t ,0; x0 , p0 ) s0+ ρ red ( s0± ;0) s0− Uˆ ref ( 0, ∆t; x0 , p0 ) s1− L i

† × sN− −1 Uˆ ref ( ( N − 1)∆t , N ∆t; x0 , p0 ) sN− e h

(2.16)

∆Φ ( x0 , p0 ; s0± ,K, s N± )

where ∆Φ is given by Eq. (2.8) but is discretized with the new time step ∆t . Besides allowing large time steps, an advantage of Eq. (2.16) compared to Eq. (2.4) is that the magnitude of the reference propagators is not independent of the solvent dynamics. This implies that a large number (often the majority) of system paths are statistically insignificant and do not need to be included. Such selection in path integral calculations can reduce the required number of terms by many orders of magnitude.42-44 In the language of nonadiabatic dynamics, state-to-state transitions are controlled by the magnitude of solvent reorganization. This is crucial for efficiency, as the interplay between nonadiabatic transitions and dissipative effects is known to play an important role in determining reaction rates.45,46 Since the system propagators are to be evaluated from truncated basis expansions,41,47,48 the use of improved time-dependent reference schemes reduces the overall phase outside of the propagators and thus leads to substantial smoothing of the QCPI integrand, allowing faster convergence of the integral with respect to solvent initial conditions. Eq. (2.16) is in these aspects similar to the expression obtained by Lambert and Makri32, but is less demanding and easier to apply.


14

Page 15 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


To implement Eq. (2.16) in its full path sum form, one propagates a solvent trajectory with the initial condition x0 , p0 following the chosen ansatz. Using the interaction potential along this trajectory, one solves the time-dependent Schrödinger equation to obtain the system propagators. Next, the full solvent-interaction trajectory is obtained for each forward-backward system path combination, along with its action integral. The difference of actions along the reference and systemdriven trajectory gives the net action ∆Φ . The phase obtained from this action difference captures the corrections to the reference approximation, Eq. (2.13).

III. Unperturbed Solvent and Self-Consistent Reference The simplest ansatz for the reference trajectories is obtained by integrating the equations of motion for the pure solvent. This can be done by turning off the interaction potential or by fixing the coordinates of the quantum system. (For example, one could place the quantum mechanical atom at the potential minimum of the reactant potential well or at the transition state.) Recent work in our group31,32 examined the accuracy of the unperturbed solvent ansatz, Eq. (2.13) (which is equivalent to the EACP approximation), as a function of propagation time on harmonic bath models. These calculations concluded that this simple treatment can lead to nearly quantitative results at high temperatures or when the system-bath interaction is weak. Under conditions that deviate significantly from this “classical” regime, the error of the classical path approximation grows over time, eventually causing the EACP approximation to yield incorrect populations.

However, even under such

unfavorable conditions, the unperturbed solvent ansatz amounts to a substantial improvement compared to results obtained using the propagator of the bare quantum system. As a result, it can be used to construct a propagator that is valid over much larger time steps than those allowed with a propagator obtained from the conventional splitting of the Hamiltonian.


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 38

To go beyond the unperturbed ansatz, one needs to obtain trajectories that take into account the interaction between system and solvent, but which do not depend on each system path. Perhaps the best way of accomplishing this task is to have the solvent follow the average system path. This prescription is known as the Ehrenfest model and can also be obtained from the time-dependent selfconsistent field (TDSCF) approximation3,4 by imposing a classical trajectory description of the solvent. The TDSCF approximation is obtained from the time-dependent variational principle49 with a product wavefunction ansatz and thus leads to optimized single-particle wavefunctions. It allows semiquantitative descriptions of energy transfer,3,4 but is unable to properly account for wavepacket splitting.50 The quantum-classical Ehrenfest version,51,52 where the solvent dynamics are obtained from an ensemble of classical trajectories that are local in space, can correct some of the main flaws of the fully quantum TDSCF treatment by avoiding the averaging with respect to the solvent in the determination of the force that drives the quantum particle. Still, the Ehrenfest model is generally inadequate for describing proton transfer events or electronic transitions, often leading to incorrect branching ratios and equilibrium populations. However, it seems an excellent candidate for obtaining a reference propagator for use in the QCPI framework, as it is likely to allow larger time steps than possible by bare system or classical path treatments. For each initial condition sampled from the appropriate phase space distribution, the Ehrenfest model integrates a classical trajectory subject to the force

−

∂ (Vb ( xt ) + Vint (sˆ, xt ) ∂xt

)

(2.17)

where Vint ( sˆ, xt ) is the average with respect to the quantum system of the interaction potential at the instantaneous position of the classical trajectory, i.e.,


16

Page 17 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(

Vint ( sˆ, xt ) = Trsys ρˆ sys (t ) Vint ( sˆ, xt )

)

(2.18)

where ρˆsys (t ) is the density matrix of the system in the Ehrenfest approximation. The latter is propagated according to the von Neumann equation

ih

d ρˆsys (t ) =  Hˆ 0 + Vint ( sˆ, xt ), ρˆsys (t )  . dt

(2.19)

IV. Application to Model Systems We illustrate the ideas presented in the last section with calculations on the frequently revisited “spin-boson” model, which consists of a two-level system (TLS) coupled to a harmonic bath.53 The quantum system, described by the Hamiltonian

H 0 = − hΩ σ x ,

(3.1)

is characterized by the tunneling splitting 2hΩ and is coupled to a dissipative bath described by the terms

 pˆ 2j 1  Hb = ∑ + mω 2j xˆ 2j , Vint = −σ z ∑ c j xˆ j .   2 j  2m j 


(3.2)

17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 38

where σ x and σ z are the standard 2 × 2 Pauli spin matrices. With appropriate initial conditions, the Hamiltonian of Equations (3.1)-(3.2) can serve as an excellent model of condensed-phase electron transfer or vibrational relaxation. In the calculations presented below we assume that the system is initially in the right-localized TLS state and that the bath is at thermal equilibrium at the temperature

k BT = β −1 . The bath frequencies and coupling coefficients are described collectively by the spectral density function54

J (ω ) =

c2

π

j δ (ω − ω j ) . ∑ 2 mω j

j

(3.3)

j

In particular, we choose the common Ohmic form,53

1 J (ω ) = π hξω e −ω /ωc , 2

(3.4)

where the dimensionless Kondo parameter ξ provides a measure of the dissipation strength (i.e., the magnitude of the system-bath coupling) and ωc gives the maximum of the spectral density. In order to apply the QCPI methodology, we approximate the bath by a sum of 60 harmonic oscillators whose frequencies and coupling coefficients are evaluated using a logarithmic discretization of the spectral density55 with a maximum frequency given by ωmax = 4ωc . In order to assess the performance of the reference trajectory schemes described in the previous section, we show in Figure 1 the average TLS position σ z (t ) obtained via a single path


18

Page 19 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


integral propagation step. The diagonal form of the initial density matrix of the system implies s0+ = s0− , while the trace operation imposes the condition sN+ = sN− . Since both endpoints of the forward and backward system paths are the same, these paths are identical within the single step treatment and ∆Φ = 0 . In this case the two single-step QCPI schemes discussed in the previous section are equivalent to the corresponding approximations, i.e., the EACP approximation given by Eq. (2.11)-(2.14) and the self-consistent (Ehrenfest) approximation. The calculations presented in Figure 1 compare results obtained with the two (unperturbed bath and Ehrenfest-based) time-dependent reference propagators, as well as the conventional propagator splitting based on the bare system Hamiltonian, for a symmetric TLS. The number of Monte Carlo samples was adjusted to produce statistical error bars comparable to or smaller than the size of the markers. This precision was achieved with 40,000 Monte Carlo points in the case of the bare (time-independent) system propagator, and 10,000-20,000 points when the time dependent reference propagators were used. The single-step QCPI results are compared to numerically exact results obtained with the iterative path integral methodology developed earlier in our group,33,34 which is based on the analytical Feynman-Vernon influence functional.56 In all cases, the results obtained using the single-step system propagator derived from the bare TLS Hamiltonian are accurate only for a small fraction of the TLS period and are qualitatively wrong after t

0.5 Ω−1 or (in the case of Fig.

1a) even earlier. Figure 1a shows the average TLS position in a high-temperature, overdamped regime, characterized by hΩβ = 0.2 , ωc = 2.5Ω and ξ = 1.2 . In this regime the behavior of the bath is almost completely classical, and the Wigner function may be replaced by the Boltzmann density. Because of the large system-bath coupling strength, the bare system propagator fails after about 0.2 Ω −1 .

This is in line with earlier QCPI calculations, where a time step ∆t = 0.125 Ω −1 was


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 38

necessary to converge the path integral.31 By contrast, the EACP approximation yields reasonable results over the entire propagation time. (The small deviations in the EACP approximation can be attributed to high-frequency bath oscillators in the tail of the spectral density, for which the bath trajectories are significantly perturbed by the system.)

With these parameters, the Ehrenfest

approximation yields results indistinguishable from those of the full path integral treatment. In Fig. 1b we show results at a lower temperature, hΩβ = 2.5 , and intermediate coupling strength characterized by ξ = 0.6 .

Even though the bath has the same frequency distribution,

( ωc = 2.5Ω ), it exhibits more pronounced quantum mechanical effects at this lower temperature. The unperturbed solvent approximation leads to accurate results for Ωt 1 , after which it is seen to exaggerate the oscillatory nature of the TLS dynamics. In this case the Ehrenfest model reproduces the underdamped oscillation of the TLS average position nearly quantitatively for all times. A strongly quantum mechanical, high-frequency, low-temperature bath, characterized by

ωc = 5 Ω , hΩβ = 5 and intermediate coupling strength ( ξ = 0.3 ), is employed in the calculations presented in Figure 1c. In this case the oscillators near the peak of the spectral density are practically at zero temperature. In this regime neither the EACP approximation nor the Ehrenfest approximation can reproduce the exact quantum mechanical results in a quantitative fashion. Specifically, the unperturbed solvent approximation exaggerates the oscillatory character of the motion because it neglects decoherence associated with quantum mechanical memory. On the other hand, the Ehrenfest approximation leads to faster damping of the TLS dynamics. Thus full QCPI calculations are necessary in this regime. Because both approximations seem adequate for times up to t

Ω −1 , either

one seems an excellent candidate for constructing a reference propagator that should allow large path integral time steps.


20

Page 21 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


An intuitive explanation of the improvement offered by the mean field (Ehrenfest) approximation emerges by examining the bath trajectories subject to a time-dependent force along particular forward-backward system paths. Figure 2 shows such trajectories in the case of a single bath degree of freedom coupled to the TLS. (For clarity, the system-bath coupling used to obtain the trajectories in Fig. 2 has been increased substantially compared to the values employed in the calculations of Fig. 1.) The momentum of these trajectories exhibits kinks where the force from the quantum system changes discontinuously.

These kinks are more prominent in lower-energy

trajectories, which are perturbed more strongly by the interaction of the bath with the system. This observation is the physical reason for the less significant role of the “back reaction” under high temperature conditions. Also shown in Fig. 2 are the trajectory of the free solvent and the Ehrenfest trajectory subject to the force from the same forward-backward system paths. It is seen that the Ehrenfest trajectory tracks the exact trajectory of the bath rather faithfully for a significant portion of the oscillation period, but eventually fails. As a consequence, system propagators that incorporate the Ehrenfest reference should be accurate for comparable time steps. We emphasize that even though Fig. 2 seems to imply that the unperturbed solvent trajectories should constitute a poor approximation, they nevertheless provide a physically sound time-dependent reference for calculation of improved system propagators that are accurate over time steps much larger than those possible when the propagator is constructed from the time-independent Hamiltonian of the bare system. Further, we note that the Ehrenfest approximation benefits significantly from the early-time localization of the quantum system imposed by the initial condition. This advantage is lost after the first time step in a multi-step path integral calculation, as the system-averaged force can even have the wrong sign for many of the system paths. As a result, not only does the Ehrenfest treatment smooth out the kinky trajectories, it also tends to overcorrect the deficiencies of the unperturbed solvent approximation, leading to deviations in the opposite direction.


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 38

In Figure 3 we present full QCPI calculations using propagators derived with the unperturbed solvent or the self-consistent (Ehrenfest) model as the reference Hamiltonian. The parameters are those employed in the calculations of Fig. 1. As can be seen in Fig. 3, the use of improved timedependent reference propagators yields converged results with large time steps. The most notable gain is again observed in the high-temperature case shown in Fig. 3a, where the Ehrenfest reference produces results that are practically exact with a single time step. Thus, a single trajectory is required for each initial condition in this case. The propagator based on the unperturbed solvent reference leads to converged QCPI results with ∆t = Ω −1 . As seen in the figure, such a time step is by far too large for the standard propagator splitting, which requires a step size ∆t = 0.125 Ω −1 for convergence. Propagation to the same final time with the time-independent system propagator would thus require 40 path integral steps and 280 trajectories for each initial condition. Under the more challenging conditions of Figures 3b and 3c, the QCPI calculations with the two time-dependent reference models converged with ∆t = 0.625 Ω −1 and ∆t = 0.75 Ω −1 respectively. For comparison, the time step required with the standard propagator discretization was ∆t = 0.25 Ω −1 . Perhaps the most significant consequence of a threefold gain in step size is that the important features of the system dynamics could be captured with only 5 path integral time steps (i.e., 210 trajectories from each initial condition), obviating the need for an iterative decomposition of the path sum.


22

Page 23 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


V. Concluding Remarks The path integral formulation of quantum-classical (or quantum-semiclassical) dynamics provides a rigorous approach to the dynamics of quantum mechanical processes in complex polyatomic environments. Recent developments indicate that QCPI calculations on many systems of interest are well within the realm of present computational power. Because the QCPI formulation is free of assumptions (besides the use of classical trajectories for the solvent), its application will produce useful benchmarks and offer valuable insights. In the present paper we presented a simple, yet powerful approach that can accelerate significantly the convergence of QCPI calculations. The idea is to incorporate a physically motivated approximation to the solvent dynamics into the propagator for the quantum system. Two such candidates were considered: the unperturbed solvent approximation and the time-dependent selfconsistent field (Ehrenfest) model. Even though both of these generally fail to capture the correct dynamics at long times, they are excellent for constructing single step QCPI propagators which are accurate for time steps that exceed considerably the step allowed by conventional time-independent propagators.

The interference among the free solvent or Ehrenfest trajectories corrects the

deficiencies of these approximations, leading to correct branching ratios and long-time populations that satisfy the detailed balance property. The use of propagators with larger time steps implies that the desired propagation time (or the decoherence time, in the case of iterative evaluation of the path sum) is reached with fewer time slices, and thus the number of classical trajectories that need to be integrated is decreased dramatically. Further advantages of time-dependent reference propagators are the dependence of quantum transitions on solvent configuration, which invites the use of path


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 38

filtering techniques, and the smoothing of the integrand achieved via the incorporation of a portion of the phase in the system propagators, which implies convergence with fewer Monte Carlo samples. The applications treated in section IV focused on the model of a dissipative two-level system, for which accurate results for comparison could be generated. In this case, the quadratic nature of the bath leads to an analytical expression for the Wigner density. Applications to complex systems will necessitate a numerical evaluation of the phase space function. This is not an easy task, but a number of approximate techniques for evaluating the Wigner transform are available.36,57-61 Also, the classical Boltzmann function should provide a sufficiently accurate alternative in many cases. We point out that it should be straightforward to incorporate the time-dependent reference schemes described in this paper in the semiclassical formulation of QCPI, which employs separate trajectories along the forward and backward branches of the time contour, along with a coherent state phase space density, 62-64

which can be evaluated via an imaginary time path integral calculation

or (approximately) by a

semiclassical treatment.65 A particularly promising implementation of QCPI is in combination with iterative Monte Carlo (IMC) algorithms.66-68 The latter incorporate powerful Metropolis Monte Carlo sampling procedures69 in iterative decompositions of the path integral. The IMC methodology is quite efficient for systems with 10-20 degrees of freedom; thus, its use within the QCPI framework should allow the treatment of polyatomic systems by full quantum mechanics, retaining the classical trajectory description of the condensed phase or biological environment. Applications of QCPI to proton transfer in solution are in progress in our group. Preliminary work suggests that the methodology described in this paper, in particular the construction of propagators based on the unperturbed solvent reference, appears quite efficient as well as practical for such calculations.


24

Page 25 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Acknowledgment

This material is based upon work supported by the National Science Foundation under Award CHE 11-12422.


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 38

References

(1)

Ehrenfest, P. Bemerkung Über Die Angenäherte Gültigkeit Der Klassischen

Mechanik Innerhalb Der Quantenmechanik. Z. Phys. 1927, 45, 455-457. (2)

Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys.

1990, 93, 1061-1071. (3)

Gerber, R. B.; Ratner, M. A. Self-Consistent-Field Methods for Vibrational

Excitations in Polyatomic Systems. Adv. Chem. Phys. 1988, 70, 97-132. (4)

Billing, G. D. Classical Path Method in Inelastic and Reactive Scattering. Int.

Rev. Phys. Chem. 1994, 13, 309. (5)

Pechukas, P. Time-Dependent Semiclassical Scattering Theory. I. Potential

Scattering. Phys. Rev. 1969, 181, 166-185. (6)

Tully, J. C.; Preston, R. K. Trajectory Surface Hopping Approach to

Nonadiabatic Molecular Collisions: The Reaction of H+ with D2. J. Chem. Phys. 1971, 55, 562-572. (7)

Herman, M. F. Dynamics by Semiclassical Methods. Ann. Rev. Phys. Chem.

1994, 45, 83. (8)

Hammes-Schiffer, S.; Tully, J. C. Proton Transfer in Solution: Molecular

Dynamics with Quantum Transitions. J. Chem. Phys. 1994, 101, 4657-4667. (9)

Hammes-Schiffer, S. Theory of Proton-Coupled Electron Transfer in Energy

Conversion Processes. Acc. Chem. Res. 2009, 42, 1881-1889.


26

Page 27 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(10)

Webster, F. J.; Schnitker, J.; Friedrichs, M. S.; Friesner, R. A.; Rossky, P. J.

Solvation Dynamics of the Hydrated Electron: A Nonadiabatic Quantum Simulation. Phys. Rev. Lett. 1991, 66, 3172-3175. (11)

Webster, F.; Rossky, P. J.; Friesner, R. A. Nonadiabatic Processes in

Condensed Matter: Semiclassical Theory and Implementation. Comput. Phys. Commun. 1991, 63, 494-522. (12)

Coker, D. F.; Xiao, L. Methods for Molecular Dynamics with Nonadiabatic

Transitions. J. Chem. Phys. 1995, 102, 496-510. (13)

Prezhdo, O. V.; Rossky, P. J. Mean Field Molecular Dynamics with Surface

Hopping. J. Chem. Phys. 1997, 107, 825-834. (14)

Ben-Nun, M.; Quenneville, J.; Martinez, T. J. Ab Initio Multiple Spawning:

Photochemistry from First Principles Molecular Dynamics. J. Phys. Chem. 2000, 104, 51615175. (15)

Burant, J. C.; Tully, J. C. Nonadiabatic Dynamics Via the Classical Limit

Schrodinger Equation. J. Chem. Phys. 2000, 112, 6097-6103. (16)

Prezhdo, O. V.; Brooksby, C. Quantum Backreaction through the Bohmian

Particle. Phys. Rev. Lett. 2001, 86, 3215-3219. (17)

Wyatt, R. E.; Lopreore, C. L.; Parlant, G. Electronic Transitions with

Quantum Trajectories. J. Chem. Phys. 2001, 114, 5113-5116. (18)

Donoso, A.; Martens, C. C. Simulation of Coherent Nonadiabatic Dynamics

Using Classical Trajectories. J. Phys. Chem. A 1998, 102, 4291.


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(19)

Page 28 of 38

Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Dynamics. J. Chem. Phys.

1999, 110, 8919-8929. (20)

Meyer, H.-D.; Miller, W. H. A Classical Analog for Electronic Degrees of

Freedom in Nonadiabatic Collision Processes. J. Chem. Phys. 1979, 70, 3214-3223. (21)

Stock, G.; Thoss, M. Semiclassical Description of Nonadiabatic Quantum

Dynamics. Phys. Rev. Lett. 1997, 78, 578-581. (22)

Miller, W. H. Electronically Nonadiabatic Dynamics Via Semiclassical Initial

Value Methods. J. Phys. Chem. 2009, 113, 1405-1415. (23)

Feynman, R. P. Space-Time Approach to Non-Relativistic Quantum

Mechanics. Rev. Mod. Phys. 1948, 20, 367-387. (24)

Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals;

McGraw-Hill: New York, 1965. (25)

Makri, N.; Thompson, K. Semiclassical Influence Functionals for Quantum

Systems in Anharmonic Environments. Chem. Phys. Lett. 1998, 291, 101-109. (26)

Thompson, K.; Makri, N. Influence Functionals with Semiclassical

Propagators in Combined Forward-Backward Time. J. Chem. Phys. 1999, 110, 1343-1353. (27)

Makri, N. Quantum Dissipative Systems: A Numerically Exact Methodology.

J. Phys. Chem. 1998, 102, 4414-4427. (28)

Makri, N. Iterative Evaluation of the Path Integral for a System Coupled to an

Anharmonic Bath. J. Chem. Phys. 1999, 111, 6164-6167. (29)

Makri, N. Path Integral Renormalization for Quantum Dissipative Dynamics

with Multiple Timescales. Mol. Phys. 2012, 110, 1001-1007.


28

Page 29 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(30)

Makri, N. Dynamics of Reduced Density Matrices: Classical Memory Vs.

Quantum Nonlocality. J. Chem. Phys. 1998, 109, 2994-2998. (31)

Lambert, R.; Makri, N. Quantum-Classical Path Integral: Classical Memory

and Weak Quantum Nonlocality. J. Chem. Phys. 2012, 137, 22A552. (32)

Lambert, R.; Makri, N. Quantum-Classical Path Integral: Numerical

Formulation. J. Chem. Phys. 2012, 137, 22A553. (33)

Makri, N.; Makarov, D. E. Tensor Multiplication for Iterative Quantum Time

Evolution of Reduced Density Matrices. I. Theory. J. Chem. Phys. 1995, 102, 4600-4610. (34)

Makri, N.; Makarov, D. E. Tensor Multiplication for Iterative Quantum Time

Evolution of Reduced Density Matrices. Ii. Numerical Methodology. J. Chem. Phys. 1995, 102, 4611-4618. (35)

Shi, Q.; Geva, E. A Relationship between Semiclassical and Centroid

Correlation Functions. J. Chem. Phys. 2003, 118, 8173-8183. (36)

Poulsen, J. A.; Nyman, G.; Rossky, P. J. Practical Evaluation of Condensed

Phase Quantum Correlation Functions: A Feynman--Kleinert Variational Linearized Path Integral Method. The Journal of Chemical Physics 2003, 119, 12179-12193. (37)

Wigner, E. J. Calculation of the Rate of Elementary Association Reactions.

Chem. Phys. 1937, 5, 720. (38)

Sun, X.; Miller, W. H. Mixed Semiclassical-Classical Approaches to the

Dynamics of Complex Molecular Systems. J. Chem. Phys. 1997, 106, 916-927. (39)

Thompson, K.; Makri, N. Rigorous Forward-Backward Semiclassical

Formulation of Many-Body Dynamics. Phys. Rev. E 1999, 59, R4729-R4732.


29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(40)

Page 30 of 38

Bacic, Z.; Light, J. C. Theoretical Methods for Rovibrational States of Floppy

Molecules. Annu. Rev. Phys. Chem. 1989, 40, 469-498. (41)

Topaler, M.; Makri, N. System-Specific Discrete Variable Representations for

Path Integral Calculations with Quasi-Adiabatic Propagators. Chem. Phys. Lett. 1993, 210, 448. (42)

Sim, E.; Makri, N. Filtered Propagator Functional for Iterative Dynamics of

Quantum Dissipative Systems. Comp. Phys. Commun. 1997, 99, 335-354. (43)

Sim, E. Quantum Dynamics for a System Coupled to Slow Baths: On-the-Fly

Filtered Propagator Method. J. Chem. Phys. 2001, 115, 4450-4456. (44)

Lambert, R.; Makri, N. Memory Path Propagator Matrix for Long-Time

Dissipative Charge Transport Dynamics. Mol. Phys. 2012, 110, 1967-1975. (45)

Frauenfelder, H.; Wolynes, P. G. Rate Theories and Puzzles of Heme Protein

Kinetics. Science 1985, 228, 337. (46)

Wolynes, P. G. Dissipation, Tunneling and Adiabaticity Criteria for Curve

Crossing Problems in the Condensed Phase. J. Chem. Phys. 1987, 86, 1957-1966. (47)

Makri, N. Improved Feynman Propagators on a Grid and Non-Adiabatic

Corrections within the Path Integral Framework. Chem. Phys. Lett. 1992, 193, 435-444. (48)

Ilk, G.; Makri, N. Real Time Path Integral Methods for a System Coupled to

an Anharmonic Bath. J. Chem. Phys. 1994, 101, 6708. (49)

Dirac, P. A. Note on the Exchange Phenomena in the Thomas Atom. Proc.

Cambridge Philos. Soc. 1930, 26.


30

Page 31 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(50)

Makri, N.; Miller, W. H. Time-Dependent Self-Consistent Field

Approximation for a Reaction Coordinate Coupled to a Harmonic Bath: Single and Multiple Configuration Treatments. J. Chem. Phys. 1987, 87, 5781-5787. (51)

Wahnstrom, G.; Carmeli, B.; Metiu, H. The Calculation of the Thermal Rate

Coefficient by a Method Combining Classical and Quantum Dynamics. J. Chem. Phys. 1988, 88, 2478-2491. (52)

Haug, K.; Metiu, H. A Test of the Possibility of Calculating Absorption

Spectra by Mixed Quantum-Classical Methods. J. Chem. Phys. 1992, 97, 4781-4791. (53)

Leggett, A. J.; Chakravarty, S.; Dorsey, A. T.; Fisher, M. P. A.; Garg, A.;

Zwerger, M. Dynamics of the Dissipative Two-State System. Rev. Mod. Phys. 1987, 59, 1-85. (54)

Caldeira, A. O.; Leggett, A. J. Path Integral Approach to Quantum Brownian

Motion. Physica A 1983, 121, 587-616. (55)

Makri, N. The Linear Response Approximation and Its Lowest Order

Corrections: An Influence Functional Approach. J. Phys. Chem. 1999, 103, 2823-2829. (56)

Feynman, R. P.; Vernon, F. L. The Theory of a General Quantum System

Interacting with a Linear Dissipative System. Annals of Physics 1963, 24, 118. (57)

Wang, H.; Sun, X.; Miller, W. H. Semiclassical Approximations for the

Calculation of Thermal Rate Constants for Chemical Reactions in Complex Molecular Systems. J. Chem. Phys. 1998, 108, 9726-9736. (58)

Shi, Q.; Geva, E. Semiclassical Theory of Vibrational Energy Relaxation in

the Condensed Phase. J. Phys. Chem. A 2003, 107, 9059-9069.


31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(59)

Page 32 of 38

Frantsuzov, P. A.; Neumaier, A.; Mandelshtam, V. A. Gaussian Resolutions

for Equilibrium Density Matrices. Chem. Phys. Lett. 2003, 381, 117. (60)

Liu, J.; Miller, W. H. Using the Thermal Gaussian Approximation for the

Boltzmann Operator in Semiclassical Initial Value Time Correlation Functions. The Journal of Chemical Physics 2006, 125, 224104. (61)

Liu, J.; Miller, W. H. A Simple Model for the Treatment of Imaginary

Frequencies in Chemical Reaction Rates and Molecular Liquids. J. Chem. Phys. 2009, 131, 074113. (62)

Jezek, E.; Makri, N. Finite Temperature Correlation Functions Via Forward-

Backward Semiclassical Dynamics. J. Phys. Chem. 2001, 105, 2851-2857. (63)

Makri, N. Monte Carlo Evaluation of Forward-Backward Semiclassical

Correlation Functions with a Quantized Coherent State Density. J. Phys. Chem. B 2002, 106, 8390-8398. (64)

Nakayama, A.; Makri, N. Simulation of Dynamical Properties in Normal and

Superfluid Helium. Proc. Nat. Acad Sci. USA 2005, 102, 4230-4234. (65)

Makri, N.; Miller, W. H. Coherent State Semiclassical Initial Value

Representation for the Boltzmann Operator in Thermal Correlation Functions. J. Chem. Phys.

2002, 116, 9207-9212. (66)

Jadhao, V.; Makri, N. Iterative Monte Carlo for Quantum Dynamics. J. Chem.

Phys. (Communication) 2008, 129, 161102. (67)

Jadhao, V.; Makri, N. Iterative Monte Carlo with Bead-Adapted Sampling for

Complex Time Correlation Functions. J. Chem. Phys. 2010, 132, 104110.


32

Page 33 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(68)

Jadhao, V.; Makri, N. Iterative Monte Carlo Path Integral with Optimal Grids

from Whole-Necklace Sampling. J. Chem. Phys. 2010, 133, 114105. (69)

Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, H.; Teller, E.

Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 10871092.


33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 38

Figure Captions

Fig. 1.

Average position of a symmetric TLS from single-step calculations. Green hollow triangles: propagator for the bare TLS. Blue solid circles: time-dependent reference propagator obtained from unperturbed solvent trajectories (EACP approximation).

Red solid squares: time-dependent reference propagator

obtained from Ehrenfest trajectories. Black line: exact results obtained via the iterative path integral methodology based on the analytical Feynman-Vernon influence functional. (a) ωc = 2.5Ω , hΩβ = 0.2 ,

ξ = 1.2 (b) ωc = 2.5Ω , hΩβ = 2.5 , ξ = 0.6 . (c) ωc = 5Ω , hΩβ = 5 , ξ = 0.3 .

Fig. 2.

Classical trajectories along arbitrarily chosen forward-backward system paths for a single-oscillator bath with ω = 2.5Ω strongly coupled to the TLS.

Blue line: free bath trajectory.

Red line: trajectory

following the Ehrenfest system force. Black line: full trajectory, following the actual force exerted by the TLS.

Fig. 3.

Average position of a symmetric TLS from QCPI calculations. Green hollow triangles: propagator for the bare TLS. Blue solid circles: time-dependent reference propagator obtained from unperturbed solvent trajectories (EACP approximation).

Red solid squares: time-dependent reference propagator obtained from Ehrenfest

trajectories. Black line: exact results obtained via the iterative path integral methodology based on the analytical Feynman-Vernon influence functional.

(a) ωc = 2.5Ω , hΩβ = 0.2 , ξ = 1.2 , Ω∆t = 1 . (b) ωc = 2.5Ω ,

hΩβ = 2.5 , ξ = 0.6 , Ω∆t = 0.625 . (c) ωc = 5Ω , hΩβ = 5 , ξ = 0.3 , Ω∆t = 0.75 .


34

Page 35 of 38

15

10

5

p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

-5

-10

-15 0

0.5

1

1.5

2

2.5

3

Ωt


35


1

1

(b)

(a) 0.5

0.5

0

0

-0.5

-0.5

-1

-1 0

1

2

3

4

5

0

Ωt

1

2

3

4

5

Ωt

1

(c) 0.5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 38

0

-0.5

-1 0

1

2

3

4

5

Ωt

Fig. 1


36

6

15

4

10

2

5

0

0

p

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


x

Page 37 of 38

-2

-5

-4

-10

-15

-6 0

0.5

1

1.5

2

2.5

3

0

0.5

1

1.5

2

2.5

3

Ωt

Ωt

Fig. 2


37


1

1

(a)

(b) 0.5

0.5

0

-0.5

0

-0.5

-1

-1 0

1

2

3

4

5

Ωt

0

1

2

Ωt

3

4

5

1

(c) 0.5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 38

0

-0.5

-1 0

1

2

Ωt

3

4

5

Fig. 3


38

Quantum-Classical Path Integral with Self ... - ACS Publications

Recommend Documents