THE JOURNAL OF CHEMICAL PHYSICS 138, 144106 (2013)

A mixed quantum-classical Liouville study of the population dynamics in a model photo-induced condensed phase electron transfer reaction Najeh Rekik,1 Chang-Yu Hsieh,2 Holly Freedman,1 and Gabriel Hanna1,a) 1 2

Department of Chemistry, University of Alberta, Edmonton, Alberta T6G 2G2, Canada Department of Chemistry, University of Toronto, Toronto, Ontario M5S 3H6, Canada

(Received 25 October 2012; accepted 19 March 2013; published online 12 April 2013) We apply two approximate solutions of the quantum-classical Liouville equation (QCLE) in the mapping representation to the simulation of the laser-induced response of a quantum subsystem coupled to a classical environment. These solutions, known as the Poisson Bracket Mapping Equation (PBME) and the Forward-Backward (FB) trajectory solutions, involve simple algorithms in which the dynamics of both the quantum and classical degrees of freedom are described in terms of continuous variables, as opposed to standard surface-hopping solutions in which the classical degrees of freedom hop between potential energy surfaces dictated by the discrete adiabatic state of the quantum subsystem. The validity of these QCLE-based solutions is tested on a non-trivial electron transfer model involving more than two quantum states, a time-dependent Hamiltonian, strong subsystembath coupling, and an initial energy shift between the donor and acceptor states that depends on the strength of the subsystem-bath coupling. In particular, we calculate the time-dependent population of the photoexcited donor state in response to an ultrafast, on-resonance pump pulse in a three-state model of an electron transfer complex that is coupled asymmetrically to a bath of harmonic oscillators through the optically dark acceptor state. Within this approach, the three-state electron transfer complex is treated quantum mechanically, while the bath oscillators are treated classically. When compared to the more accurate QCLE-based surface-hopping solution and to the numerically exact quantum results, we find that the PBME solution is not capable of qualitatively capturing the population dynamics, whereas the FB solution is. However, when the subsystem-bath coupling is decreased (which also decreases the initial energy shift between the donor and acceptor states) or the initial shift is removed altogether, both the PBME and FB results agree better with the QCLE-based surface-hopping results. These findings highlight the challenges posed by various conditions such as a time-dependent external field, the strength of the subsystem-bath coupling, and the degree of asymmetry on the accuracy of the PBME and FB algorithms. © 2013 American Institute of Physics. [http://dx.doi.org/10.1063/1.4799272] I. INTRODUCTION

Nonadiabatic dynamics is a quantum phenomenon which occurs in systems that interact sufficiently with their environments to cause a breakdown of the Born-Oppenheimer approximation. It is involved in many chemical processes such as proton/electron transfer in solution and biological systems, and the response of molecules to radiation fields and their subsequent relaxation. The development of general theoretical methods for studying such phenomena in the condensed phase remains an open problem since full quantum mechanical simulations scale exponentially with the number of degrees of freedom (DOF) and, thus, become prohibitively expensive for condensed phase systems. The last several decades have witnessed the development of an array of methods for describing the quantum dynamics of these systems, including mean-field and surface-hopping schemes,1–7 semi-classical path integral methods,8–22 and algorithms based on the quantum-classical Liouville equation (QCLE).23–35 For systems of considerable size, those methods which treat a few crucial DOF (e.g., a) [email protected]

0021-9606/2013/138(14)/144106/12/$30.00

transferring protons and electrons, high-frequency vibrational modes, photoexcited electronic DOF of a chromophore, exciton states of a light harvesting system) fully quantum mechanically and their environment (e.g., remainder of chromophore, protein backbone, solvent) at a lower level of approximation, can provide computationally feasible and accurate algorithms for simulating nonadiabatic processes. For example, many surface-hopping schemes, mean-field methods,17, 36, 37 and some semi-classical path integral formulations adopt a mixed quantum-classical approach in which a few DOF of primary interest are evolved quantum mechanically, while the rest are evolved according to Newton’s equations of motion. In this work, we rely on treatments based on the QCLE. The QCLE has been shown to yield excellent agreement with numerically exact quantum results for a variety of model systems, attesting to its ability to accurately describe nonadiabatic dynamics and to properly account for quantum decoherence. A number of different algorithms, whose details depend on the basis chosen to represent the quantum subsystem, have been constructed for its solution.23, 24, 30, 33, 35, 38–41 When the QCLE is represented in the adiabatic basis, its solution can be obtained from an ensemble of surface-hopping

138, 144106-1

© 2013 American Institute of Physics

144106-2

Rekik et al.

trajectories, which incorporate quantum coherence effects in a rigorous fashion through evolution segments on the mean of two adiabatic surfaces. Unfortunately, in the case of complex systems, this approach becomes more demanding due to the increase in computational time required to average out the highly oscillatory terms which enter into the calculation of an observable.42 Inspired by some semi-classical approaches,17 Kim et al. used the mapping basis9, 11, 43 to represent the discrete quantum DOF in the QCLE in terms of continuous variables, giving rise to a trajectory description for all DOF in the system.33–35, 44 The Liouvillian in the resulting equation of motion is composed of a Poisson bracket term that acts in the phase space of both the quantum and classical DOF, and a more complex term responsible for higher order correlations between the subsystem and the bath. If this latter term is dropped, one arrives at a simple evolution equation for the coupled dynamics of the quantum and classical DOF, known as the Poisson Bracket Mapping Equation (PBME), which can be solved by propagating an ensemble of independent trajectories under Hamiltonian dynamics. The PBME solution has been shown to be accurate for some model systems,33, 44, 45 but for others it can lead to instabilities.34 With the mapping representation as a starting point, the coherent state basis46 has also been used to derive an approximate solution of the QCLE in terms of forward-backward (FB) propagators.35 This solution gives rise to a simple coupled dynamics in which the coherent state variables are propagated in forward and backward trajectories, and the bath variables are propagated using the average potential determined from these forward and backward trajectories. In this sense, the two sets of trajectories of the coherent state variables are connected by the evolution of the bath, which is consistent with the fact that the full solution of the QCLE in the mapping basis can be expressed in terms of an ensemble of entangled trajectories.34 This algorithm has been tested on a host of simple model systems including the two-level symmetric and asymmetric spin-boson models, exhibiting a range of agreements with the numerically exact quantum results from moderate to excellent.56 In light of the above, the mapping basis could provide an alternative route for deriving faster, more robust QCLEbased algorithms that do not involve surface-hopping, which would be particularly useful in the simulation of more complex systems. Accordingly, the main goal of this paper is to test the accuracy of the PBME and FB approaches for simulating the dynamics of quantum subsystems that interact with time-dependent fields and that are strongly coupled to their environments, which is relevant to the study of nonadiabatic dynamics in nonlinear spectroscopic experiments and in photoinduced processes such as excited-state proton transfer in chromophores and electron transfer reactions in photosynthetic reaction centres. More specifically, we apply these methods to the calculation of the time-dependent population of the donor state in a simple model for a condensed phase electron transfer reaction initiated by a laser pulse. The numerically exact quantum result for the population dynamics, obtained via the self-consistent hybrid (SCH) method47, 48 is reported in Ref. 49, and provides us with a benchmark for comparison. Since the full surface-hopping solution of the

J. Chem. Phys. 138, 144106 (2013)

QCLE is known to be more accurate than both PBME and FB, we also compute the donor-state population via the sequential short-time propagation algorithm for comparison with both the PBME/FB and exact results.38 The paper is organized as follows: Section II outlines the dynamical approaches taken for treating a mixed quantumclassical system driven by a classical electric field. In Sec. III, we describe the simple model for a photoinduced electron transfer reaction in the condensed phase. Section IV presents the results and discussion for the application of these quantum-classical methods to the model system. Concluding remarks are made in Sec. V. II. MIXED QUANTUM-CLASSICAL LIOUVILLE DYNAMICS OF A SYSTEM DRIVEN BY A CLASSICAL ELECTRIC FIELD

Consider a quantum subsystem driven by a classical electric field that is coupled to a classical environment, whose time-dependent Hamiltonian is given by P2 pˆ 2 ˆ + Vˆc (q, ˆ R) + Wˆ (q, ˆ t) Hˆ (t) = + Ve (R) + + Vˆs (q) 2M 2m ˆ p, ˆ R, P) + Wˆ (q, ˆ t) ≡ Hˆ M (q, ˆ p) ˆ + Vˆc (q, ˆ R) + Wˆ (q, ˆ t) ≡ He (R, P) + hˆ s (q, ˆ ≡ He (R, P) + h(R, t),

(1)

ˆ and pˆ are the vectors of masses, positions, and where m, q, momenta of the n photoactive DOF (i.e., the quantum subsystem), respectively; M, R, and P are the vectors of masses, positions, and momenta of the N photoinactive DOF (i.e., the classical environment), respectively; Hˆ M = Hˆ − Wˆ is the pˆ 2 P2 material Hamiltonian; He = 2M + Ve and hˆ s = 2m + Vˆs are the environment and subsystem Hamiltonians, respectively; Vˆc is the subsystem-environment coupling potential energy and hˆ = hˆ s + Vˆc + Wˆ [throughout this paper, operators are ˆ and, for the remainder of the paper, all veccapped (e.g., A) tors except those corresponding to masses, positions, and momenta will be boldfaced (e.g., A)]. The field-subsystem interˆ t), is assumed to be given by action term, Wˆ (q, ˆ t) = −μ( ˆ · E(t), Wˆ (q, ˆ q)

(2)

ˆ is the subsystem’s dipole moment operator, and where μ( ˆ q) the incident electric field is generally expressed in terms of a sum of Np coherent laser pulses E(t) =

Np 

Aα fα (t − tα ) exp{i[kα r − ωα (t − tα )]} + c.c.,

α=1

(3) where each laser pulse is characterized by its wave vector kα , leading frequency ωα , and pulse envelope fα (t − tα ), and Aα is the product of the polarization vector and the overall amplitude of the pulse. The envelope function is taken to be a Gaussian centered at time tα ,    4 ln 2 (t − tα )2 , (4) exp −4 ln 2 fα (t − tα ) = π τα2 τα2 where τ α is the full-width at half-maximum of the pulse.

144106-3

Rekik et al.

J. Chem. Phys. 138, 144106 (2013)

In general, we are interested in the quantum-classical dynamics of an observable, which can be described by the partial ˆ over Wigner transform50 of its corresponding operator, A(t), the environmental DOF,  ˆ ˆ AW (X, t) = dZeiP ·Z/¯ R − Z/2|A(t)|R + Z/2, (5) where X = (R, P). Within this description, the expectation value of Aˆ W is given by  ˆ A(t) = T rs dX ρˆW (X)Aˆ W (X, t), (6) where Trs denotes a trace over the subsystem DOF and the partial Wigner transform of the density operator, ρ, ˆ is given by    1 N ˆ + Z/2. dZeiP ·Z/¯ R − Z/2|ρ|R ρˆW (X) = 2π ¯ (7) When m/M  1, the QCLE can be used to accurately describe √ the dynamics of Aˆ W (X, t) to first order in m/M,28, 29, 32 ∂ ˆ i AW (X, t) = [Hˆ W (X, t), Aˆ W (X, t)] ∂t ¯ 1 − ({Hˆ W (X, t), Aˆ W (X, t)} 2 − {Aˆ W (X, t), Hˆ W (X, t)}).

¯  ∂hλλ 8 λλ ∂R   ∂ ∂ ∂ ∂ ∂ + × Am (t). ∂rλ ∂rλ ∂pλ ∂pλ ∂P

≡ − {Hm (t), Am (t)}X,x +

(10)

In this equation, the mapping analogue of an operator Aˆ W (X) is given by 

Aλλ (11) Am (x, X) = W (X)cλλ (x), λλ

where x = (r, p) = (r1 , . . . , rn , p1 , . . . , pn ) denotes the positions and momenta of the subsystem mapping variables and the subscripts/superscripts λλ denote a representation in the subsystem basis, {|λ; λ = 1, . . . , n}, which is defined by the eigenvalue problem hˆ s |λ = λ |λ. Here, cλλ (x) is given by44 1 [rλ rλ + pλ pλ + i(rλ pλ − rλ pλ ) − ¯δλλ ]. 2¯ (12) Consequently, the mapping Hamiltonian may be written as 1  λλ

h (R, t) Hm (x, X) = He (X) + 2¯ λλ

cλλ (x) =

× (rλ rλ + pλ pλ − ¯δλλ ),

(13)

where

ˆ + Vˆc (q, ˆ R) + Wˆ (q, ˆ t)|λ  hλλ (R, t) = λ|pˆ 2 /2m + Vˆs (q)

(8)

Here, [···] is the commutator and {···} is the Poisson bracket given by {Aˆ W (X, t), Bˆ W (X, t)} → → ← − − ← − − = Aˆ W (X, t)( ∇ P ∇ R − ∇ R ∇ P )Bˆ W (X, t) ≡ Aˆ W (X, t)Bˆ W (X, t), (9) ← − − → → ← − − − → ← − where  = ( ∇ P ∇ R − ∇ R ∇ P ), and ∇ R/P and ∇ R/P correspond to taking the gradient with respect to R/P of the term to the right and left, respectively. A. Quantum-classical Liouville dynamics in the mapping basis

In this section, we provide a summary of the main theoretical results needed for solving the QCLE and calculating expectation values of operators in the mapping basis. More detailed discussions are provided in Refs. 33 and 44. In the mapping basis, the QCLE takes the following form:33 ∂ Am (X, x, t) ∂t   1  λλ

∂ ∂ = pλ Am (t) h − rλ ¯ λλ

∂rλ

∂pλ

  ∂Hm ∂ P ∂ − Am (t) + M ∂R ∂R ∂P 

 ∂ ∂ ∂ ¯  ∂hλλ ∂ ∂ + + Am (t) 8 λλ ∂R ∂rλ ∂rλ ∂pλ ∂pλ ∂P





= λ δλλ + Vcλλ (R) + W λλ (t),

(14)



with W λλ (t) = −μˆ λλ · E(t) and it has been assumed that



hλλ = hλ λ . It should be noted that the Poisson bracket, {Am (t), Bm (t)}X, x , acts in the space of both the classical bath phase space and quantum mapping variables. If the last term of Eq. (10) is dropped, Am (t) evolves according to a simple set of Hamiltonian equations of motion33 ∂Hm drλ 1 = = hλλ (R(t), t)pλ (t) dt ∂pλ ¯ λ

dpλ ∂Hm 1 =− =− hλλ (R(t), t)rλ (t) dt ∂rλ ¯ λ

(15)

dR(t) ∂Hm P (t) dP (t) ∂H = = , =− . dt ∂P (t) M dt ∂R(t) Thus, the coupled dynamics of the mapping and bath variables can be easily simulated in terms of Newtonian trajectories. An algorithm for numerically solving this coupled set of differential equations is derived in the Appendix. This approximate form of Eq. (10) (containing only the Poisson bracket term) is known as the PBME. On the other hand, if the complicated last term of Eq. (10) (containing products of bath and mapping derivatives) is retained, a more complex algorithm would be required for solving this equation. The results of previous PBME simulations of the spin-boson model,33 simple curve crossing models,44 and the Fenna-Mathews-Olsen complex,45 agree well with the exact results, indicating that the effect of the last term is small for these systems. However, in the case of simple conical intersection and collinear reactive collision models, the PBME results do not agree as

144106-4

Rekik et al.

J. Chem. Phys. 138, 144106 (2013)

well with the full QCLE and exact results.34 In Ref. 44, it is argued that the effects of this extra term may be small if the subsystem-bath coupling is weak or if the bath is large and only a part of it is directly coupled to the subsystem. Thus, the validity of the PBME varies from system to system. ˆ may Finally, the expectation value of an observable A(t) 33 be written in terms of the mapping basis as  A(t) = dxdXAm (x, X, t)ρ˜m (x, X), (16) where ρ˜m (x, X) =

1  λ λ gλλ (x)ρW (X). (2π ¯)n λλ

(17)

and applying the effective evolution operator to the coherent state and bath phase space variables leads to   d 2 z1 d 2 z

1 (X, t) = mλ |z1 z1 |mλ  Aλλ W n πn π μμ

 μμ

× z1 (t1 )|mμ AW (Xt )mμ |z1 (t) .

Evaluating the overlaps according to mλ |z = zλ e−|z| /2 , expressing the above equation in terms of x = (q, p) variables, and using the fact that ν (qν2 + pν2 ) is conserved under coherent state dynamics, leads to35 

dxdx φ(x)φ(x ) Aλλ (X, t) = W 2

μμ

In the above equation, gλλ (x) is given by44 2n+1 −x 2 /¯ gλλ (x) = e ¯   ¯ × rλ rλ − i(rλ pλ − rλ pλ ) + pλ pλ − δλλ . 2 (18)

In this section, we provide a summary of the main theoretical results needed for the FB trajectory solution of the QCLE in terms of the coherent state basis. Using the coherent state basis, the matrix elements of Aˆ W (X, t) in the subsystem basis may be expressed as M  d 2 zi d 2 zi

Aλλ (X, t) = mλ |z1 z1 |mλ  W n πn π μμ

i=1



× z1 (t1 )|z2 ei t2 Le (Xt1 ,z2 ,z2 )/2 z2 | . . .   μμ

×AW (Xt1 ) . . . |z2  z2 |z1 (t1 ) , (19)

where |mλ  is the eigenstate of n fictitious harmonic oscillators (with occupation numbers 0 or 1 such that |mλ  = |01 , . . . , coherent state46 at 1λ , . . . , 0n ),17, 51 |zi  is the n-dimensional √ time step i with eigenvalue z = (q + ip)/ 2¯, and ti = ti − ti−1 = τ for all i with t0 = 0 and tM = t. Here, q = (q1 , . . . , qn ) and p = (p1 , . . . , pn ) are the average positions and momenta of the harmonic oscillators in the state |z, respecˆ tively; e.g., q = z|q|z. In the above equation, the effective operator iLe (X, z, z ) is given by ∂ ∂Ve (X, z, z ) ∂ P · − · , M ∂R ∂R ∂P

1 μμ

(qλ + ipλ )(qλ − ipλ )AW (Xt ) 2¯ 1 × (qμ (t) − ipμ (t))(qμ (t) + ipμ (t)), 2¯

×



(22)

where φ(x) = (2π ¯)−n e− ν (qν +pν )/2¯ . The equations of motion governing the coupled evolution of the bath and coherent state variables in this equation are 2

2

∂Hcl (R, P , q, p) dpμ ∂Hcl (R, P , q, p) dqμ = =− , , dt ∂pμ dt ∂qμ

B. Quantum-classical Liouville dynamics in the coherent state basis

iLe (X, z, z ) =

(21)

(20)

where Ve (X, z, z ) = [Vcl (R, z) + Vcl (R, z )]/2 and Vcl (R, z)



ˆ + Vcλλ (R)zλ∗ zλ . As a result, eiLe (X,z,z )τ = Ve (R) − Trs h(R) Aˆ W (X) = Aˆ W (Xτ ). If the phase space coordinates of two coherent states are substantially different, then the overlap between them decays quickly. In such a case, it is assumed that z1 (t1 )|z2  ≈ π n δ(z2 − z1 (t1 )) and z2 |z1 (t1 ) ≈ π n δ(z2 − z1 (t1 )), for example. Then, performing the integrals over the coherent state variables zi and zi for {i, i } ≥ 2

∂Hcl (R, P , q , p ) dpμ

∂Hcl (R, P , q , p ) dqμ

= =− , , dt ∂pμ

dt ∂qμ

P dP ∂He (R, P , q, p, q , p ) dR = , =− , dt M dt ∂R where

(23)

He (R, P , q, p, q , p ) =

1 [Hcl (R, P , q, p) + Hcl (R, P , q , p )] 2

(24)

with ˆ Hcl (X, x) = He (X) − Trs h(R) 1  λλ

+ h (qλ qλ + pλ pλ ). 2¯ λλ

(25)

Solving these equations gives rise to a simple dynamics in which the forward and backward trajectories of the coherent state variables are propagated forward in time, while the bath coordinates evolve under the influence of the mean potential that depends on these forward and backward trajectories. It should be noted that the two sets of coherent state variables are only coupled via the bath coordinates. ˆ Finally, the expectation value of an observable A(t) is given by 

λ λ dXAλλ A(t) = (X). (26) W (X, t)ρ λλ

C. Quantum-classical Liouville dynamics in the adiabatic state basis

The subsystem DOF may be represented using the states of an adiabatic basis, |α;R, which are the solutions of

Rekik et al.

144106-5

J. Chem. Phys. 138, 144106 (2013)

ˆ h(R)|α; R = Eα (R)|α; R. In this basis, the QCLE takes the following form:28  ∂Aαα ββ

W (X, t) =i Lαα ,ββ AW (X, t), ∂t ββ



∂ 1

∂ P

 · + FWα + FWα · , M ∂R 2 ∂P

(28)

(29)

which, when the quantum indices are different, yields classical propagation via Hellmann-Feynman forces, FWα ˆ (q,R) ˆ = −α; R| ∂ Vc∂R |α; R, on mean surfaces accompanied by quantum mechanical phase oscillations with frequency ωαα = (Eα − Eα )/¯. The diagonal terms correspond to adiabatic evolution on single potential energy surfaces. The second term   P ∂ 1 Jαα ,ββ (t) = − · dαβ 1 + Sαβ · δα β

M 2 ∂P   1 ∗ P ∂ ∗ · d 1 + Sα β · δαβ − M αβ 2 ∂P i



− [W αβ (t)δα β − W β α (t)δαβ ] ¯

(30)

is responsible for nonadiabatic transitions and associated changes in the momentum of the bath in order to conserve enαβ P ergy. The quantity Sαβ is defined as Sαβ = FWα δαβ − FW ( M · αβ P dαβ )−1 = Eαβ dαβ ( M · dαβ )−1 , where FW are the off-diagonal matrix elements of the force and dαβ is the nonadiabatic cou∂ |β; R. The nonadiapling matrix element, dαβ = α; R| ∂R batic coupling matrix elements in this operator account for nonradiative nonadiabatic transitions, while the bath momentum derivative accounts for the energy transfer involved in the subsystem state change. The matrix element of the fieldmatter interaction term, W αβ (t) = −μαβ (R) · E(t), is responsible for radiative nonadiabatic transitions. The solution of the QCLE in the adiabatic basis may be obtained via a decomposition of the propagator into short-time segments.38 The quantum-classical evolution of

Aαα W (X, t) can be written as 

α αN

(ei Lt )αα ,αN αN AWN ˆ

(X) .

(31)

αN αN

If the time interval t is divided into N segments such that the j th segment has length tj = tj − tj − 1 = t (which may be chosen to be either equal to or an integer multiple of the time step), we have

N  ˆ 

α α

Aαα (ei L tj )αj −1 αj −1 ,αj αj AWN N (X). W (X, t) = (α1 α1 )...(αN αN )

iLαj −1 α

j −1

tj

× (δαj −1 αj δαj −1 αj + tJαj −1 αj −1 ,αj αj ), (33)

The first term contains the classical Liouville operator,



ˆ

iωαj −1 α

iLαα ,ββ = (iωαα + iLαα )δαβ δα β − Jαα ,ββ .

Aαα W (X, t) =

(ei L tj )αj −1 αj −1 ,αj αj ≈ Wαj −1 αj −1 (tj −1 , tj )e

(27)

where the QCL evolution operator is

iLαα =

The short-time propagator can be written as

j =1

(32)

tj

j −1 where Wαj −1 αj −1 (tj −1 , tj ) = e is the phase factor associated with that time segment. The algorithm prescribed by Eq. (32) corresponds to a surface-hopping solution of the QCLE. The short-time segments involve evolution along the surface (αα ), which may be adiabatic (when α = α ) or the arithmetic mean of two adiabatic surfaces (when α = α ), governed by the propagator e−i(ωαα +Lαα )t . These trajectory segments are interrupted by nonadiabatic transitions due to the operator J , which causes the system to undergo a transition to a new surface (or mean surface), followed by subsequent evolution on this surface. It should be noted that the algorithm described in Ref. 38 was slightly modified in order to account for the nonadiabatic transitions due to the external field. Within this approach, the expectation value of an observˆ is given by able A(t) 

α α A(t) = (34) dXAαα W (X, t)ρW (X).

αα

III. MODEL OF PHOTOINDUCED ELECTRON TRANSFER IN THE CONDENSED PHASE

In this work, we apply the above techniques for simulating the field-driven population dynamics in a model developed in Ref. 49 for a photoinduced electron transfer (ET) reaction in the condensed phase. In Ref. 49, the authors computed the pump-probe signal for this model based on a nonperturbative approach for describing nonlinear optical response and a SCH method47, 48 for simulating the quantum dynamics of the system, which is in principle numerically exact. The model subsystem is composed of three electronic states: a ground state |g, a photoinduced excited state |d corresponding to the donor of the ET reaction (prepared by an excitation from the ground state by a laser field), and an optically dark charge transfer state |a corresponding to the acceptor of the ET reaction. The model bath (used to mimic the effect of a polar solvent on the ET reaction) is composed of a collection of independent harmonic oscillators that are linearly coupled only to the acceptor state of the subsystem, rendering this model asymmetric. In our mixed quantum-classical approach, the subsystem is treated quantum mechanically, whereas the harmonic oscillator bath is treated classically. The Hamiltonian for this model system is given by Hˆ (t) = |g g g| + |d d d| + |a a a|+ (|da|+|ad|) ⎡ 2 ⎤  N  1 ⎣Pj2 + ωj2 Rj + |a 2cj a| ⎦ − μ·E(t), ˆ + 2 j =1 ωj2 (35) where g , d , and a are the energies of the ground, excited donor, and excited acceptor electronic states, respectively,

144106-6

Rekik et al.

J. Chem. Phys. 138, 144106 (2013)

is the donor-acceptor electronic coupling, Rj and Pj are the mass-weighted coordinate and momentum, respectively, of the j th bath mode with frequency ωj and coupling constant cj , and Np = 1 in E(t). The subsystem dipole moment operator, which couples states |g and |d, is given by μˆ = |gμgd d| + |dμdg g|,

(36)

where μgd = μdg is the electronic transition dipole moment. It should be noted that in this model both and μgd are assumed to be independent of the classical DOF. From Eq. (13), the mapping form of the Hamiltonian is given by Hm (t) =

N   g 2 1  2 rg + pg2 − ¯ Pj + ωj2 Rj2 + 2 j =1 2¯

β¯ω

where ZM = Tr(e−β HM ) is the partition function of the material system and   2 ¯ −x 2 /¯ 2 2 e (44) ρ˜sm (x) = 4 3 rg + pg − ¯ π 2 with x 2 = x · x = λ (rλ2 + pλ2 ). The expectation value in Eq. (43) is determined by sampling the initial bath and mapping variables from their appropriate Gaussian distributions, reweighting by ¯42π 3 (rg2 + pg2 − ¯2 ), and simulating the time evolutions of rd (t) and pd (t) according to Eq. (15) with the mapping Hamiltonian in Eq. (37). According to the FB approach, the time-dependent population of the donor state, Pd (t), can be computed using Eqs. (22) and (26),  1 Pd (t − t1 ) = 2 dXdxdx ρe (X)φ(x)φ(x ) 4¯ ˆ

(37)

j

(38)

where He = 12 j (Pj2 + ωj2 Rj2 ) and I is the identity matrix. The linear coupling between the subsystem and environment is characterized by its spectral density J (ω) c2 = π2 j ωjj δ(ω − ωj ), which is chosen to be Debye in form, namely,52 (39)

where λD is the bath reorganization energy and ωD is the characteristic frequency. Following the bath discretization procedureoutlined in Refs. 47 and 48, one can show that cj = λD tan−1 (ωmax /ωD )/(π N)ωj and ωj = tan (j tan −1 (ωmax /ωD )/N)ωD , where ωmax is the maximum frequency of the spectral density and N is the number of oscillators.53 The initial density matrix is assumed to be uncorrelated, with the subsystem in its ground state and the bath in thermal equilibrium ρW (0) = ρs (0)ρe (X).

with u

j = uj coth uj and uj = 2 j . With the aid of Eqs. (16)–(18), we can compute the timedependent population of the donor state, Pd (t), via the PBME approach according to

(43)

+ He I,

λD  ωωD J (ω) = , 2 2 j ω2 + ωD

and the Wigner distribution of the bath, ρ e (X), is given by54

 2  N βωj β Pj 1 2 2 ρe (X) = exp −

(42) + ωj Rj 2π u

j uj 2 2 j =1

Pd (t − t1 ) =

In the subsystem basis {|g, |d, |a}, Hˆ (t) has the following matrix form: ⎛ ⎞

g −μgd E (t) 0 ⎜ ⎟

d Hˆ (t) = ⎝−μdg E(t) cj2 ⎠ 0

a +2 j Rj cj +2 j ω2 ˆ ≡ h(R) + He I,

(41)

1 ˆ Tr[e−β HM Uˆ † (t)|dd|Uˆ (t)] ZM  = dxdXρe (X)ρ˜sm (x)(rd (t)2 + pd (t)2 − ¯),



d 2 rd + pd2 − ¯ + 2¯ ⎞ ⎛   2cj2 1 ⎝ ⎠ + 2Rj cj +

a + 2 2¯ ω j j j

 × ra2 + pa2 − ¯ + (rd ra + pd pa ) ¯ 1 − [(μgd + μdg )E(t)(rg rd + pg pd )]. ¯

Here, the subsystem density matrix, ρ s (0), is ⎛ ⎞ 1 0 0 ρs (0) = ⎝ 0 0 0 ⎠ , 0 0 0

(40)

× (qg + ipg )(qg − ipg )(qd (t) − ipd (t))(qd (t) + ipd (t)).

(45)

This expectation value is determined by sampling the initial bath and mapping variables from their appropriate Gaussian distributions, reweighting by (qg + ipg )(qg − ipg ), and simulating the time evolutions of the coherent state variables qd (t), pd (t), qd (t), and pd (t) according to Eq. (23) with the Hamiltonian in Eq. (25) and the definitions of He and hˆ given above. Using the surface-hopping approach, the time-dependent population of the donor state, Pd (t), can be computed using Eqs. (32) and (34),  Pd (t − t1 ) = dXρe (X)cαd (t)cαd (t), (46) where cαd (t) = α; R|d(t). This expectation value is determined by sampling the initial bath variables from its Gaussian distribution and simulating the time evolutions of the cαd (t)’s according to Eq. (32).

144106-7

Rekik et al.

J. Chem. Phys. 138, 144106 (2013)

IV. RESULTS AND DISCUSSION

Population / TA signal

We now present the results of the population dynamics of the donor state in the model electron transfer reaction described above. The parameters used in our simulations were taken from Ref. 49. The electronic energies of the ground, excited donor, and excited acceptor states are g = 0 cm−1 , d = 13 000 cm−1 , and a = 13 000 cm−1 , respectively, and the donor-acceptor coupling is = 50 cm−1 . The characteristic frequency and reorganization energy of the bath, dictated by the Debye spectral density in Eq. (39), are ωD = 50 cm−1 and λD = 500 cm−1 , respectively. The number of harmonic oscillators needed to represent the condensed phase environment is N = 20 and the temperature is T = 300 K. The laser pump pulse has a carrier frequency of ω = 13 000 cm−1 (in resonance with the electronictransition), a duration of τ = 50 fs, and a strength of μgd |A| 4 ln 2/(π τ 2 ) = 50 cm−1 . For convenience, we have used dimensionless variables and parameters with time scaled by ωD in all of our simulations. The dimensionless time step used for integrating the equations of motion (according to the algorithm in the Appendix) is t = 2.5 × 10−5 . Converged results were obtained with 2 × 105 trajectories in the case of PBME and FB and 3 × 107 trajectories in the case of surface-hopping. Figure 1 contains the main results of this paper for the time-dependent population of the photoexcited donor state, Pd (t). We first consider the population dynamics when = 0 cm−1 (upper panel of Fig. 1), which corresponds to no electron transfer and, in turn, no subsystem-bath coupling. In this

0.06 0.05 0.04 0.03

TA signal−Exact Population−PBME Population−FB Population−SH

0.02 0.01 0 0.07

Population, Pd

0.06 0.05 0.04 0.03

Exact, Δ=50 cm−1 PBME, Δ=50 cm−1 FB, Δ=50 cm−1 SH, Δ=50 cm−1 AD, Δ=50 cm−1

0.02 0.01 0 −100

0

100 200 300 400 500 600 700 800

Time, t (fs) FIG. 1. (Upper panel) Comparison between the transient absorption (TA) signal from Ref. 49 and the PBME (red), FB (green), and surface-hopping (blue) calculations of the time-dependent population of the donor state, Pd (t), in the absence of electron transfer (i.e., = 0 cm−1 ). (Lower panel) The PBME (red), FB (green), surface-hopping (blue), adiabatic (pink), and exact (black) time-dependent populations of the donor state, Pd (t), in the presence of electron transfer (i.e., = 50 cm−1 ).

case, PBME, FB, and surface-hopping yield identical results with the population increasing to a plateau value of ∼6.2% due to photoexcitation by the pump pulse. For comparison, the exact transient absorption (TA) signal for this model, computed in Ref. 49, is shown since it monitors the population dynamics.60 It should be noted that the magnitude of the TA signal was rescaled by the authors to coincide with that of the population of the photoexcited donor state. From Fig. 1, we see that the QCLE-based methods and the exact result have very similar profiles and plateau at the same value. This is to be expected since when the subsystem is decoupled from the bath, the excess coupling term in the full mapping QCLE does not contribute, the delta function approximation made in the derivation of the FB solution becomes exact, and the surfacehopping solution is exact. When = 50 cm−1 , electron transfer occurs and, in turn, the subsystem and bath couple. In this case, based on the exact result from Ref. 49, the population of the donor state should decay monotonically to the acceptor state after the initial rapid increase due to the pump pulse (see lower panel of Fig. 1). However, in the case of PBME, the population continues to increase after the initial photoexcitation. An examination of the dynamics of the other states (not shown) reveals that this nonphysical increase is due to the population of the acceptor state becoming negative. Thus, under these conditions, PBME is not capable of qualitatively capturing the population dynamics. In an effort to improve this result, we also employed the traceless form of the Hamiltonian according to Ref. 34 (which, if used, has been shown to tame dynamical instabilities associated with approximate solutions of the QCLE), but saw that the result actually becomes worse, exhibiting a more pronounced increase in the population (result not shown). On the other hand, FB is capable of qualitatively capturing the decay, but only by ∼0.4% (cf. ∼2% in the case of the exact result). It should be noted that both PBME and FB are not capable of capturing the decrease in the population maximum at t ≈ 60 fs (caused by population transfer to the acceptor state, while the pulse is on) in going from = 0 cm−1 to = 50 cm−1 . Nevertheless, these results suggest that replacing the overlap integral by a delta function (see text above Eq. (21)) to arrive at the FB solution is less harsh of an approximation than dropping the extra coupling term in the QCLE to arrive at the PBME. In the case of surface-hopping, not only does it capture this decrease in the population maximum, but it also captures the 2% decay. For comparison, we also computed the adiabatic result, which was obtained by prohibiting nonradiative nonadiabatic transitions, in order to give us a sense of the magnitude of the nonadiabatic effects. As can be seen, the adiabatic and surface-hopping results are in good agreement until the end of the pulse, but begin to diverge after that, indicating that nonradiative nonadiabatic transitions are crucial for accurately capturing the population transfer to the acceptor state. At the end of this section, we discuss potential reasons for the failure of PBME and the inability of FB to quantitatively capture the decay. Previous calculations on symmetric and asymmetric spin-boson models have shown that both PBME and FB perform well in the case of the symmetric models, but less so in the case of the asymmetric models. The 2 j cj2 /ωj2 |aa|

Rekik et al.

J. Chem. Phys. 138, 144106 (2013)

0.06

0.06

0.05

0.05

Population, Pd

Population, Pd

144106-8

0.04 0.03 0.02 0.01 0 −100

0.04 0.03 0.02 0.01

0

100 200 300 400 500 600 700 800

Time, t (fs)

0 −100

0

100 200 300 400 500 600 700 800

Time, t (fs)

FIG. 2. The time-dependent population of the donor state, Pd (t), generated without the initial energy shift between the donor and acceptor states, calculated via the PBME, FB, and surface-hopping approaches. The colour scheme is the same as in Fig. 1.

FIG. 3. The time-dependent population of the donor state, Pd (t), for two values of the bath reorganization energy [λD = 50 cm−1 (solid) and λD = 250 cm−1 (dashed)], calculated via the PBME, FB, and surface-hopping approaches. The colour scheme is the same as in Fig. 1.

term in the Hamiltonian introduces an asymmetry of ≈10 between the donor and acceptor states at t = 0 fs, and thus it is not surprising that these methods do not yield accurate results. This may be due to the fact that, in the limit of strong subsystem-bath coupling, the magnitude of this term can be large relative to that of the 2 j Rj cj |aa| term (which is initially zero, on average) during the ultrafast pulse, thereby making it more difficult to correctly capture the stabilization of the acceptor state by the bath, which leads to electron transfer. However, if one drops this term, the resulting Hamiltonian becomes symmetric initially (on average), but as time goes on the asymmetry builds due to the fact that the bath is only coupled to the acceptor state. The results for this initially symmetric Hamiltonian are presented in Fig. 2, where we see that after the initial photoexcitation the populations decay by ∼2.4% and ∼1% in the cases of PBME and FB, respectively (cf. ∼−1% and ∼0.4%, respectively, in the initially asymmetric case). These results are more consistent with the surfacehopping result, which exhibits a ∼2.8% population decay (cf. ∼2% decay in the initially asymmetric case). In Fig. 3, we show Pd (t) for three values of the bath reorganization energy, λD = 50 cm−1 , 250 cm−1 , and 500 cm−1 . As the bath reorganization energy decreases, the electron transfer to the acceptor state occurs more readily and, consequently, we would expect a more pronounced decay in the population of the donor state. In the case of PBME, a reduction of λD from 500 cm−1 (or λD β = 2.4) to 250 cm−1 (or λD β = 1.2) leads to a smaller increase in the donor-state population after the photoexcitation is over. When λD is further reduced to 50 cm−1 (or λD β = 0.24), the population decays (as expected) to ∼2.6%. In the case of FB, as expected, the decay becomes more pronounced as the reorganization energy is reduced. When compared to the surface-hopping results, we see that the PBME and FB results improve as λD is reduced, with decent agreement between the three methods found for λD = 50 cm−1 . These results suggest that FB performs better than PBME over a wider range of subsystem-bath couplings for this particular model. In several previous studies of standard model systems, the use of focused sampling18, 55 of the initial conditions of the quantum subsystem variables has substantially re-

duced the number of trajectories required for numerical convergence.18, 33, 55 However, in two recent studies,55, 56 it was explained that focused sampling is a good approximation when the model is either weakly nonadiabatic in nature or similar to a symmetric spin-boson model where the subsystem states are degenerate in energy. Unfortunately, both requirements are not satisfied in our case. First of all, the radiative transition is a highly nonadiabatic event and, second, due to the energy shift and the fact that only the acceptor state couples to the bath, the system resembles a generalized asymmetric spin-boson model. Since these sufficient conditions are not satisfied in our model, we have no theoretical basis to expect the results with focused sampling to be good. As will be shown below, the FB result is reasonably good, whereas the PBME result is not. In Fig. 4, we present the donor-state population results obtained by setting the initial √ values of the √ mapping and coherent state variables to {rg = 0.75, pg = 0.75, rd = 0, pd = 0, ra = 0, pa = 0} and {qg = 1, pg = 1, qd = 0, pd = 0, qa = 0, pa = 0, qg

= 1, pg = 1, qd = 0, pd = 0, qa = 0, pa = 0}, respectively. As can be seen, the use of focused initial conditions yields unphysical results in the case of PBME, with the population starting at ∼66.6% and ending at ∼60.9% (see upper panel), whereas in the case of FB, the agreement with the original result is fairly good (see lower panel). The drastic difference between the PBME results for focused and Gaussian sampling is a consequence of the form of the observable, namely, the −¯ term in Eq. (43). Due to this term, at t = 0, the normalized population (in dimensionless units) for each trajectory will, according to Eq. (43), always be (1)(−1) Pd (0) = ≈ 0.67. Pg (0) + Pd (0) + Pa (0) (1)(0.5 − 1 − 1)

(47)

In the above equation, the −¯ term prevents the unnormalized donor and acceptor populations from being positive definite. It should be noted that for population differences (e.g., Pd (t) − Pg (t)), the −¯ term cancels and accurate PBME results can be obtained in some cases.33 In the case of FB, the observable (see Eq. (45)) does not contain such a −¯ term, and therefore physically reasonable results can be obtained when focused initial conditions are used.

144106-9

Rekik et al.

J. Chem. Phys. 138, 144106 (2013)

0.7

Population, Pd

0.6 Focused sampling Gaussian sampling

0.5 0.4 0.3 0.2 0.1 0 −100

0

100 200 300 400 500 600 700 800

0.07

Population, Pd

0.06 0.05 Focused sampling Gaussian sampling

0.04 0.03 0.02 0.01 0 −100

0

100 200 300 400 500 600 700 800

Time, t (fs) FIG. 4. (Upper panel) The effects of focused versus Gaussian sampling on the time-dependent populations of the donor state, Pd (t), in the presence of electron transfer (i.e., = 50 cm−1 ), calculated via the PBME approach. (Lower panel) The effects of focused versus Gaussian sampling on the timedependent populations of the donor state, Pd (t), in the presence of electron transfer (i.e., = 50 cm−1 ), calculated via the FB approach.

The failure of PBME and the inability of FB to quantitatively capture the population decay for this model (see lower panel of Fig. 1) is in part due to an interplay between two factors: (1) a fast, radiative population transfer from |g to |d, and (2) a large initial energy gap (i.e., on average, Eda = Hdd − Haa ≈ 10 ) between the donor and acceptor states. By comparing to the more accurate surface-hopping calculations, we can qualitatively understand the effects of these factors in the following manner. The population dynamics during the relatively short pulse period, as captured by PBME and FB, essentially solely involves transfer between the ground and donor states, since Eda shifts the acceptor state out of resonance with the donor state and thereby prevents subsequent population transfer to the acceptor state on the timescale of the pulse. Thus, at the end of the pulse period, PBME and FB predict populations localized in the ground and donor states. On the other hand, the surface hopping results reveal that the radiative transitions lead mostly to population transfer between the ground and first-excited adiabatic states. Due to Eda being large on average, the ground and first-excited adiabatic states are essentially orthogonal to the acceptor state during the initial stage of the pulse period. Up to this point, the PBME and FB results are in good agreement with the surfacehopping and exact results. However, as observed in individual surface hopping trajectories, at later points during the pulse period, the bath is occasionally able to sufficiently stabilize the energy of the acceptor state to cause avoided crossings

between the first-excited and second-excited adiabatic states. At an avoided crossing, the surface-hopping algorithm decides whether or not the system hops. If the system remains in the first-excited adiabatic state, there is an abrupt change in the first-excited adiabatic state from being mostly donor in character to being mostly acceptor in character, causing the acceptor-state population to increase. The existence of these avoided crossings leads to the decrease in the maximum of the donor-state population in the = 50 cm−1 surface-hopping result in Fig. 1 compared to the = 0 cm−1 result, as noted above. Since the surface-hopping result is in good agreement with the adiabatic result up to the end of the pulse period, this suggests that, up to this point, the dynamics is primarily adiabatic with the exception of radiative transitions. PBME and FB, which employ generalized mean potentials that depend on all of the quantum mapping/coherent-state variables to propagate the bath coordinates, do not accurately capture these abrupt changes in the potential surface during the short pulse period. After the pulse, in the case of surface hopping, the bath coordinates either evolve on the first-excited adiabatic surface or hop between the first-excited, mean of the first- and second-excited, and second-excited adiabatic surfaces. Upon averaging over an ensemble of such trajectories, we see a decrease in the donor-state population and an increase in the acceptor-state population. As mentioned earlier, these nonadiabatic transitions are required to fully capture the extent of the decay of the donor-state population. Therefore, the divergence between the PBME/FB and surface-hopping results as time goes by is partially due to the inadequacy of the mean potentials provided by these methods to yield an accurate picture of the nonadiabatic dynamics. Other potential reasons for the breakdown of the PBME and FB methods may be related to the form of the observable and the nature of the initial sampling. A detailed study aimed at comprehensively understanding how these various factors contribute to the breakdown would be worthwhile.

V. CONCLUDING REMARKS

The PBME and FB trajectory solutions, two approximate solutions of the QCLE in the mapping representation, yield simple coupled dynamics of the quantum and classical DOF in terms of continuous variables, thereby making them potentially viable alternatives to standard surface-hopping solutions. In order to determine the extent to which these approximations are generally valid, these algorithms must first be tested on model systems containing multiple quantum states, time-dependent Hamiltonians, strong subsystem-bath couplings, and asymmetries between states. Thus, in this paper, we apply these approximate QCLE-based techniques to the simulation of the laser-induced population dynamics of a photoexcited donor state in a three-state model of an electron transfer complex that is coupled asymmetrically to its condensed phase environment through the acceptor state, and has an initial energy shift, Eda , between the donor and acceptor states that depends on the strength of the subsystem-bath coupling.

144106-10

Rekik et al.

Our results show that, under conditions of strong subsystem-bath coupling (which, in turn, leads to a large initial  Eda  value), the PBME approach is not able to capture the physically expected decay of the donor-state population, whereas the FB approach does but not to the full extent as the exact result. The PBME and FB results are also benchmarked against the full QCLE surface-hopping solution, which is in good agreement with the exact result. By comparing with the adiabatic result, we are able to ascertain that both adiabatic dynamics and nonadiabatic transitions play important roles in reproducing the exact result, both of which are not adequately captured by the mean potentials inherent to PBME and FB. Under very weak subsystem-bath coupling conditions (e.g., λD β = 0.24), the PBME and FB solutions yield population decays that are similar to the more accurate surface-hopping result. This can be traced back to the fact that their underlying approximations improve as the strength of the subsystem-bath coupling decreases. Overall, we found that FB performs better than PBME over a wider range of subsystem-bath couplings, but, in general, both methods work best in low subsystembath coupling situations (i.e., λD β ≤ 1). For this particular model, decreasing the bath reorganization energy not only lowers the subsystem-bath coupling strength, but also reduces the initial asymmetry between the donor and acceptor levels (i.e., initial value of  Eda ). When the initial asymmetry is removed by omitting the 2 j cj2 /ωj2 |aa| term altogether from the Hamiltonian, the PBME and FB solutions yield population decays that are more consistent with the surfacehopping result. In both situations (i.e., weak subsystem-bath coupling and  Eda  = 0 initially), the improvement in the results is likely due to a minimization of the initial asymmetry. Given their attractive features (i.e., ease of implementation, good convergence properties), the PBME and FB algorithms may be useful in nonperturbative simulations of multidimensional optical spectra of model and realistic mixed quantumclassical systems that exhibit weak chromophore-bath coupling strengths. Such applications are currently underway in our group and will be reported in future publications.

ACKNOWLEDGMENTS

This work was supported by grants from the University of Alberta and the Natural Sciences and Engineering Research Council of Canada (NSERC). We would like to thank H. Wang and M. Thoss for helpful discussions.

J. Chem. Phys. 138, 144106 (2013)

formal solution of this equation is given by  t  dt iL0m (t ) Am (0) ≡ Um (t, 0) Am (0), Am (t) = T exp 0

(A2) where Um (t, 0) is the time-ordered propagator. Dividing the time interval t into M intervals of length t = tj − tj − 1 , Um (t, 0) can be written as a time-ordered concatenation of singleinterval propagators with time steps of length t. One can then write a single-interval propagator in terms of the following exponential operator:57

  Um ( t, 0) = exp t iL0m (t = 0) + I , (A3) where the super-operator I is a differential operator I=

F (t)e tI G(t) = F (t + t)G(t),

(A5)

where F(t) and G(t) are time-dependent operators or functions. As can be seen, I acts on objects to its left, but not to its right.57 When t is small, one can apply a symmetric Trotter factorization58 to Um (t, 0) to yield       t t 0 I exp tiLm (t = 0) exp I Um ( t, 0) = exp 2 2   = exp ( tI) exp tiL0m ( t/2) . (A6) Even though exp ( tI) no longer acts on exp ( tiL0m ( t/2)) in the second line of Eq. (A6), we retain it when building the full time-ordered concatenation of propagators for evolving the system to time t. The action of the single-interval propagator Um ( t, 0) generates trajectory segments following Hamiltonian dynamics according to Eq. (15). In order to accomplish this, we derive a low-order numerical integration scheme based on the factorization of the time-dependent propagator in Eq. (A6).57, 59 We start by decomposing the time-dependent operator in Eq. (A1) as follows: L0m (t) = L1 + L2 + L3 (t) + L4 (t),

(A7)

where P ∂ , M ∂R

(A8)

∂Hm ∂ , ∂R ∂P

(A9)

∂ 1  λλ

h (R, t) pλ

, ¯ λλ m ∂rλ

(A10)

L2 = −

In this Appendix, we derive the numerical integration scheme we used for solving the PBME, ∂ Am (X, x, t) ≈ − {Hm (t), Am (t)}X,x ≡ iL0m (t)Am (t), ∂t (A1) where iL0m (t) = − {Hm (t), }X,x is the evolution operator that generates the dynamics dictated by the Poisson bracket. The

(A4)

which acts according to

L1 = APPENDIX: NUMERICAL INTEGRATION OF EQUATIONS OF MOTION

← − ∂ , ∂t

L3 (t) =

L4 (t) = −

∂ 1  λλ

hm (R, t) rλ

. ¯ λλ

∂pλ

(A11)

144106-11

Rekik et al.

J. Chem. Phys. 138, 144106 (2013)

Using Eq. (A6), the single-interval propagator can be factorized as follows:  t I exp{ t[L1 + L2 + L3 (t) Um ( t, 0) = exp 2   t I + L4 (t)]} exp 2     t t I exp L2 exp{ t[L1 + L3 (t) = exp 2 2     t t L2 exp I + L4 (t)]} exp 2 2       t t t = exp I exp L2 exp L4 (t) 2 2 2

lution, with the mapping variables replaced by the coherent state variables.



× exp { t [L1 + L3 (t)]}   t L4 (t) × exp 2   τ ! t L2 exp I × exp 2 2   t L2 = exp ( tI) exp 2    t t L4 t + × exp 2 2  #  " t × exp t L1 + L3 t + 2      t t t L4 t + L2 . × exp exp 2 2 2 (A12) By applying the individual exponential propagators from right to left in Eq. (A12), the following integration algorithm can be written down: t Fm (0) 2 t  λλ

h (R(0), t/2) rλ (0) pλ ( t/2) = pλ (0) − 2¯ λ m P ( t/2) = P (0) +

P ( t/2) (A13) M τ  λλ

rλ ( t) = rλ (0) + h (R( t), t/2) pλ ( t/2) ¯ λ m

R ( t) = R (0) + t

pλ ( t) = pλ ( t/2) −

t  λλ

h (R( t), t/2) rλ ( t) 2¯ λ m

P ( t) = P ( t/2) +

t Fm ( t). 2

Using an integration step t = 2.5 × 10−5 in scaled units, this algorithm works well for the calculation of the trajectories needed for the population simulations. It should be noted that an analogous integration algorithm is used for the FB so-

1 J.

C. Tully, J. Chem. Phys. 93, 1061 (1990). C. Tully, in Modern Methods for Multidimensional Dynamics Computations in Chemistry, edited by D. L. Thompson (World Scientific, New York, 1998), p. 34. 3 F. Webster, E. T. Wang, P. J. Rossky, and R. A. Friesner, J. Chem. Phys. 100, 4835 (1994). 4 O. V. Prezhdo, J. Chem. Phys. 111, 8366 (1999). 5 A. W. Jasper, C. Zhu, S. Nangia, and D. G. Truhlar, Faraday Discuss. 127, 1 (2004). 6 M. J. Bedard-Hearn, R. E. Larsen, and B. J. Schwartz, J. Chem. Phys. 123, 234106 (2005). 7 S. A. Fischer, C. T. Chapman, and X. Li, J. Chem. Phys. 135, 144102 (2011). 8 M. P. Miller and C. W. McCurdy, J. Chem. Phys. 69, 5163 (1978). 9 H. D. Meyer and W. H. Miller, J. Chem. Phys. 70, 3214 (1979). 10 X. Sun, H. Wang, and W. H. Miller, J. Chem. Phys. 109, 7064 (1998). 11 G. Stock and M. Thoss, Phys. Rev. Lett. 78, 578 (1997). 12 U. Müller and G. Stock, J. Chem. Phys. 108, 7516 (1998). 13 M. Thoss and G. Stock, Phys. Rev. A 59, 64 (1999). 14 K. Thompson and N. Makri, Phys. Rev. E 59, R4729 (1999). 15 W. H. Miller, J. Phys. Chem. A 105, 2942 (2001). 16 M. Thoss and H. Wang, Annu. Rev. Phys. Chem. 55, 299 (2004). 17 G. Stock and M. Thoss, Adv. Chem. Phys. 131, 243 (2005). 18 S. Bonella and D. F. Coker, J. Chem. Phys. 118, 4370 (2003). 19 S. Bonella and D. F. Coker, J. Chem. Phys. 122, 194102 (2005). 20 S. Bonella, D. Montemayor, and D. F. Coker, Proc. Natl. Acad. Sci. U.S.A. 102, 6715 (2005). 21 E. Dunkel, S. Bonella, and D. F. Coker, J. Chem. Phys. 129, 114106 (2008). 22 P. Huo and D. F. Coker, J. Chem. Phys. 135, 201101 (2011). 23 C. C. Martens and J. Fang, J. Chem. Phys. 106, 4918 (1997). 24 A. Donoso and C. C. Martens, J. Chem. Phys. 102, 4291 (1998). 25 A. Donoso and C. C. Martens, J. Chem. Phys. 112, 3980 (2000). 26 A. Donoso, D. Kohen, and C. C. Martens, J. Chem. Phys. 112, 7345 (2000). 27 M. Santer, U. Manthe, and G. Stock, J. Chem. Phys. 114, 2001 (2001). 28 R. Kapral and G. Ciccotti, J. Chem. Phys. 110, 8919 (1999). 29 S. Nielsen, R. Kapral, and G. Ciccotti, J. Chem. Phys. 112, 6543 (2000). 30 C. Wan and J. Schofield, J. Chem. Phys. 113, 7047 (2000). 31 C. Wan and J. Schofield, J. Chem. Phys. 112, 4447 (2000). 32 R. Kapral, Annu. Rev. Phys. Chem. 57, 129 (2006). 33 H. Kim, A. Nassimi, and R. Kapral, J. Chem. Phys. 129, 084102 (2008). 34 A. Kelly, R. van Zon, J. M. Schofield, and R. Kapral, J. Chem. Phys. 136, 084101 (2012). 35 C.-Y. Hsieh and R. Kapral, J. Chem. Phys. 137, 22A507 (2012). 36 J. C. Tully, in Classical and Quantum Dynamics in Condensed Phase Simulations, edited by B. J. Berne, G. Ciccotti, and D. F. Coker (World Scientific, Singapore, 1998). 37 M. F. Herman, Annu. Rev. Phys. Chem. 45, 83 (1994). 38 D. MacKernan, R. Kapral, and G. Ciccotti, J. Phys.: Condens. Matter 14, 9069 (2002). 39 D. MacKernan, G. Ciccotti, and R. Kapral, J. Phys. Chem. B 112, 424 (2008). 40 I. Horenko, C. Salzmann, B. Schmidt, and C. Schutte, J. Chem. Phys. 117, 11075 (2002). 41 C. Wan and J. Schofield, J. Chem. Phys. 116, 494 (2002). 42 G. Hanna and R. Kapral, J. Chem. Phys. 122, 244505 (2005). 43 G. Stock, J. Chem. Phys. 103, 1561 (1995). 44 A. Nassimi, S. Bonella, and R. Kapral, J. Chem. Phys. 133, 134115 (2010). 45 A. Kelly and Y. M. Rhee, J. Phys. Chem. Lett. 2, 808 (2011). 46 R. Glauber, Phys. Rev. 131, 2766 (1963). 47 H. Wang, M. Thoss, and W. H. Miller, J. Chem. Phys. 115, 2979 (2001). 48 H. B. Wang and M. Thoss, J. Phys. Chem. A 107, 2126 (2003). 49 H. B. Wang and M. Thoss, Chem. Phys. Lett. 389, 43 (2004). 50 E. Wigner, Phys. Rev. 40, 749 (1932). 51 J. Schwinger, in Quantum Theory of Angular Momentum, edited by L. C. Biedenharn, and H. V. Dam (Academic, New York, 1965), p. 229. 52 H. Wang, X. Song, D. Chandler, and W. H. Miller, J. Chem. Phys. 110, 4828 (1999). 53 P. Huo and D. F. Coker, Mol. Phys. 110, 1035 (2012). 2 J.

144106-12

Rekik et al.

J. Chem. Phys. 138, 144106 (2013)

54 H.

59 A.

55 P.

60 The

Kim and R. Kapral, Chem. Phys. Lett. 423, 76 (2006). Huo and D. F. Coker, J. Chem. Phys. 137, 22A535 (2012). 56 C.-Y. Hsieh and R. Kapral, J. Chem. Phys. 138, 134110 (2013). 57 M. Suzuki, Proc. Jpn. Acad., Ser. B: Phys. Biol. Sci. 69, 161 (1993). 58 H. F. Trotter, Proc. Am. Math. Soc. 10, 545 (1959).

Ricci and G. Ciccotti, Mol. Phys. 101, 1927 (2003). comparison to the transient absorption data is made since the population data are not available in Ref. 49. Since these two observables are different, the exact results are not expected to be identical, so they should only be qualitatively compared.

The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp