Calculation of absorption spectra involving multiple excited states: Approximate methods based on the mixed quantum classical Liouville equation Shuming Bai, Weiwei Xie, Lili Zhu, and Qiang Shi Citation: The Journal of Chemical Physics 140, 084105 (2014); doi: 10.1063/1.4866367 View online: http://dx.doi.org/10.1063/1.4866367 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/8?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Analysis of the forward-backward trajectory solution for the mixed quantum-classical Liouville equation J. Chem. Phys. 138, 134110 (2013); 10.1063/1.4798221 Analysis of the quantum-classical Liouville equation in the mapping basis J. Chem. Phys. 133, 134115 (2010); 10.1063/1.3480018 Dynamics of coupled Bohmian and phase-space variables: A moment approach to mixed quantum-classical dynamics J. Chem. Phys. 122, 094103 (2005); 10.1063/1.1856462 A derivation of the mixed quantum-classical Liouville equation from the influence functional formalism J. Chem. Phys. 121, 3393 (2004); 10.1063/1.1771641 Fully adaptive propagation of the quantum-classical Liouville equation J. Chem. Phys. 120, 8913 (2004); 10.1063/1.1691015
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
THE JOURNAL OF CHEMICAL PHYSICS 140, 084105 (2014)
Calculation of absorption spectra involving multiple excited states: Approximate methods based on the mixed quantum classical Liouville equation Shuming Bai, Weiwei Xie, Lili Zhu, and Qiang Shia) Beijing National Laboratory for Molecular Sciences, State Key Laboratory for Structural Chemistry of Unstable and Stable Species, Institute of Chemistry, Chinese Academy of Sciences, Zhongguancun, Beijing 100190, China
(Received 16 November 2013; accepted 10 February 2014; published online 26 February 2014) We investigate the calculation of absorption spectra based on the mixed quantum classical Liouville equation (MQCL) methods. It has been shown previously that, for a single excited state, the averaged classical dynamics approach to calculate the linear and nonlinear spectroscopy can be derived using the MQCL formalism. This work focuses on problems involving multiple coupled excited state surfaces, such as in molecular aggregates and in the cases of coupled electronic states. A new equation of motion to calculate the dipole-dipole correlation functions within the MQCL formalism is first presented. Two approximate methods are then proposed to solve the resulted equations of motion. The first approximation results in a mean field approach, where the nuclear dynamics is governed by averaged forces depending on the instantaneous electronic states. A modification to the mean field approach based on first order moment expansion is also proposed. Numerical examples including calculation of the absorption spectra of Frenkel exciton models of molecular aggregates, and the pyrazine molecule are presented. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4866367] I. INTRODUCTION
Mixed quantum-classical treatments of the molecular dynamics are widely employed in theoretical studies of problems involving electronic transitions.1–3 In multi-dimensional systems, the mixed quantum-classical schemes are computationally more efficient than full quantum treatments such as the multi-configuration time dependent Hartree (MCTDH) method,4, 5 yet are reasonably accurate to describe the nonadiabatic dynamics. Another advantage of the mixed quantum-classical simulations is that, they can be effectively combined with the on-the-fly quantum chemistry calculations of forces and electronic couplings, without the need to fit multi-dimensional potential energy surfaces. An important application of mixed quantum classical simulations is in calculations of the linear and nonlinear electronic and vibrational spectroscopic signals, which are important in revealing the structure and dynamics of complex molecules and condensed phase systems.6, 7 Here, we restrict our discussion of mixed quantum classical methods to those treating the electronic degrees of freedom (or the high frequency vibrational modes of interest) quantum mechanically, while the (other) nuclear degrees of freedom classically. This kind of studies of the absorption spectra can date back to the Kubo-Anderson theory of line shapes,8, 9 where the quantum system is subject to stochastic fluctuations caused by a classical environment. Later, semiclassical schemes are adopted to study various types of linear and nonlinear spectroscopy, based on calculating the optical response functions.6 In current studies, a typical approach is based on the time a) Electronic mail: [email protected]
0021-9606/2014/140(8)/084105/10/$30.00
dependent nuclear dynamics obtained from equilibrium classical molecular dynamics simulations performed on the ground state potential surface. Following the same nomenclature as in Ref. 10, we will refer to this approach as the dynamical classical (DCL) method. The DCL method has now been widely applied in studying coupled quantum subsystems, such as vibrational spectra of liquids11–14 and biological molecules,15, 16 and electronic spectra of molecular aggregates.17–19 A major problem of the DCL method is that, the energy fluctuations of the quantum system are described only by classical dynamics on the ground state, i.e., the excited state energy surface only enters the calculation of the instantaneous transition energies, and does not have any effects on the nuclear dynamics. Drawback of such treatment has been analyzed previously for spectra of a single chromophore,20, 21 and recently been tested in the case of absorption spectra of model dimers.22 Knoester, Jansen, and co-workers have also shown that the DCL method gives incorrect dynamics in simulating the two-dimensional (2D) spectra of molecular aggregates, and proposed new methods based on the Ehrenfest method23 and surface hopping method24 to calculate the third order optical response functions. A method that improves over the DCL approach is to run simulation on the average potential energy surface of the ground and excited states, which is often termed as the ACL (averaged classical) approach.10, 20 (If combined with the Wigner distribution of the classical degrees of freedom for initial sampling, it is called the Wigner averaged classical [WACL] method.20, 25 ) The WACL method happens to be exact for displaced harmonic-oscillator models, and was shown to work well in general anharmonic systems.20, 21 But
140, 084105-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
084105-2
Bai et al.
J. Chem. Phys. 140, 084105 (2014)
this method of propagation on the average potential energy surface is only derived for two surface (i.e., the ground and one excited state) problems. When transitions to multiple coupled electronic states are involved, such as in the case of conical intersections,26 or in the case of coupled chromophores in molecular aggregates,27 there is no unique average surface for calculating the optical response functions. Therefore, further development of mixed quantum classical methods to calculate the optical response functions involving multiple excited states is still needed. In this work, we focus on the problem of absorption spectra, where the dipole-dipole correlation function is needed to calculate the spectra. We will show that, by using the general mixed quantum classical Liouville equation (MQCL) formalism, a set of equations can be obtained to calculate the dipoledipole correlation functions in the multiple surface case. This set of equations reduce to the WACL approach in the single excited state case. However, in the case of multiple coupled excited states, solving this set of multidimensional partial differential equations is not a simple task as in the case of a single excited state. We then propose two approximate methods based on trajectory descriptions of MQCL equation to calculate the dipole-dipole correlation functions. The first one is a mean field method where the average forces are calculated based on the instantaneous electronic states. A modified mean field approach is then obtained using the mean field trajectory as a reference. The new methods are tested on calculations of Frenkel exciton models of molecular aggregates, as well as 4-mode and 24-mode models of pyrazine.26 The remaining sections of this paper are arranged as follows. In Secs. II A and II B, we use the MQCL equations to derive the expression of mixed quantum-classical method for the dipole-dipole correlation function in calculating absorption spectra. In Sec. II C, we propose the mean field method based on the MQCL equation. A modified mean field approach is then derived in Sec. II D. In Sec. III, the absorption spectra of model dimers, the Fenna-Matthews-Olson (FMO) complex, and the pyrazine molecule are calculated using different methods and the results are compared. The conclusions and discussions are made in Sec. IV. II. THEORY A. The mixed quantum-classical Liouville equation
Consider a general Hamiltonian describing a quantum subsystem coupled to a group of heavy particles that are to be treated classically, 2 pˆ 2 Pˆj ˆ ˆ Q), Hˆ = + Vˆ (q, + 2μ 2Mj j
the quantum Liouville equation.29, 31, 33 First, the quantum Liouville equation for the density operator of the total system ρ T is given by i ∂ρT = − [Hˆ , ρT ]. ∂t ¯
(2)
The partial Wigner transform of ρ T is defined as 1 i T ρw (Q, P) = de− ¯ P · Q + /2|ρˆT |Q − /2, N (2π ¯) (3) where the integration is performed for the degrees of freedom that are to be treated classically. We note that ρw defined in Eq. (3) is still a density operator of the quantum subsystem. A direct derivation of the MQCL utilizes the following formally exact equation for the partial Wigner transform of ˆ 6, 34 the products of two operators Aˆ and B: ˆ w = Aw exp i¯ Bw = Bw exp − i¯ Aw , (Aˆ B) 2 2 (4) where Aw and Bw are the partial Wigner transforms of the opˆ is the Poisson bracket. For a pair of phase erators Aˆ and B, space distributions Aw (Q, P) and Bw (Q, P), is defined as ∂Aw ∂Bw ∂Aw ∂Bw . Aw Bw = {Aw , Bw } ≡ − ∂Qj ∂Pj ∂Pj ∂Qj j (5) By doing partial Wigner transform of the quantum Liouville equations in Eq. (2), and applying the classical limit ¯ → 0 to Eq. (4), the MQCL is given by29, 31, 32 i ∂ρw 1 = − [Hˆ w , ρw ] + [{Hˆ w , ρw } − {ρw , Hˆ w }], ∂t ¯ 2
(6)
where Hˆ w is the partial Wigner transform of the total Hamiltonian Hˆ . The above Eq. (6) is still in the operator form for the quantum subsystem. To perform the actual simulations, one will need to choose a basis set for the quantum degrees of freedom. Very often, equations of motion in the adiabatic basis are used.29, 31, 33 For the problems considered in this work, i.e., the Frenkel exciton model for molecular aggregates, and the case with a conical intersection, the calculations are performed in the diabatic basis. We note that extension of the current work to the adiabatic basis can be rather straightforward.
(1)
ˆ qˆ are the momentum and coordinate of the quantum where p, ˆ Q ˆ are those of the heavy particles, and μ and subsystem, P, ˆ is the coupling beˆ Q) Mj are the corresponding masses. Vˆ (q, tween the system and bath. For simplicity, only one degree of freedom for the quantum subsystem is considered. The MQCL equations can be derived through many different ways.28–32 For the reason of completeness, we briefly sketch the derivation using the partial Wigner transform of
B. Absorption spectra from the MQCL equation
In this subsection, we show how the MQCL equations can be applied to calculation of absorption spectra in a problem with several excited electronic states. To this end, we consider a model where the electronic degree of freedom is treated quantum mechanically, and the nuclear degrees of freedom are treated classically. The total molecular Hamiltonian plus the interaction with a classical
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
084105-3
Bai et al.
J. Chem. Phys. 140, 084105 (2014)
dynamics of the N coherence terms |mg| are decoupled from other terms of the density matrix including |gg|, |gm|, and |mn|. We further define ρwm (P, Q; t) = m|(μρg )w (P, Q; t)|g, the correlation function in Eq. (10) can then be calculated as N 1 C(t) = μ(t)μ(0)g = μm ρwm (Q, P; t). dPdQ (2π ¯)N m=1
electromagnetic (EM) field can be written as Hˆ = Hˆ g (P, Q)|gg| +
N
Hˆ m (P, Q)|mm|
m=1
+
N
Hˆ mn (P, Q)|mn|
m=n=1
− E(t) ·
N
μm (Q)[|mg| + |gm|].
m=1
In the above equation, |g and |m stand for the ground and the mth excited states, Q and P are the Cartesian coordinates and momenta of the molecule system, E(t) is the external EM field, and μm (Q) = μgm (Q) is the transition dipole operator, Hmn is the coupling between the mth and nth electronic states. The total Hamiltonian in Eq. (7) has neglected the coupling between ground and excited states, which is in general a good approximation and used in the previous works,10, 20 since the energy difference between ground and excited states is usually large.26, 35 When the EM field is weak, the absorption spectra can be calculated using the dipole-dipole correlation function defined as6 +∞ 1 I (ω) ∝ dteiωt μ(t) · μ(0)g , (8) 2π −∞ where μ = N m=1 μm (|gm| + |mg|) is the total dipole operator, the subscript g in the average indicates that the initial state to calculate the correlation function is equilibrated on the ground electronic state, ρ g = e−βHg /Tre−βHg (β = 1/kB T is the inverse temperature), and Hmol is the molecular Hamiltonian without the EM field, Hˆ mol = Hˆ g (P, Q)|gg| +
N
Hˆ mn (P, Q)|mn|,
(11)
(7)
(9)
Now we derive the equation of motion for ρwm . Using the general MQCL equation in Eq. (6), and since the coupling terms between ground state and excited states are zero, the dynamics of ρwm (Q, P) is governed by d m i ρ (Q, P; t) = − [Hw , ρw (Q, P; t)]mg dt w ¯ 1 ∂Hmn ∂ρwn ∂ρwn ∂Hmn + · − · 2 n ∂Q ∂P ∂Q ∂P 1 ∂Hg ∂ρwn ∂ρ n ∂Hg + · − w · . (12) 2 ∂Q ∂P ∂Q ∂P In the diabatic basis, the off-diagonal Hmn terms do P2 not depend on momentum. We also assume that Hg = 2M P2 + Vg (Q), Hm = 2M + Vm (Q). The N ρwm (Q, P; t) terms now form a closed equation, d m i ρ (Q, P; t) = − [Vm (Q) − Vg (Q)]ρwm (Q, P) dt w ¯ i − Hmn ρwn (Q, P) ¯ n=m 1 ∂Hmn ∂ρwn + {Hm , ρwm (P, Q)}+ · , 2 n=m ∂Q ∂P (13)
m,n=1
where we have used the notations Hˆ mm = Hˆ m when m = n. We then investigate how the absorption spectra can be calculated using the MQCL formalism. We work with the Condon approximation by assuming that the transition dipole moment does not depend on the nuclear degrees of freedom, while extensions to the non-Condon cases are straightforward. For simplicity, we will also drop the orientation dependence of μ in further derivations. The first step is to write the correlation function in Eq. (8) in the form of partial Wigner transformation, 1 dPdQμw (μρg )w (Q, P; t), μ(t)μ(0)g = TrS (2π ¯)N (10) N where μw = μ = m=1 μm (|gm|+|mg|), (μρg )w (Q, P; t) is the partial Wigner transform of the operator ˆ ˆ e−i H t/¯ (μˆ ρˆg )ei H t/¯ , which can be calculated using the MQCL approximation in Eq. (4). The initial state to the partial Wigner transform calculate N −βHg μ e /Tre−βHg |mg|. Since (μρg )w (P, Q; t) is m=1 m there is no direct coupling between the ground and excited states in the molecular Hamiltonian in Eq. (9), the
where the effective Hamiltonian Hm for the mth surface P2 is defined as Hm = 2M + V˜m (Q), with V˜m (Q) = 12 [Vg (Q) + Vm (Q)], and {Hm , ρwm (P, Q)} = −
∂ V˜m ∂ρwm P ∂ρwm · + · . M ∂Q ∂Q ∂P
(14)
The four terms on the right-hand side of Eq. (13) have different effects on the evolution of ρwm : (1) The first term gives a phase factor related to energy difference Vm (Q) − Vg (Q); (2) In the second term, the different ρwm components are coupled by the off-diagonal terms Hmn of the Hamiltonian; (3) The Poisson bracket in the third term indicates that the distribution ρwm (Q, P) evolves on the potential energy surface defined by V˜m , which is the average of the ground state potential energy surface Vg and the mth excited state Vm ; (4) The fourth term represents the “off-diagonal forces” originated from the derivative of the Hmn (Q) terms. Equation (13) consists of a set of partial differential equations for the Wigner distribution functions, which is generally hard to solve in the multi-dimensional case. However, if there is only one excited state, the correlation function can be calculated rather straightforwardly. In such cases, there is only
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
084105-4
Bai et al.
J. Chem. Phys. 140, 084105 (2014)
one term ρ = ρ eg , and the correlation function in Eq. (8) is 1 C(t) ≈ (15) dQdPμρw (Q, P; t); (2π ¯)N and the equation of motion for ρw (P, Q; t) is given by d i ρw (Q, P; t) = − [Ve (Q) − Vg (Q)]ρw (Q, P; t) dt ¯ + {Hav , ρw (Q, P; t)}, P2 2M
(16)
Ve (Q)+Vg (Q) 2
+ is the classical Hamiltonian where Hav = on the averaged potential energy surface. It is easy to find that t 1 i C(t) ≈ dQdPμρw (Q, P; 0)e− ¯ 0 dτ U (Q(τ )) , N (2π ¯) (17) where the dynamics is calculated on the average potential energy surface of the ground and excited states, and U(Q) is the difference between the excited state and ground state potential. This is just the WACL approach obtained in previous studies.20, 25 By neglecting the off-diagonal force terms in Eq. (13), and assuming that ρwm (P, Q; t) evolves on the ground state potential energy surface rather than V˜m , it can be shown that the correlation function in Eq. (8) can be calculated by the commonly used DCL method for multiple excited states11–14 ∞ I (ω) ∝ Re dte−iωt μm (0)Fmn (t)μn (t), (18) 0
mn
where the terms Fmn are calculated as dF (t) = iF (t) (t), (19) dt with the initial condition that Fmn (0) = δ mn , and the matrix is give by ¯ mn (t) = (Ei − Eg )(t)δmn + Hmn (t)(1 − δmn ).
(20)
The average is taken over the Boltzmann distribution of nuclear degrees of freedom on the ground state, and μ(t), (t) are now classical variables, whose dynamics are governed by the ground state nuclear Hamiltonian. C. The mean field approximation
As stated above, the goal is to derive a trajectory based method to solve Eqs. (11) and (13). Although this can be done easily for the single excited state problem as shown in Eq. (17), the same method can not be extended directly to the case of multiple coupled excited surfaces. In literature, there are several numerical algorithms developed to solve the MQCL equations based on surface-hopping type trajectory simulations.31, 33, 36–38 However, these methods usually suffer from large numerical noise at longer simulation times, which will cause noise problems in calculating the Fourier transform in Eq. (8). In this section, we derive a mean field approximation to solve the MQCL equation in Eq. (13). According to the previous work,28, 39, 40 the mean field approximation for the general MQCL equation in Eq. (6) is briefly presented in the Appendix. As shown below, the derivation of the mean
field approximation for the dipole-dipole correlation functions in Eq. (13) uses a similar idea, but with several technical differences. In the Appendix, the expectation values of classical variables Q and P are defined using ρ mm , which is the population on each state [see Eqs. (A2) and (A3)]. However, since in Eq. (13), only off-diagonal terms ρwm (P, Q; t) = m|(μρg )w (P, Q; t)|g are involved, a new definition of the mean field positions and momenta is needed. Fortunately, we can derive the following normalizing equation about the ρwm terms from Eq. (13): 1 ρwm ρwm∗ = C. (21) dPdQ N (2π ¯) m To develop a method with trajectory description, the following integral of ρwm (Q, P; t) is defined: ρ¯m (t) = dQdPρwm (Q, P; t), (22) and by using Eq. (21), we define the averages of the positions Q and momenta P as 1 ¯ = ρwm (Q, P)ρwm∗ (Q, P)Q, (23) Q dQdP (2π ¯)N C m P¯ =
1 (2π ¯)N C
dQdP
ρwm (Q, P)ρwm∗ (Q, P)P,
(24)
m
where the constant C is defined in Eq. (21). So the difference between the current derivation and the conventional mean field approximation in the Appendix is that, instead of using the averages of Q and P defined using the diagonal elements of the density matrix, Eqs. (23) and (24) have used a “wave function” like definition of the averages. The rest of the derivation is then very similar to the conventional mean field approximation. We use the following ansatz for ρ j (Q, P): ¯ ¯ P − P(t)). ρwm (Q, P) = ρwm (t)δ(Q − Q(t),
(25)
Together with Eq. (13) and the definitions in Eqs. (22)– (24), the equations of the motion which constitute the mean field approximation to calculate the dipole-dipole correlation functions can be obtained as follows: i i d ¯ ¯ ¯ ρ¯n , (26) ρ¯m = − [Vm (Q)−V Hmn (Q) g (Q)]ρ¯m − dt ¯ ¯ n=m and d ¯ P¯ , Q= dt M −1 ¯ d ¯ 1 ∂Vg (Q) 1 2 |ρ¯m | − P = Fmf = − ¯ dt 2 ∂Q 2 m
¯ ∂Hmn (Q) . × Re ρ¯m∗ ρ¯n ¯ ∂Q n
(27)
(28)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
084105-5
Bai et al.
J. Chem. Phys. 140, 084105 (2014)
D. The modified mean field method
The shortcomings of the mean field dynamics have been widely investigated.3 Now we work on possible extensions of the mean field approach. We note that a major problem of the mean field ansatz is that, in the multiple surface case, trajectories associated with different ρwm (Q, P) should propagate on different potential surfaces (see Eq. (14)). The mean field ansatz then becomes problematic, since even though the initial positions and momenta are the same on different surfaces, they soon become different evolving under the influence of different forces. Nevertheless, we can still take the mean field dynamics as a reference, and investigate further modifications. To derive the modified mean field method, we assume that each ρwm term is slightly different from the mean field ansatz in Eq. (25). So we can still use the mean field trajec¯ ¯ as the zeroth order approximatory defined by Q(t) and P(t) tion, and obtain the correction terms though expansions of the phase space variables P, Q, and Hmn (P, Q), near the reference ¯ ¯ trajectory Q(t), P(t). The following terms are defined as deviations of the Q and P expectation values from the mean field ansatz: ¯ m (Q, P), Xm = dQdP(Q − Q)ρ (29) w m =
¯ m (Q, P). dQdP(P − P)ρ w
(30)
In this study, we only include the first order moments for simplicity, and high order moments are not considered. Several other expansions of the phase space average are needed in deriving the modified mean field equation. Calculation of the dQdPHmn (Q)ρwn (Q, P) term can be approximated as dQdPHmn (Q)ρwn (Q, P) ≈
¯ n (Q, P) dQdPHmn (Q)ρ w
¯ ∂Hmn (Q) ¯ n (Q, P) (Q − Q)ρ w ¯ ∂Q ¯ ¯ ρ¯n + ∂Hmn (Q) · Xm . = Hmn (Q) (31) ¯ ∂Q Terms like dQdPHmn (Q)ρwn (Q, P)Q can be approximated in the following way: dQdPHmn (Q)ρwn (Q, P)Q +
dQdP
≈
¯ Qρ ¯ n (Q, P) dQdPHmn (Q) w +
¯ ¯ n (Q, P) ¯ Q ¯ ∂Hmn (Q) (Q− Q)ρ dQdP Hmn (Q)+ w ¯ ∂Q
¯ Q ¯ ρ¯n + Hmn (Q)X ¯ m+Q ¯ = Hmn (Q)
¯ ∂Hmn (Q) · Xm . ¯ ∂Q
(32)
Using the approximations similar to the ones in the above two equations, further derivation of the modified mean field
equations can be obtained by taking proper phase space averages over Eq. (13). For example, the dynamics for ρ¯m is obtained by combining the integration over Eqs. (13) and (31): i d ¯ − Vg (Q)] ¯ ρ¯m − i ¯ ρ¯n ρ¯m = − [Vm (Q) Hmn (Q) dt ¯ ¯ n=m
¯ − Vg (Q)] ¯ i ∂[Vm (Q) − ¯ ¯ ∂Q
¯ i ∂Hmn (Q) · Xn , · Xm − ¯ ¯ n=m ∂Q
(33)
which adds correction terms involving X to the original mean field result in Eq. (26). Equations (27) and (28) define the reference mean field trajectory and stay the same. The following equation of motion for Xm can be derived by multiplying Q to Eq. (13) and doing the integration, and utilizing Eq. (32): i d m i ¯ − Vg (Q)]X ¯ ¯ n. Xm = − [Vm (Q) Hmn (Q)X m− dt M ¯ ¯ n=m (34) The equation of motion for involves second order derivatives of Hmn (Q), which is given by d i i ¯ − Vg (Q)] ¯ m = − [Vm (Q) Hmn n − ρ¯m Fmf m− dt ¯ ¯ n=m −
¯ + Vg (Q)] ¯ ¯ 1 ∂[Vm (Q) 1 ∂Hmn (Q) ρ¯m − ρ¯n ¯ ¯ 2 2 n=m ∂ Q ∂Q
−
¯ ¯ ¯ 1 ∂ 2 [Vm (Q)+V 1 ∂ 2 Hmn (Q) g (Q)] ·X − ·Xn , m 2 2 ¯ ¯ 2 2 n=m ∂ Q ∂Q (35)
where Fmf is the expression of force in mean field approximation, as in Eq. (28); and order derivative terms the second mn (Q) m ρw (Q, P). come from integral like dQdP ∂H∂Q We note that this modified mean field method is similar to the modified Ehrenfest dynamics obtained by Horsfield and co-workers,41, 42 which is also an equation of motion for first order moments, and has been used for calculations of coupled electron-ion dynamics in inelastic quantum transport. III. RESULTS
We now apply the above methods to calculate the dipoledipole correlation functions and absorption spectra of several model systems, which include (1) a model of molecular aggregate dimers; (2) a model of the FMO complex; and (3) models of the pyrazine molecule with S1 -S2 conical intersection, including the 4-mode and 24-mode models. The results calculated using the mean field and modified mean field methods are compared with benchmark calculations using the Hierarchical equations of motion (HEOM) method,43, 44 and the MCTDH method.26, 45, 46 In the simulations, Eqs. (26), (27), and (28) are solved for the mean field method. The initial positions and momenta of the trajectory are sampled from the Wigner distribution on
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
Bai et al.
J. Chem. Phys. 140, 084105 (2014)
A. Absorption spectra of molecular aggregates
We consider in this subsection the Frenkel exciton model of N two-level molecules coupled to a phonon bath, where the total Hamiltonian is written as H = He + Hph + He−ph ,
(36)
where He describes the electronic degrees of freedom, Hph describes the nuclear degrees of freedom, and He−ph describes their couplings. To calculate the absorption spectra, only the ground state and the single excited states are considered: He = 0 |00|+
N
m |mm|+
N N
Jmn (|mn|+|nm|),
m=1 n
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
THE JOURNAL OF CHEMICAL PHYSICS 140, 084105 (2014)
Calculation of absorption spectra involving multiple excited states: Approximate methods based on the mixed quantum classical Liouville equation Shuming Bai, Weiwei Xie, Lili Zhu, and Qiang Shia) Beijing National Laboratory for Molecular Sciences, State Key Laboratory for Structural Chemistry of Unstable and Stable Species, Institute of Chemistry, Chinese Academy of Sciences, Zhongguancun, Beijing 100190, China
(Received 16 November 2013; accepted 10 February 2014; published online 26 February 2014) We investigate the calculation of absorption spectra based on the mixed quantum classical Liouville equation (MQCL) methods. It has been shown previously that, for a single excited state, the averaged classical dynamics approach to calculate the linear and nonlinear spectroscopy can be derived using the MQCL formalism. This work focuses on problems involving multiple coupled excited state surfaces, such as in molecular aggregates and in the cases of coupled electronic states. A new equation of motion to calculate the dipole-dipole correlation functions within the MQCL formalism is first presented. Two approximate methods are then proposed to solve the resulted equations of motion. The first approximation results in a mean field approach, where the nuclear dynamics is governed by averaged forces depending on the instantaneous electronic states. A modification to the mean field approach based on first order moment expansion is also proposed. Numerical examples including calculation of the absorption spectra of Frenkel exciton models of molecular aggregates, and the pyrazine molecule are presented. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4866367] I. INTRODUCTION
Mixed quantum-classical treatments of the molecular dynamics are widely employed in theoretical studies of problems involving electronic transitions.1–3 In multi-dimensional systems, the mixed quantum-classical schemes are computationally more efficient than full quantum treatments such as the multi-configuration time dependent Hartree (MCTDH) method,4, 5 yet are reasonably accurate to describe the nonadiabatic dynamics. Another advantage of the mixed quantum-classical simulations is that, they can be effectively combined with the on-the-fly quantum chemistry calculations of forces and electronic couplings, without the need to fit multi-dimensional potential energy surfaces. An important application of mixed quantum classical simulations is in calculations of the linear and nonlinear electronic and vibrational spectroscopic signals, which are important in revealing the structure and dynamics of complex molecules and condensed phase systems.6, 7 Here, we restrict our discussion of mixed quantum classical methods to those treating the electronic degrees of freedom (or the high frequency vibrational modes of interest) quantum mechanically, while the (other) nuclear degrees of freedom classically. This kind of studies of the absorption spectra can date back to the Kubo-Anderson theory of line shapes,8, 9 where the quantum system is subject to stochastic fluctuations caused by a classical environment. Later, semiclassical schemes are adopted to study various types of linear and nonlinear spectroscopy, based on calculating the optical response functions.6 In current studies, a typical approach is based on the time a) Electronic mail: [email protected]
0021-9606/2014/140(8)/084105/10/$30.00
dependent nuclear dynamics obtained from equilibrium classical molecular dynamics simulations performed on the ground state potential surface. Following the same nomenclature as in Ref. 10, we will refer to this approach as the dynamical classical (DCL) method. The DCL method has now been widely applied in studying coupled quantum subsystems, such as vibrational spectra of liquids11–14 and biological molecules,15, 16 and electronic spectra of molecular aggregates.17–19 A major problem of the DCL method is that, the energy fluctuations of the quantum system are described only by classical dynamics on the ground state, i.e., the excited state energy surface only enters the calculation of the instantaneous transition energies, and does not have any effects on the nuclear dynamics. Drawback of such treatment has been analyzed previously for spectra of a single chromophore,20, 21 and recently been tested in the case of absorption spectra of model dimers.22 Knoester, Jansen, and co-workers have also shown that the DCL method gives incorrect dynamics in simulating the two-dimensional (2D) spectra of molecular aggregates, and proposed new methods based on the Ehrenfest method23 and surface hopping method24 to calculate the third order optical response functions. A method that improves over the DCL approach is to run simulation on the average potential energy surface of the ground and excited states, which is often termed as the ACL (averaged classical) approach.10, 20 (If combined with the Wigner distribution of the classical degrees of freedom for initial sampling, it is called the Wigner averaged classical [WACL] method.20, 25 ) The WACL method happens to be exact for displaced harmonic-oscillator models, and was shown to work well in general anharmonic systems.20, 21 But
140, 084105-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
084105-2
Bai et al.
J. Chem. Phys. 140, 084105 (2014)
this method of propagation on the average potential energy surface is only derived for two surface (i.e., the ground and one excited state) problems. When transitions to multiple coupled electronic states are involved, such as in the case of conical intersections,26 or in the case of coupled chromophores in molecular aggregates,27 there is no unique average surface for calculating the optical response functions. Therefore, further development of mixed quantum classical methods to calculate the optical response functions involving multiple excited states is still needed. In this work, we focus on the problem of absorption spectra, where the dipole-dipole correlation function is needed to calculate the spectra. We will show that, by using the general mixed quantum classical Liouville equation (MQCL) formalism, a set of equations can be obtained to calculate the dipoledipole correlation functions in the multiple surface case. This set of equations reduce to the WACL approach in the single excited state case. However, in the case of multiple coupled excited states, solving this set of multidimensional partial differential equations is not a simple task as in the case of a single excited state. We then propose two approximate methods based on trajectory descriptions of MQCL equation to calculate the dipole-dipole correlation functions. The first one is a mean field method where the average forces are calculated based on the instantaneous electronic states. A modified mean field approach is then obtained using the mean field trajectory as a reference. The new methods are tested on calculations of Frenkel exciton models of molecular aggregates, as well as 4-mode and 24-mode models of pyrazine.26 The remaining sections of this paper are arranged as follows. In Secs. II A and II B, we use the MQCL equations to derive the expression of mixed quantum-classical method for the dipole-dipole correlation function in calculating absorption spectra. In Sec. II C, we propose the mean field method based on the MQCL equation. A modified mean field approach is then derived in Sec. II D. In Sec. III, the absorption spectra of model dimers, the Fenna-Matthews-Olson (FMO) complex, and the pyrazine molecule are calculated using different methods and the results are compared. The conclusions and discussions are made in Sec. IV. II. THEORY A. The mixed quantum-classical Liouville equation
Consider a general Hamiltonian describing a quantum subsystem coupled to a group of heavy particles that are to be treated classically, 2 pˆ 2 Pˆj ˆ ˆ Q), Hˆ = + Vˆ (q, + 2μ 2Mj j
the quantum Liouville equation.29, 31, 33 First, the quantum Liouville equation for the density operator of the total system ρ T is given by i ∂ρT = − [Hˆ , ρT ]. ∂t ¯
(2)
The partial Wigner transform of ρ T is defined as 1 i T ρw (Q, P) = de− ¯ P · Q + /2|ρˆT |Q − /2, N (2π ¯) (3) where the integration is performed for the degrees of freedom that are to be treated classically. We note that ρw defined in Eq. (3) is still a density operator of the quantum subsystem. A direct derivation of the MQCL utilizes the following formally exact equation for the partial Wigner transform of ˆ 6, 34 the products of two operators Aˆ and B: ˆ w = Aw exp i¯ Bw = Bw exp − i¯ Aw , (Aˆ B) 2 2 (4) where Aw and Bw are the partial Wigner transforms of the opˆ is the Poisson bracket. For a pair of phase erators Aˆ and B, space distributions Aw (Q, P) and Bw (Q, P), is defined as ∂Aw ∂Bw ∂Aw ∂Bw . Aw Bw = {Aw , Bw } ≡ − ∂Qj ∂Pj ∂Pj ∂Qj j (5) By doing partial Wigner transform of the quantum Liouville equations in Eq. (2), and applying the classical limit ¯ → 0 to Eq. (4), the MQCL is given by29, 31, 32 i ∂ρw 1 = − [Hˆ w , ρw ] + [{Hˆ w , ρw } − {ρw , Hˆ w }], ∂t ¯ 2
(6)
where Hˆ w is the partial Wigner transform of the total Hamiltonian Hˆ . The above Eq. (6) is still in the operator form for the quantum subsystem. To perform the actual simulations, one will need to choose a basis set for the quantum degrees of freedom. Very often, equations of motion in the adiabatic basis are used.29, 31, 33 For the problems considered in this work, i.e., the Frenkel exciton model for molecular aggregates, and the case with a conical intersection, the calculations are performed in the diabatic basis. We note that extension of the current work to the adiabatic basis can be rather straightforward.
(1)
ˆ qˆ are the momentum and coordinate of the quantum where p, ˆ Q ˆ are those of the heavy particles, and μ and subsystem, P, ˆ is the coupling beˆ Q) Mj are the corresponding masses. Vˆ (q, tween the system and bath. For simplicity, only one degree of freedom for the quantum subsystem is considered. The MQCL equations can be derived through many different ways.28–32 For the reason of completeness, we briefly sketch the derivation using the partial Wigner transform of
B. Absorption spectra from the MQCL equation
In this subsection, we show how the MQCL equations can be applied to calculation of absorption spectra in a problem with several excited electronic states. To this end, we consider a model where the electronic degree of freedom is treated quantum mechanically, and the nuclear degrees of freedom are treated classically. The total molecular Hamiltonian plus the interaction with a classical
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
084105-3
Bai et al.
J. Chem. Phys. 140, 084105 (2014)
dynamics of the N coherence terms |mg| are decoupled from other terms of the density matrix including |gg|, |gm|, and |mn|. We further define ρwm (P, Q; t) = m|(μρg )w (P, Q; t)|g, the correlation function in Eq. (10) can then be calculated as N 1 C(t) = μ(t)μ(0)g = μm ρwm (Q, P; t). dPdQ (2π ¯)N m=1
electromagnetic (EM) field can be written as Hˆ = Hˆ g (P, Q)|gg| +
N
Hˆ m (P, Q)|mm|
m=1
+
N
Hˆ mn (P, Q)|mn|
m=n=1
− E(t) ·
N
μm (Q)[|mg| + |gm|].
m=1
In the above equation, |g and |m stand for the ground and the mth excited states, Q and P are the Cartesian coordinates and momenta of the molecule system, E(t) is the external EM field, and μm (Q) = μgm (Q) is the transition dipole operator, Hmn is the coupling between the mth and nth electronic states. The total Hamiltonian in Eq. (7) has neglected the coupling between ground and excited states, which is in general a good approximation and used in the previous works,10, 20 since the energy difference between ground and excited states is usually large.26, 35 When the EM field is weak, the absorption spectra can be calculated using the dipole-dipole correlation function defined as6 +∞ 1 I (ω) ∝ dteiωt μ(t) · μ(0)g , (8) 2π −∞ where μ = N m=1 μm (|gm| + |mg|) is the total dipole operator, the subscript g in the average indicates that the initial state to calculate the correlation function is equilibrated on the ground electronic state, ρ g = e−βHg /Tre−βHg (β = 1/kB T is the inverse temperature), and Hmol is the molecular Hamiltonian without the EM field, Hˆ mol = Hˆ g (P, Q)|gg| +
N
Hˆ mn (P, Q)|mn|,
(11)
(7)
(9)
Now we derive the equation of motion for ρwm . Using the general MQCL equation in Eq. (6), and since the coupling terms between ground state and excited states are zero, the dynamics of ρwm (Q, P) is governed by d m i ρ (Q, P; t) = − [Hw , ρw (Q, P; t)]mg dt w ¯ 1 ∂Hmn ∂ρwn ∂ρwn ∂Hmn + · − · 2 n ∂Q ∂P ∂Q ∂P 1 ∂Hg ∂ρwn ∂ρ n ∂Hg + · − w · . (12) 2 ∂Q ∂P ∂Q ∂P In the diabatic basis, the off-diagonal Hmn terms do P2 not depend on momentum. We also assume that Hg = 2M P2 + Vg (Q), Hm = 2M + Vm (Q). The N ρwm (Q, P; t) terms now form a closed equation, d m i ρ (Q, P; t) = − [Vm (Q) − Vg (Q)]ρwm (Q, P) dt w ¯ i − Hmn ρwn (Q, P) ¯ n=m 1 ∂Hmn ∂ρwn + {Hm , ρwm (P, Q)}+ · , 2 n=m ∂Q ∂P (13)
m,n=1
where we have used the notations Hˆ mm = Hˆ m when m = n. We then investigate how the absorption spectra can be calculated using the MQCL formalism. We work with the Condon approximation by assuming that the transition dipole moment does not depend on the nuclear degrees of freedom, while extensions to the non-Condon cases are straightforward. For simplicity, we will also drop the orientation dependence of μ in further derivations. The first step is to write the correlation function in Eq. (8) in the form of partial Wigner transformation, 1 dPdQμw (μρg )w (Q, P; t), μ(t)μ(0)g = TrS (2π ¯)N (10) N where μw = μ = m=1 μm (|gm|+|mg|), (μρg )w (Q, P; t) is the partial Wigner transform of the operator ˆ ˆ e−i H t/¯ (μˆ ρˆg )ei H t/¯ , which can be calculated using the MQCL approximation in Eq. (4). The initial state to the partial Wigner transform calculate N −βHg μ e /Tre−βHg |mg|. Since (μρg )w (P, Q; t) is m=1 m there is no direct coupling between the ground and excited states in the molecular Hamiltonian in Eq. (9), the
where the effective Hamiltonian Hm for the mth surface P2 is defined as Hm = 2M + V˜m (Q), with V˜m (Q) = 12 [Vg (Q) + Vm (Q)], and {Hm , ρwm (P, Q)} = −
∂ V˜m ∂ρwm P ∂ρwm · + · . M ∂Q ∂Q ∂P
(14)
The four terms on the right-hand side of Eq. (13) have different effects on the evolution of ρwm : (1) The first term gives a phase factor related to energy difference Vm (Q) − Vg (Q); (2) In the second term, the different ρwm components are coupled by the off-diagonal terms Hmn of the Hamiltonian; (3) The Poisson bracket in the third term indicates that the distribution ρwm (Q, P) evolves on the potential energy surface defined by V˜m , which is the average of the ground state potential energy surface Vg and the mth excited state Vm ; (4) The fourth term represents the “off-diagonal forces” originated from the derivative of the Hmn (Q) terms. Equation (13) consists of a set of partial differential equations for the Wigner distribution functions, which is generally hard to solve in the multi-dimensional case. However, if there is only one excited state, the correlation function can be calculated rather straightforwardly. In such cases, there is only
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
084105-4
Bai et al.
J. Chem. Phys. 140, 084105 (2014)
one term ρ = ρ eg , and the correlation function in Eq. (8) is 1 C(t) ≈ (15) dQdPμρw (Q, P; t); (2π ¯)N and the equation of motion for ρw (P, Q; t) is given by d i ρw (Q, P; t) = − [Ve (Q) − Vg (Q)]ρw (Q, P; t) dt ¯ + {Hav , ρw (Q, P; t)}, P2 2M
(16)
Ve (Q)+Vg (Q) 2
+ is the classical Hamiltonian where Hav = on the averaged potential energy surface. It is easy to find that t 1 i C(t) ≈ dQdPμρw (Q, P; 0)e− ¯ 0 dτ U (Q(τ )) , N (2π ¯) (17) where the dynamics is calculated on the average potential energy surface of the ground and excited states, and U(Q) is the difference between the excited state and ground state potential. This is just the WACL approach obtained in previous studies.20, 25 By neglecting the off-diagonal force terms in Eq. (13), and assuming that ρwm (P, Q; t) evolves on the ground state potential energy surface rather than V˜m , it can be shown that the correlation function in Eq. (8) can be calculated by the commonly used DCL method for multiple excited states11–14 ∞ I (ω) ∝ Re dte−iωt μm (0)Fmn (t)μn (t), (18) 0
mn
where the terms Fmn are calculated as dF (t) = iF (t) (t), (19) dt with the initial condition that Fmn (0) = δ mn , and the matrix is give by ¯ mn (t) = (Ei − Eg )(t)δmn + Hmn (t)(1 − δmn ).
(20)
The average is taken over the Boltzmann distribution of nuclear degrees of freedom on the ground state, and μ(t), (t) are now classical variables, whose dynamics are governed by the ground state nuclear Hamiltonian. C. The mean field approximation
As stated above, the goal is to derive a trajectory based method to solve Eqs. (11) and (13). Although this can be done easily for the single excited state problem as shown in Eq. (17), the same method can not be extended directly to the case of multiple coupled excited surfaces. In literature, there are several numerical algorithms developed to solve the MQCL equations based on surface-hopping type trajectory simulations.31, 33, 36–38 However, these methods usually suffer from large numerical noise at longer simulation times, which will cause noise problems in calculating the Fourier transform in Eq. (8). In this section, we derive a mean field approximation to solve the MQCL equation in Eq. (13). According to the previous work,28, 39, 40 the mean field approximation for the general MQCL equation in Eq. (6) is briefly presented in the Appendix. As shown below, the derivation of the mean
field approximation for the dipole-dipole correlation functions in Eq. (13) uses a similar idea, but with several technical differences. In the Appendix, the expectation values of classical variables Q and P are defined using ρ mm , which is the population on each state [see Eqs. (A2) and (A3)]. However, since in Eq. (13), only off-diagonal terms ρwm (P, Q; t) = m|(μρg )w (P, Q; t)|g are involved, a new definition of the mean field positions and momenta is needed. Fortunately, we can derive the following normalizing equation about the ρwm terms from Eq. (13): 1 ρwm ρwm∗ = C. (21) dPdQ N (2π ¯) m To develop a method with trajectory description, the following integral of ρwm (Q, P; t) is defined: ρ¯m (t) = dQdPρwm (Q, P; t), (22) and by using Eq. (21), we define the averages of the positions Q and momenta P as 1 ¯ = ρwm (Q, P)ρwm∗ (Q, P)Q, (23) Q dQdP (2π ¯)N C m P¯ =
1 (2π ¯)N C
dQdP
ρwm (Q, P)ρwm∗ (Q, P)P,
(24)
m
where the constant C is defined in Eq. (21). So the difference between the current derivation and the conventional mean field approximation in the Appendix is that, instead of using the averages of Q and P defined using the diagonal elements of the density matrix, Eqs. (23) and (24) have used a “wave function” like definition of the averages. The rest of the derivation is then very similar to the conventional mean field approximation. We use the following ansatz for ρ j (Q, P): ¯ ¯ P − P(t)). ρwm (Q, P) = ρwm (t)δ(Q − Q(t),
(25)
Together with Eq. (13) and the definitions in Eqs. (22)– (24), the equations of the motion which constitute the mean field approximation to calculate the dipole-dipole correlation functions can be obtained as follows: i i d ¯ ¯ ¯ ρ¯n , (26) ρ¯m = − [Vm (Q)−V Hmn (Q) g (Q)]ρ¯m − dt ¯ ¯ n=m and d ¯ P¯ , Q= dt M −1 ¯ d ¯ 1 ∂Vg (Q) 1 2 |ρ¯m | − P = Fmf = − ¯ dt 2 ∂Q 2 m
¯ ∂Hmn (Q) . × Re ρ¯m∗ ρ¯n ¯ ∂Q n
(27)
(28)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
084105-5
Bai et al.
J. Chem. Phys. 140, 084105 (2014)
D. The modified mean field method
The shortcomings of the mean field dynamics have been widely investigated.3 Now we work on possible extensions of the mean field approach. We note that a major problem of the mean field ansatz is that, in the multiple surface case, trajectories associated with different ρwm (Q, P) should propagate on different potential surfaces (see Eq. (14)). The mean field ansatz then becomes problematic, since even though the initial positions and momenta are the same on different surfaces, they soon become different evolving under the influence of different forces. Nevertheless, we can still take the mean field dynamics as a reference, and investigate further modifications. To derive the modified mean field method, we assume that each ρwm term is slightly different from the mean field ansatz in Eq. (25). So we can still use the mean field trajec¯ ¯ as the zeroth order approximatory defined by Q(t) and P(t) tion, and obtain the correction terms though expansions of the phase space variables P, Q, and Hmn (P, Q), near the reference ¯ ¯ trajectory Q(t), P(t). The following terms are defined as deviations of the Q and P expectation values from the mean field ansatz: ¯ m (Q, P), Xm = dQdP(Q − Q)ρ (29) w m =
¯ m (Q, P). dQdP(P − P)ρ w
(30)
In this study, we only include the first order moments for simplicity, and high order moments are not considered. Several other expansions of the phase space average are needed in deriving the modified mean field equation. Calculation of the dQdPHmn (Q)ρwn (Q, P) term can be approximated as dQdPHmn (Q)ρwn (Q, P) ≈
¯ n (Q, P) dQdPHmn (Q)ρ w
¯ ∂Hmn (Q) ¯ n (Q, P) (Q − Q)ρ w ¯ ∂Q ¯ ¯ ρ¯n + ∂Hmn (Q) · Xm . = Hmn (Q) (31) ¯ ∂Q Terms like dQdPHmn (Q)ρwn (Q, P)Q can be approximated in the following way: dQdPHmn (Q)ρwn (Q, P)Q +
dQdP
≈
¯ Qρ ¯ n (Q, P) dQdPHmn (Q) w +
¯ ¯ n (Q, P) ¯ Q ¯ ∂Hmn (Q) (Q− Q)ρ dQdP Hmn (Q)+ w ¯ ∂Q
¯ Q ¯ ρ¯n + Hmn (Q)X ¯ m+Q ¯ = Hmn (Q)
¯ ∂Hmn (Q) · Xm . ¯ ∂Q
(32)
Using the approximations similar to the ones in the above two equations, further derivation of the modified mean field
equations can be obtained by taking proper phase space averages over Eq. (13). For example, the dynamics for ρ¯m is obtained by combining the integration over Eqs. (13) and (31): i d ¯ − Vg (Q)] ¯ ρ¯m − i ¯ ρ¯n ρ¯m = − [Vm (Q) Hmn (Q) dt ¯ ¯ n=m
¯ − Vg (Q)] ¯ i ∂[Vm (Q) − ¯ ¯ ∂Q
¯ i ∂Hmn (Q) · Xn , · Xm − ¯ ¯ n=m ∂Q
(33)
which adds correction terms involving X to the original mean field result in Eq. (26). Equations (27) and (28) define the reference mean field trajectory and stay the same. The following equation of motion for Xm can be derived by multiplying Q to Eq. (13) and doing the integration, and utilizing Eq. (32): i d m i ¯ − Vg (Q)]X ¯ ¯ n. Xm = − [Vm (Q) Hmn (Q)X m− dt M ¯ ¯ n=m (34) The equation of motion for involves second order derivatives of Hmn (Q), which is given by d i i ¯ − Vg (Q)] ¯ m = − [Vm (Q) Hmn n − ρ¯m Fmf m− dt ¯ ¯ n=m −
¯ + Vg (Q)] ¯ ¯ 1 ∂[Vm (Q) 1 ∂Hmn (Q) ρ¯m − ρ¯n ¯ ¯ 2 2 n=m ∂ Q ∂Q
−
¯ ¯ ¯ 1 ∂ 2 [Vm (Q)+V 1 ∂ 2 Hmn (Q) g (Q)] ·X − ·Xn , m 2 2 ¯ ¯ 2 2 n=m ∂ Q ∂Q (35)
where Fmf is the expression of force in mean field approximation, as in Eq. (28); and order derivative terms the second mn (Q) m ρw (Q, P). come from integral like dQdP ∂H∂Q We note that this modified mean field method is similar to the modified Ehrenfest dynamics obtained by Horsfield and co-workers,41, 42 which is also an equation of motion for first order moments, and has been used for calculations of coupled electron-ion dynamics in inelastic quantum transport. III. RESULTS
We now apply the above methods to calculate the dipoledipole correlation functions and absorption spectra of several model systems, which include (1) a model of molecular aggregate dimers; (2) a model of the FMO complex; and (3) models of the pyrazine molecule with S1 -S2 conical intersection, including the 4-mode and 24-mode models. The results calculated using the mean field and modified mean field methods are compared with benchmark calculations using the Hierarchical equations of motion (HEOM) method,43, 44 and the MCTDH method.26, 45, 46 In the simulations, Eqs. (26), (27), and (28) are solved for the mean field method. The initial positions and momenta of the trajectory are sampled from the Wigner distribution on
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 15:04:10
Bai et al.
J. Chem. Phys. 140, 084105 (2014)
A. Absorption spectra of molecular aggregates
We consider in this subsection the Frenkel exciton model of N two-level molecules coupled to a phonon bath, where the total Hamiltonian is written as H = He + Hph + He−ph ,
(36)
where He describes the electronic degrees of freedom, Hph describes the nuclear degrees of freedom, and He−ph describes their couplings. To calculate the absorption spectra, only the ground state and the single excited states are considered: He = 0 |00|+
N
m |mm|+
N N
Jmn (|mn|+|nm|),
m=1 n
