Electronically Nonadiabatic Dynamics in Complex Molecular Systems

Article pubs.acs.org/JPCA

Electronically Nonadiabatic Dynamics in Complex Molecular Systems: An Efficient and Accurate Semiclassical Solution Guohua Tao* School of Advanced Materials, Peking University, Shenzhen Graduate School, Shenzhen, China 518055 S Supporting Information *

ABSTRACT: Chemical reaction dynamics is always a central theme in chemistry research. In many important chemical processes, reaction dynamics is electronically nonadiabatic, i.e., dynamics involves coupled multiple electronic states. We demonstrate in this paper that a semiclassical (SC) treatment based on an initial value representation methodology and a classical mapping formalism for the electronic degrees of freedom is now able to provide a rigorous and practical solution to electronically nonadiabatic dynamics in complex molecular systems. The key component of this treatment is to incorporate a correlated importance sampling protocol in nonadiabatic SC calculations, which results in a speedup factor of 100 or more in comparison with that using the standard sampling approach. This is illustrated by application to a two-state model coupled with up to 10 nuclear bath modes for a benchmark nonadiabatic excitation energy transfer problem. This work provides great opportunities for the effectively theoretical investigations on reaction mechanisms in complex molecular systems, in which electronically nonadiabatic dynamics plays an importance role.

■

methodology15,16 integrates both electronic and nuclear motions into a consistent dynamical framework, and quantum coherent dynamics is inherently incorporated.17 The main challenge is how this SC methodology can be implemented efficiently to electronically nonadiabatic dynamics in complex molecular systems. In this paper we report how the recently developed correlated importance sampling method18 makes this SC solution to high-dimensional nonadiabatic dynamics practical. In comparison with the standard sampling method, a speedup factor of 100 or more has been achieved in the new nonadiabatic SC calculations for a benchmark two-state system coupled with up to 10 nuclear DOF to model the excitation energy transfer in photosynthesis systems.19 This model photosynthetic system has been studied by a variety of theoretical methods including the reduced hierarchy equation method,20 linearized SC-IVR (LSC-IVR),21 partial linearized density matrix approach,22 and other methods which are mostly derived from density matrix formalism.23−25

INTRODUCTION The well-known Bohn−Oppenheimer (BO) approximation,1 which assumes that nuclei evolve along a single potential energy surface that is defined by fast-moving electrons, is a milestone in the roadmap for understanding chemical reactions at a molecular level. However, in many fields such as photochemical reactions,2 nanoscale electronics,3 and chemical dynamics in biological systems,4,5 BO approximation breaks down and nonadiabatic transitions occur between different electronic states, or different BO potential energy surfaces.6,7 Nonadiabatic effects here play a crucial role in characterizing many important features in chemical reactions, such as reaction rates and branching ratios. An ideal general theoretical method for treating electronically nonadiabatic chemical dynamics should be able to treat both electronic and nuclear dynamics consistently. However, in prevalent mixed quantum-classical simulations,8−10 only a few degrees of freedom (DOF) are treated quantum mechanically while nuclear dynamics is usually treated classically. The inconsistency between the quantum and classical treatments may results in unphysical or inaccurate description of the dynamics, e.g., Ehrenfest dynamics may predict a wrong asymptotic behavior,11 and a priori knowledge of the quantum decay time is normally needed in the surface hopping algorithms.12 Alternatively the classical mapping formalism, i.e., the Meyer−Miller−Stock−Thoss (MMST) theory,13,14 exactly maps discrete quantum electronic states into continuous classical degrees of freedom, i.e., each electronic state is mapped into a classical harmonic oscillator. Combined with the MMST theory, a semiclassical (SC) approach based on an initial value representation (IVR) © 2013 American Chemical Society

■

THEORETICAL METHODS

To proceed, we start from the theoretical foundation of the MMST and SC-IVR methods. In terms of Cartesian coordinates and momenta, the MMST nuclear-electronic Hamiltonian is given by Received: May 17, 2013 Revised: June 22, 2013 Published: June 25, 2013 5821

dx.doi.org/10.1021/jp404856p | J. Phys. Chem. A 2013, 117, 5821−5825

The Journal of Physical Chemistry A

Article

n

H(x , p, Q , P) =

∑ 1 (xi 2 + pi 2 2

i=1

n

+

C LSC(t ) = (2π ℏ)−F

− 1)Hii(P, Q)

(6a)

n

∑ ∑

(xixj + pp )Hij(P, Q) i j

where Aw and Bw are the Wigner functions of the corresponding operators, defined by

(1)

i=1 j=i+1

Ow (q 0, p0) =

where n is the number of electronic states, Hii and Hij are the diabatic Hamiltonian matrix elements (which are in general functions of nuclear coordinates and momenta Q, P), and xi and pi are the Cartesian variables for the ith electronic DOF, respectively. Classical trajectories for the nuclear and electronic DOF are generated, as usual, by Hamilton’s equations ∂H , Q̇ = ∂P

∂H Ṗ = − , ∂Q

ẋ =

∂H , ∂P

ṗ = −

∂H ∂x

ρMC (p0 , q 0, p′0 , q′0 ) ∝ |⟨p0 , q 0|Âs |p′0 , q′0 ⟩⟨P0̅ Q̅ 0|Âb |P0̅ Q̅ 0⟩|

∫ dq 0 ∫ dp0Ct(p0, q 0)eiS (p ,q )/ℏ|pqt t⟩ t

0

0

⟨p0q 0|

(4a)

where F is the DOF of the system, |p0q0⟩ and |ptqt⟩ are the coherent states with (p0,q0) and (pt,qt) the initial and timeevolved classical phase points, respectively, St(p0,q0) is the classical action (the time integral of the Lagrangian) along the classical path, and Ct(p0,q0) a pre-exponential factor involving the monodromy matrix elements

ρMC (p0 , q 0, p′0 , q′0 ) ∝ |⟨p0 , q 0|Â e |p′0 , q′0 ⟩⟨P0̅ Q̅ 0|Â n |P0̅ Q̅ 0⟩|

− iℏγ

∂q t ∂p0

γ

1/2

⎤1/2 i −1/2 ∂pt −1/2⎞⎥ ⎟⎟ + γ γ ℏ ∂q 0 ⎠⎥⎦ (4b)

2⎤ ⎛ Pki 2 cki ⎞ ⎥ 1 ⎟ + mkiωki 2⎜Q ki − 2 mkiωki 2 ⎠ ⎥⎦ ⎝ ki = 1 ⎣ 2mki n N ⎡ ⎤ Pkj 2 1 +∑∑⎢ + mkjωkj 2Q kj 2 ⎥ ⎥⎦ 2 ⎣ 2mkj j ≠ i kj = 1 ⎢ N

By substituting eq 4 for the two time evolution operators in eq 3 we obtain the full SC-IVR correlation function as the following double phase space average CSC(t ) = (2π ℏ)−2F

Hii = εi 0 +

∫ dq 0dp0 ∫ dq′0 dp′0 Ct(p0, q 0)Ct(p′0

, q′0 )*ei[St(p0, q 0) − St(p ′0 , q ′0 )]/ ℏ⟨p0q 0|Â |p′0 q′0 ⟩ ⟨p′t q′t |B|̂ pq ⟩ t t

(8)

where Â e and Â n are the factorized electronic and nuclear part of operator Â , the phase point (p0, q0) includes both electronic and nuclear DOF. We demonstrate in the following how this correlated importance sampling protocol provides an efficient and accurate SC solution to high-dimensional nonadiabatic dynamics by considering a benchmark photosynthetic model system.20 The Frenkel exciton Hamiltonian for this model system is given by

⎡ ⎛ ∂p ∂q 1 Ct(p0 , q 0) = det⎢ ⎜⎜γ 1/2 t γ −1/2 + γ −1/2 t γ 1/2 ⎢⎣ 2 ⎝ ∂q 0 ∂p0 1/2

(7)

,where P̅0 = (P0 + P′0)/2, Q̅ 0 = (Q0 + Q′0)/2 are the transformed bath DOF and here we assume the operator Â can be factorized into a product of a system part and a bath part. Note the phase point (p0,q0) includes both the system and the bath DOF, and the path correlation between the two initial phase points for the bath DOF is incorporated in ⟨p0,q0|Â s| p′0,q′0⟩. In a separate paper,30 we show that how this timeindependent correlated importance sampling approach performs in comparison with the previously proposed timedependent importance sampling method.31 We now generalize the correlated importance sampling to electronically nonadiabatic dynamics, by treating the electronic-nuclear coupling as a system-bath model. The correlated importance sampling function for the SC-IVR calculations of nonadiabatic dynamics is therefore given by

(3)

here Â and B̂ are operators for general quantities of interest. Considering the extreme difficulty in the exact quantum calculations for complex systems, approximations are normally made to eq 3. The SC-IVR theory approximates the time evolution operators by a phase space average over the initial phase points of classical paths (trajectories); e.g., the coherent state IVR due to Herman and Kluk27,28 is ̂

0

where O = A, B. Although LSC-IVR has been applied to systems with hundreds of DOF,21 the more accurate SC-IVR was previously thought to work only for a few DOF in practice. To facilitate the implementation of SC-IVR, we proposed a very efficient correlated importance sampling approach,18 which introduces the path correlation in the sampling function of the initial phase points:

(2)

̂ ̂ CQM(t ) = tr[Â eiHt / ℏBê −iHt / ℏ]

∫ dΔqe−iΔq·p ⟨q 0 + Δq/2|Ô|q 0 − Δq/2⟩ (6b)

with the Hamiltonian given by eq 1. These are the equations of motions in the SC-IVR calculations of dynamical properties. The dynamical properties of complex molecular systems can be evaluated, in principle, by the corresponding quantum mechanical time correlation function26

e−iHt / ℏ = (2π ℏ)−F

∫ dq 0 ∫ dp0Aw (q 0, p0)Bw (q 0, p0)

⎡

∑ ⎢⎢

Hij = Jij = constant (9)

(5)

here n and N are the total number of bacteriochlorophyll sites and the number of the bath modes coupled with each site, respectively; mki, Pki, Qki ωki, cki are the mass, momentum, coordinate, frequency, and coupling strength, respectively, for

over the initial phase points of a pair of paths. Further linearized approximation to the SC-IVR leads to LSC-IVR or the classical Wigner model for time correlation functions:29 5822



Article

the kth harmonic phonon bath mode coupled with ith site. Note here each electronic state couples with independent nuclear bath, and we only consider a two-site model in this work. The energy levels of the two sites are ε10 = 100 cm−1, ε20 = 0 cm−1, the site−site electronic coupling J12 = 100 cm−1, and the bath temperature is 300 K. The spectra density of phonon modes (for all sites) takes the Debye form, ρ(ω) = 2λ

ωτc−1 ω 2 + τc−2

(10a)

where the reorganization energy λ represents the electronicnuclear coupling strength, and τc = 500 fs is the characteristic time of the phonon bath. We discretize the continuous phonon spectra density eq 10a as done in ref 32. In brief, we have ⎛ arctan(ωm /ωc) ⎞ j⎟ ωj = ωc tan⎜ ⎝ ⎠ N cj = ωj

4λarctan( −ωm /ωc) πN

(10b)

Figure 1. Excitation energy transfer dynamics of a two-state model coupled with uncorrelated harmonic phonon bath. The population of each site (solid lines for site 1 and dashed lines for site 2) is calculated by using the standard and the correlated importance sampling methods: (a) no bath, (b) 2 bath modes for each state, (c) 5 bath modes for each state.

where N is the number of bath modes, ωc is the characteristic frequency of the bath, and ωm = 100ωc is the maximum frequency. The electronic excitation population functions for each site can be evaluated in terms of time correlation function by eqs 3−6, ̂

̂

becomes less pronounced. In panel b there are 2 bath modes coupled to each electronic state, the correlated importance sampling approach provides nicely converged results, while even with 10 times more trajectories used, the standard sampling method performs poor at longer times. By examining the standard deviations of these calculations in Figure 2, it appears that 109 trajectories in the standard sampling protocol result in a comparable standard deviation with that from 107 trajectories with use of the correlated sampling method, which indicates that 2 orders of more

̂

−βHb iHt / ℏ ̂ Pi(t ) = tr[ρ (0)e e ρiî e−iHt / ℏ]

with the operators Â = ρ̂(0)e matrix element is

(11)

−βĤ b

, B̂ = ρ̂ii. Here the density

ρ̂ii = |Φ⟩⟨Φ| i i

(12a)

where the electronic-oscillator wave function for electronic state i is given by n

Φi(x) = ⟨x|Φ⟩ i =

2 ⎛ 1 ⎞n /4 ⎜ ⎟ 2 xi ∏ e−1/2xj ⎝π ⎠

j=1

(12b)

and the corresponding coherent state wave function is n

Φi(x , p) = ⟨x , p|Φ⟩ i =

2 2 2 (xi − ipi ) ∏ e−1/4(xj + pj − 2ixjpj ) 2 j=1

(12c)

The nonadiabatic dynamics is integrated by using the fourth order Runge−Kutta method with a time step dt = 50.0 au. Monodromy matrices are also integrated in the same way, from which the SC prefactor is calculated by using eq 4.

■

RESULTS AND DISCUSSIONS Figure 1 illustrates the excitation energy transfer calculated by SC-IVR with use of different sampling methods. For comparison, we also show in panel a the results for the pure electronic two-state model, i.e., with no bath, by using the SCIVR (SC), and the quantum discrete variable representation method33 (QM). Both methods provide identical transfer dynamics. In the SC-IVR calculations, 108 trajectories are used for the correlated importance sampling protocol and 109 for the standard sampling method (105 for the no bath case). The electronic-nuclear coupling strength λ = J12. As more bath modes are introduced, the Rabi oscillation (coherent dynamics)

Figure 2. Same system as in Figure 1, but comparing the results of the standard deviation of the population functions (site 1 only) in the SCIVR calculations. The numbers of trajectories used are given in the legends. 5823



Article

trajectories are required for the standard sampling to reach a comparable convergence with the correlated sampling. As the system becomes more complex, for the 5 bath modes per state case, the advantages of using the correlated sampling approach are more appreciated: the results from the standard sampling method with use of even 109 trajectories are far from converged. The corresponding standard deviation calculations in Figure 2 imply that the speedup factor of using the correlated sampling is about 1000. The high-dimensional nonadiabatic calculations were also performed for a variety of reorganization energies (electronicnuclear coupling strength). Figure 3 illustrates the effect of

Figure 4. Same system as in Figure 3, but comparing the standard deviation of the population functions (site 1 only) by using the standard and the correlated importance sampling methods.

effect of electronic-nuclear coupling strength on the computational efficiency is also investigated. The success of this SC treatment for a benchmark two-state model coupled with up to 10 nuclear DOF paves the way for future work on exploring chemical reaction mechanisms in complex molecular systems, i.e., in clusters, solutions, or even materials and biological systems.

■

Figure 3. The site population functions for different electronic-nuclear coupling strengths (reorganization energy). The results are for the 5 bath modes/state case using the correlated importance sampling. (a) λ = J12/50, (b) λ = J12/5, and (c) λ = J12.

ASSOCIATED CONTENT

S Supporting Information *

The results of the standard deviation of the population functions for site 2. This material is available free of charge via the Internet at http://pubs.acs.org.

■

electronic-nuclear coupling strength on the transfer dynamics for the 5 bath modes per state case. As the electronic-nuclear coupling increases, it shows a greater degree of quantum decoherence (the magnitude of Rabi oscillation becomes smaller). At the same time, the coupled Hamiltonian dynamics becomes highly nonlinear therefore the performance of the SCIVR calculations degrades a bit, especially at long times, as shown in Figure 4 by comparing the standard deviations of the correlation function for different coupling strength parameters in the SC-IVR calculations. Figure 4 demonstrates clearly again the advantage of using the correlated importance sampling over the standard one.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interests.

■

ACKNOWLEDGMENTS This work was supported by Peking University Shenzhen Graduate School, Shenzhen Engineering Lab of Thin Film Materials of Solar Cells, and the Shenzhen Strategic Emerging Industry Development Funds (Grant No. JCYJ20120829170028565).

■

■

CONCLUSIONS Electronically nonadiabatic dynamics involves the evolution of nuclear motions on coupled multiple BO potential energy surfaces. In this work we show that with the aid of a correlated importance sampling approach, the SC-IVR in the MMST framework provides an efficient and accurate solution to electronically nonadiabatic dynamics in complex molecular systems. The inclusion of path correlation in the bath DOF in the sampling protocol requires 100 or more times fewer trajectories than that using the standard sampling method to obtain converged results in nonadiabatic SC calculations. The

REFERENCES

(1) Born, M.; Oppenheimer, E. Zur Quantentheorie der Molekeln. Ann. Phys. 1927, 84, 457−484. (2) Zewail, A. H. Femtochemistry: Atomic-Scale Dynamics of the Chemical Bond. J. Phys. Chem. A 2000, 104, 5660−5694. (3) Nitzan, A.; Ratner, M. A. Electron Transport in Molecular Wire Junctions. Science 2003, 300, 1384−1389. (4) Marcus, R. A.; Sutin, N. Electron Transfers in Chemistry and Biology. Biochim. Biophys. Acta 1985, 811, 265−322. (5) Nagel, Z. D.; Klinman, J. P. Tunneling and Dynamics in Enzymatic Hydride Transfer. Chem. Rev. 2010, 110, PR41−PR67.

5824



Article

(6) Jasper, A. W.; Zhu, C.; Nangia, S.; Truhlar, D. G. Non-Adiabatic Effects in Chemical Dynamics. Faraday Discuss. 2004, 127, 1−22. (7) Tully, J. C. Perspective: Nonadiabatic Dynamics Theory. J. Chem. Phys. 2012, 137, 22A301. (8) Billing, G. D. On the Applicability of the Classical Trajectory Equations in Inelastic Scattering Theory. Chem. Phys. Lett. 1975, 30, 391−393. (9) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061−1071. (10) Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab Initio Ehrenfest Dynamics. J. Chem. Phys. 2005, 123, 084106. (11) Miller, W. H. Electronically Nonadiabatic Dynamics via Semiclassical Initial Value Methods. J. Phys. Chem. A 2009, 113, 1405−1415. (12) Schwartz, B. J.; Bittner, E. R.; Prezhdo, O. V.; Rossky, P. J. Quantum Decoherence and the Isotope Effect in Condensed Phase Nonadiabatic Molecular Dynamics Simulations. J. Chem. Phys. 1994, 104, 5942−5955. (13) Meyer, H. D.; Miller, W. H. A Classical Analog for Electronic Degrees of Freedom in Nonadiabatic Collision Processes. J. Chem. Phys. 1979, 70, 3214−3223. (14) Stock, G.; Thoss, M. Semiclassical Description of Nonadiabatic Quantum Dynamics. Phys. Rev. Lett. 1997, 78, 578−581. (15) Miller, W. H. Spiers Memorial Lecture: Quantum and Semiclassical Theory of Chemical Reaction Rates. Faraday Discuss. 1998, 110, 1−21. (16) Miller, W. H. The Semiclassical Initial Value Representation: A Potentially Practical Way for Adding Quantum Effects to Classical Molecular Dynamics Simulations. J. Phys. Chem. A 2001, 105, 2942− 2955. (17) Miller, W. H. Perspective: Quantum or Classical Coherence? J. Chem. Phys. 2012, 136, 210901. (18) Pan, F.; Tao, G. Importance Sampling Including Path Correlation in Semiclassical Initial Value Representation Calculations for Time Correlation Functions. J. Chem. Phys. 2013, 138, 091101. (19) Engel, G. S.; Calhoun, T. R.; Read, E. L.; Ahn, T. K.; Mancal, T.; Cheng, Y. C.; Blankenship, R. E.; Fleming, G. R. Evidence for Wavelike Energy Transfer through Quantum Coherence in Photosynthetic Systems. Nature 2007, 446, 782−786. (20) Ishizaki, A.; Fleming, G. R. Unified Treatment of Quantum Coherent and Incoherent Hopping Dynamics in Electronic Energy Transfer: Reduced Hierarchy Equation Approach. J. Chem. Phys. 2009, 130, 234111. (21) Tao, G.; Miller, W. H. Semiclassical Description of Electronic Excitation Population Transfer in a Model Photosynthetic System. J. Phys. Chem. Lett. 2010, 1, 891−894. (22) Huo, P.; Coker, D. F. Partial Linearized Density Matrix Dynamics for Dissipative, Non-Adiabatic Quantum Evolution. J. Chem. Phys. 2011, 135, 201101. (23) Jang, S. Theory of Multichromophoric Coherent Resonance Energy Transfer: A Polaronic Quantum Master Equation Approach. J. Chem. Phys. 2011, 135, 034105. (24) Kelly, A.; Rhee, Y. M. Mixed Quantum-Classical Description of Excitation Energy Transfer in a Model Fenna−Matthews−Olsen Complex. J. Phys. Chem. Lett. 2011, 2, 808−812. (25) Berkelbach, T. C.; Markland, T. E.; Reichman, D. R. Mixed Quantum-Classical Description of Excitation Energy Transfer in a Model Fenna−Matthews−Olsen Complex. J. Chem. Phys. 2012, 136, 084104. (26) Berne, B. J.; Harp, G. D. On the Calculation of Time Correlation Functions. Adv. Chem. Phys. 1970, 17, 63−227. (27) Herman, M. F.; Kluk, E. A Semiclasical Justification for the Use of Non-Spreading Wavepackets in Dynamics Calculations. Chem. Phys. 1984, 91, 27−34. (28) Kluk, E.; Herman, M. F.; Davis, H. L. Comparison of the Propagation of Semiclassical Frozen Gaussian Wave Functions with Quantum Propagation for a Highly Excited Anharmonic Oscillator. J. Chem. Phys. 1986, 84, 326−334.

(29) Sun, X.; Miller, W. H. Forward−Backward Initial Value Representation for Semiclassical Time Correlation Functions. J. Chem. Phys. 1999, 110, 6635−6644. (30) Tao G. Importance Sampling in Semiclassical Initial Value Representation Calculations for Time Correlation Functions: TimeDependent or Time-Independent? Submitted to J. Chem. Phys. (31) Tao, G.; Miller, W. H. Time-Dependent Importance Sampling in Semiclassical Initial Value Representation Calculations for Time Correlation Functions. J. Chem. Phys. 2011, 135, 024104. (32) Wang, H.; Thoss, M.; Sorge, K. L.; Gelabert, R.; Giménez, X.; Miller, W. H. Semiclassical Description of Quantum Coherence Effects and Their Quenching: A Forward−Backward Initial Value Representation Study. J. Chem. Phys. 2001, 114, 2562−2571. (33) Colbert, D. T.; Miller, W. H. A Novel Discrete Variable Representation for Quantum Mechanical Reactive Scattering via the SMatrix Kohn Method. J. Chem. Phys. 1992, 96, 1982−1991.

5825


Electronically Nonadiabatic Dynamics in Complex Molecular Systems

Recommend Documents