LETTER pubs.acs.org/JPCL
Mixed Quantum-Classical Description of Excitation Energy Transfer in a Model FennaMatthewsOlsen Complex Aaron Kelly* and Young Min Rhee* Department of Chemistry, Pohang University of Science and Technology (POSTECH), Pohang 790-784, Korea
bS Supporting Information ABSTRACT: The FennaMatthewsOlsen (FMO) complex has recently become a paradigmatic model system in terms of understanding the long-lived electronic quantum coherence that has been experimentally observed in photosynthetic systems. In this article we investigate the quantum dynamics of an FMO model within a mixed quantum-classical dynamics approach known as the Poisson bracket mapping equation (PBME), and explore the consequences of adopting this approximate description by simulating population transfer and electronic coherence. The results obtained via the Poisson bracket mapping formalism are explicitly compared with a selection of recent results on the FMO complex. The PBME is shown to be in excellent agreement with benchmark computational results at physiological temperature, and thus provides a computationally efficient and physically consistent algorithm that may be easily integrated in all-atom molecular dynamics simulations. SECTION: Dynamics, Clusters, Excited States
P
hotosynthesis typically involves electronic excitation energy transfer (EET) between many spatially separated chlorophyll molecules that are embedded within a protein matrix. In a large class of systems, this process occurs at ambient temperature. Hence, it is quite remarkable that recent experiments have shown that EET in photosynthetic (and other light-harvesting) systems exhibits coherent or wave-like characteristics that persist under such conditions for hundreds of femtoseconds.13 Reliable and tractable approximate quantum theories, and associated numerical simulation methods, are thus needed in order to understand the underlying effects that preserve the delicate quantum coherence exhibited by these systems, and to assess what role, if any, quantum coherence plays in the overall function of the photosynthetic process. In this letter we consider a quantum dynamics scheme described by an approximate form of the quantum-classical Liouville equation (QCLE)4 cast in the mapping basis, known as the Poisson bracket mapping equation (PBME).5,6 The QCLE describes the dynamics of a quantum subsystem coupled to a classical bath, and provides a basis for the construction of methods to simulate quantum evolution in large complex systems of this type. As well, it has been shown to yield an accurate description of the dynamics for many condensed phase systems that can be approximated as a subset of quantum degrees of freedom interacting with an almost-classical environment.7,8 In fact, the QCLE is an exact quantum description in the case of an arbitrary (n-level) quantum subsystem bilinearly coupled to a harmonic environment. The mapping basis formalism itself has been used by a number of r 2011 American Chemical Society
researchers in the development of various approximate quantum mechanical treatments over the years,914 and it entails an exact mapping of discrete quantum states onto continuous variables. In quantum-classical systems, such a mapping leads to phase-spacelike evolution equations for both quantum and classical degrees of freedom, and consequently its simulation is simple and given in terms of Newtonian trajectories. Comparisons of the solutions of the PBME with exact quantum results for the spin-boson model and other simple curve crossing models have shown that it provides a quantitatively accurate description of the dynamics for these systems.5,6 In this letter we explore the dynamics prescribed by the PBME in a model for the seven-state FennaMatthews Olsen (FMO) complex.1517 In order to accurately model quantum dynamics in systems of the photosynthetic type, one requires a physical theory which lies beyond the perturbative and Markovian regimes.18 The fact that the environmental coupling strengths are of the same order of magnitude as the intersite couplings in these systems negates the possibility of applying a perturbative approach, and the relaxation time of the environment is too long to justify a Markovian approximation. Additionally, in order to treat more realistic systems, a theory which may be easily extended to all-atom molecular dynamics (MD) simulations is also desirable. Indeed, Received: January 13, 2011 Accepted: March 10, 2011 Published: March 16, 2011 808
dx.doi.org/10.1021/jz200059t | J. Phys. Chem. Lett. 2011, 2, 808–812
The Journal of Physical Chemistry Letters
LETTER
where {Am,Bm(t)}x,X denotes a Poisson bracket in the full mapping-bath phase space of the system. By turning to an approximate form of QCL dynamics5 that neglects the second term on the right of eq 4, the PBME is obtained:
the QCLE and PBME fulfill all of these above requirements. Since the PBME requires significantly less computational effort to simulate than the QCLE, it is a much more attractive option in terms of extensions to larger numerical simulations. However, since the PBME is an approximation to the full QCLE, its validity, properties and applicability to various model systems must first be carefully assessed. Interestingly, the PBME formalism leads to the same set of evolution equations as the linearized semiclassical initial value representation (LSC-IVR) of Miller and co-workers.11,12 However, as these formalisms are derived in very different manners, and numerical simulations of the observable quantities also contain various differences, the validity check for the PBME method should be performed separately from the LSC-IVR case.19 Explicit comparisons of these two different approaches could also be useful in determining the relative strengths and weaknesses of each method. We also offer comparisons of the PBME simulation results with corresponding outcomes from other nonperturbative, non-Markovian theories17,20 in order to assess the quality of the PBME approximation. From this, we find that the PBME simulations yield excellent agreement with benchmark results for the FMO complex at room temperature. In performing mixed quantum-classical calculations, we are interested in computating observables such as electronic state populations and coherence as a function of time, for systems which can be partitioned into a quantum subsystem and a classical-like environment. When projected onto the mapping basis, the expression for the expectation value of a general ^W(R,P) over an initial nonequilibrium density matrix observable B F^W(X) takes the following form:6 Z BðtÞ ¼ dX dxBm ðx, X, tÞ~ F m ðx, XÞ ð1Þ
d Bm ðx, X, tÞ ¼ iL 0m Bm ðx, X, tÞ dt
Here, the approximate Liouville operator denotes iL 0m= {Hm, 3 }x,X. The neglected term has been shown to be equal to one-quarter of the back reaction of the quantum subsystem on the bath, which vanishes in the case of adiabatic dynamics.6 For observables that do not depend explicitly on the environmental momenta, the error induced by neglecting this part of the Liouville operator enters the time-dependent quantity indirectly through the density, and accrues as time proceeds.6 Hence, one expects that the PBME formulation may indeed capture shorttime dynamics accurately in instances where the back-reaction term is small, but that longer time-scale phenomena, such as the establishment of thermal equilibrium, may not be accurately obtained. Although estimates of the influence of this term on average values can be made, no easily implementable and rigorous simulation algorithms for the full QCLE in the mapping basis, eq 4, have yet been constructed. The task of computing time-dependent average values using ~m(x,X), and the PBME consists of two parts: sampling from F evaluating the dynamics governed by the classical-like propaga0 tor, eiL mt. In order to simulate the dynamics generated by iL 0m, an integration scheme motivated by the separation of time-scales between the electronic and nuclear motions in the problem was employed together with a Trotter-type decomposition of the system propagator. The integration method is time-reversible, symplectic, and provides an analytic solution for the quantum subsystem degrees of freedom (unpublished results). The initial ~m(x,X) and the integration scheme mentioned sampling of F above are both key differences between the approach described here and the LSC-IVR approach. In this study, the dissipative exciton Hamiltonian16,17 model for the FMO complex is used. We assume initial conditions of the FeynmanVernon type, where the initial density may be decomposed into electronic subsystem and environmental components, F^W(X) = F^sFW,e(X). We consider various initial states for the electronic subsystem where a single excitation is localized on one of the chromophoric sites, i.e., F^s = |λæÆλ|. We also assume that the phonon bath associated with each site, λ, is initially in thermal equilibrium, " # Y βωl β 2 ðλÞ 00 ðPl þ ω2l Rl2 Þ ð6Þ FW , e ðXÞ ¼ 00 exp 2πul ul l
0
λλ N where ~ Fm(x,X) = [1/(2πh h) ]∑λλ0 gλλ0 FW (X), and
gλλ0 ðxÞ ¼
2N þ 1 x2 =p e p p rλ rλ0 þ pλ pλ0 þ iðrλ pλ0 pλ rλ0 Þ δλλ0 ð2Þ 2 λλ0
∑
The subscript W indicates the partial Wigner transform with respect to the environmental degrees of freedom, whose phase space coordinates are denoted by X = (R,P), and the phase space coordinates x = (r,p) and state indices (λ,λ0 ) refer to the quantum ^W(R,P) subsystem. The mapping representation of the operator B may be written as Bm ðx, XÞ ¼
ð5Þ
∑
0 1 B λλ W ðXÞ 2p λλ0
where u00l = (βhhωl/2) coth(βhhωl/2), and FW,e(X) = ΠλF(λ) W,e(X). All simulations reported employ a discretized spectral density of the Debye-type12 with τc = 50 fs, and λR = 35 cm1. Figure 1 displays the time evolution of the initial state population when the excitation is initially located at site 1 at T = 77 K. For comparison, the figure also provides the benchmark results from the reduced hierarchy equation (RHE) formalism,17 which is exact for this particular model of the FMO system, and from the iterative linearized density matrix (ILDM) scheme,20 which is more reliable but computationally much more demanding than PBME, together with approximate LSC-IVR19 and the linearized approach to nonadiabatic dynamics with mapping (LAND-map) results.13 In this case, the
ðrλ rλ0 þ pλ pλ0 þ iðrλ pλ0 pλ rλ0 Þ pδλλ0 Þ ð3Þ ^ where Bλλ W (X) are the matrix elements of BW in the subsystem basis. In this basis, the QCLE is5,6 0
d Bm ðx, X, tÞ ¼ fHm , Bm ðtÞgx, X dt 0 ! p ∂hλλ0 @ ∂ ∂ ∂ ∂ ∂ Bm ðtÞ þ þ 8 λλ0 ∂R ∂rλ0 ∂rλ ∂pλ0 ∂pλ 3 ∂P
∑
ð4Þ 809
dx.doi.org/10.1021/jz200059t |J. Phys. Chem. Lett. 2011, 2, 808–812
The Journal of Physical Chemistry Letters
LETTER
Figure 1. Initial state population as a function of time for the FMO model with Debye bath at T = 77 K. The initial state corresponds to the exciton localized at site 1.
Figure 2. Selected concurrences, Cij = 2|Fij|, as a function of time for the FMO model with Debye bath at T = 77 K. The initial state corresponds to the exciton localized at site 1.
PBME does not provide a very accurate estimate of the initial state population as given by the RHE and ILDM results, but performs approximately as well as the LAND-map technique. The LSC-IVR results indeed provide a better estimate of the initial state population than both the PBME and LAND-map results. However, the PBME captures the time scale and size of the oscillations in the initial state population extremely well. The LAND-map approach also seems to reproduce the correct oscillation frequency, but the size of the oscillations are smaller than in the RHE and ILDM curves. The LSC-IVR result does not correctly capture these features, and becomes out of phase with the other curves after approximately 300 fs. Additionally, the LSC-IVR results for this parameter set, presented in ref 19, showed negative populations for some of the sites in the model. In contrast, PBME simulations for the seven-state FMO model produced nonnegative site populations both at 77 and 300 K (see Figure S1 in Supporting Information). In fact, this different behavior is likely related to the characteristics of the PBME formulation. As explained previously, due to the ignored term from eq 4, part of the dissipative back reaction from the quantum subsystem to the bath is ignored.6 This artifact will exaggerate the effect of the bath-tosubsystem fluctuation, which will eventually drive the elements of the quantum subsystem in the long time limit to populate according to the equilibrium condition at infinite temperature. Thus, care must be taken in analyzing the data regarding population evolution with the PBME method. Nonetheless, the exhibition of correct oscillatory behavior from PBME simulation is encouraging, especially for its future application toward the understanding of coherence with realistic simulations, where the description of short time (up to ∼1 ps) oscillations is crucial. In addition to the time evolution of the diagonal elements of the subsystem density matrix shown above, the time evolution of selected off-diagonal elements is also presented. In Figure 2 the concurrence Cij, which combines both the real and imaginary parts of the density elements in the subsystem basis, is shown. Although quantum coherence and concurrece are, in general, distinctly different quantities, the two are intimately related as the coherence in the site basis is necessary and sufficient to produce concurrence.21 Since the adopted photosynthetic model is restricted to the manifold of single-excitations, the concurrence, which is typically used as a measure of bipartite quantum
entanglement,21,22 can be expressed as Cij = 2|Fij| in the site-basis. Even though a rigorous discussion of quantum coherence may only be attained by quantitatively considering the entire density matrix, in the adopted photosynthetic model the concurrence represents the presence of sitesite superposition states, which are at least qualitatively related to the coherent behavior of the subsystem.21 At the very least, monitoring the behavior of off-diagonal elements with PBME formalism in relation to other methods will provide additional information regarding the performance of PBME. As is shown in Figure 2, the PBME results, at T = 77 K, seem to match the benchmark ILDM results with excellent agreement at short times, up to approximately 100 fs. At longer times the shape of the ILDM curves are well reproduced by the PBME results, but differences in the magnitudes are rather large. This trend is also in accord with the behavior of the population growth shown in Figure 1. At physiological temperature, T = 300 K, the quality of the PBME results increases significantly, and very good agreement with the benchmark RHE and ILDM results is achieved. In fact, Figure 3 shows that the PBME results for various state populations are essentially equivalent in quality to the ILDM results, and reproduces RHE results very closely. However, the PBME simulation algorithm is much less computationally intensive than that of the ILDM, and will require orders of magnitude fewer trajectories to obtain converged results. For example, the PBME results shown in Figure 1 required an ensemble size of N = 5 105 trajectories, whereas corresponding ILDM results from ref 20 required N = 4 108. The ILDM has previously been shown20 to yield excellent agreement with the RHE in the treatment of the electronic concurrences. This is also the case for the PBME, as shown in Figure 4 where various sitesite concurrences are compared with the corresponding ILDM results. Again, excellent agreement is obtained while requiring less computational effort. In conclusion, we have presented the results of quantumclassical dynamics simulations for the seven-state FMO model in the context of the PBME approximation and have provided comparisons to the results of other recent nonperturbative, nonMarkovian techniques. This approach can indeed yield reasonable and trustworthy physical results for this important model system, particularly in the subpicosecond regime at physiological temperatures. This will be important for future work on more realistic model systems, as the utility of more reliable approaches 810
dx.doi.org/10.1021/jz200059t |J. Phys. Chem. Lett. 2011, 2, 808–812
The Journal of Physical Chemistry Letters
LETTER
’ ACKNOWLEDGMENT We are grateful to the research groups of Profs. David Coker, Graham Fleming, and William H. Miller for providing their data for comparison. We are also grateful to Prof. Raymond Kapral for his helpful comments on an earlier draft of this work. This work was financially supported by the Basic Science Research Program (Grant No. 2009-0067085) through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science, and Technology. Supercomputer time from the Korea Institute of Science and Technology Information (KISTI) under Grant No. KSC-2010-C1-0042 is also thankfully acknowledged. ’ REFERENCES (1) 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. (2) Lee, H.; Cheng, Y.-C.; Fleming, G. R. Coherence Dynamics in Photosynthesis: Protein Protection of Electronic Coherence. Science 2007, 316, 1462–1465. (3) Collini, E.; Scholes, G. D. Coherent Intrachain Energy Migration in a Conjugated Polymer at Room Temperature. Science 2009, 323, 369–373. (4) Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Dyanmics. J. Chem. Phys. 1999, 110, 8919–8929. (5) Kim, H.; Nassimi, A.; Kapral, R. Quantum-Classical Liouville Dynamics in the Mapping Basis. J. Chem. Phys. 2008, 129, 084102/ 1–084102/6. (6) Nassimi, A.; Bonella, S.; Kapral, R. Analysis of the QuantumClassical Liouville Equation in the Mapping Basis. J. Chem. Phys. 2010, 133, 134115/1–134115/11. (7) Kapral, R. Progress in the Theory of Mixed Quantum-Classical Dynamics. Annu. Rev. Phys. Chem. 2006, 57, 129–157. (8) Grunwald, R.; Kelly, A.; Kapral, R. Quantum Dynamics in Almost Classical Environments. In Energy Transfer Dynamics in Biomaterial Systems; Burghardt, I., May, V., Micha, D., Bittner, E., Eds.; Springer: Berlin, 2009; pp 383413. (9) 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. (10) Stock, G.; Thoss, M. Semiclassical Description of Nonadiabatic Quantum Dynamics. Phys. Rev. Lett. 1997, 78, 578–581. (11) Sun, X.; Wang, H.; Miller, W. H. Semiclassical Theory of Electronically Nonadiabatic Dynamics: Results of a Linearized Approximation to the Initial Value Representation. J. Chem. Phys. 1998, 109, 7064–7074. (12) Wang, H.; Song, X.; Chandler, D.; Miller, W. H. Semiclassical Study of Electronically Nonadiabatic Dynamics in the Condensed Phase: Spin-Boson Problem with Debye Spectral Density. J. Chem. Phys. 1999, 110, 4828–4840. (13) Bonella, S.; Coker, D. F. LAND-map, A Linearized Approach to Nonadiabatic Dynamics Using the Mapping Formalism. J. Chem. Phys. 2005, 122, 194102/1–194102/13. (14) Dunkel, E. R.; Bonella, S.; Coker, D. F. Iterative Linearized Approach to Nonadiabatic Dynamics. J. Chem. Phys. 2008, 129, 114106/ 1–114106/15. (15) Fenna, R. E.; Matthews, B. W. Chlorophyll Arrangement in a Bactriochlorophyll Protein from Chlorobium Limicola. Nature 1975, 258, 573–577. (16) Adolphs, J.; Renger, T. How Proteins Trigger Excitation Energy Transfer in the FMO Complex of Green Sulfur Bacteria. Biophys. J. 2004, 91, 2778–2797. (17) Ishizaki, A.; Fleming, G. R. Theoretical Examination of Quantum Coherence in a Photosynthetic System at Physiological Temperature. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 17255–17260.
Figure 3. Exciton site populations as a function of time for the FMO model with Debye bath at T = 300 K. The initial state corresponds to the exciton localized at site 6.
Figure 4. Concurrences, Cij = 2|Fij|, as a function of time for the FMO model with Debye bath at T = 300 K: (a) C45, (b) C56, (c) C46, (d) C67. The initial state corresponds to the exciton localized at site 6. Full lines correspond to the benchmark ILDM results, and dashed lines are from PBME simulations.
such as RHE are currently limited to the adopted simplified photosynthetic model. These results thus provide validation for possibility of incorporating the PBME into larger (potentially allatom) MD simulations of photosynthetic and other light-harvesting reaction systems under physiological conditions, where the computational efficiency of the PBME formalism would prove to be particularly desirable.
’ ASSOCIATED CONTENT
bS
Supporting Information. Population development on seven sites from PBME simulations. This material is available free of charge via the Internet at http://pubs.acs.org.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected] (A.K.);
[email protected] (Y.M.R.). 811
dx.doi.org/10.1021/jz200059t |J. Phys. Chem. Lett. 2011, 2, 808–812
The Journal of Physical Chemistry Letters
LETTER
(18) 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/1–234111/10. (19) 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. (20) Huo, P.; Coker, D. F. Iterative Linearized Density Matrix Propagation for Modelling Coherent Excitation Energy Transfer in Photosynthetic Light Harvesting. J. Chem. Phys. 2010, 133, 184108/ 1–184108/12. (21) Sarovar, M.; Ishizaki, A; Fleming, G. R.; Whaley, K. B. Quantum Entanglement in Photosynthetic Light-Harvesting Complexes. Nat. Phys. 2010, 6, 462–467. (22) Hill, S.; Wootters, W. K. Entanglement of a Pair of Quantum Bits. Phys. Rev. Lett. 1997, 78, 5022–5025.
812
dx.doi.org/10.1021/jz200059t |J. Phys. Chem. Lett. 2011, 2, 808–812
