Energy transfer properties and absorption spectra of the FMO complex: from exact PIMC calculations to TCL master equations

0.79 MB·English

by  Piet Schijven

Energy transfer properties and absorption spectra of
the FMO complex: from exact PIMC calculations to
3
1
0
2
TCL master equations
n
a
J

9
Piet Schijven,∗ Lothar Mühlbacher,∗ and Oliver Mülken∗
]
h
p
- PhysikalischesInstitut,UniversitätFreiburg,Hermann-Herder-Strasse3,79104Freiburg,
m
e
Germany
h
c
.
s
c E-mail: [email protected]; [email protected];
i
s
y [email protected]
h
p
[

2
v
9
3
8
0
.
1
0
3
1
:
v
i
X
r
a
∗Towhomcorrespondenceshouldbeaddressed
1
Abstract
Weinvestigate theexcitonic energy transfer(EET)intheFenna-Matthews-Olsen complex
andobtainthelinearabsorptionspectrum(at300K)byaphenomenologicaltime-convolutionless
(TCL)masterequation whichisvalidated byutilizing PathIntegralMonteCarlo(PIMC)sim-
ulations. ByapplyingMarcus’theoryforchoosingtheproperLindbladoperatorsforthelong-
time incoherent hopping process and using local non-Markovian dephasing rates, our model
shows very good agreement with the PIMC results for EET. It also correctly reproduces the
linearabsorption spectrum thatisfoundinexperiment, withoutusinganyfittingparameters.
Introduction
Since the seminal experiments on the Fenna-Matthews-Olsen (FMO) complex performed by the
Fleminggroupin20071 andsubsequentlyconfirmed bytheEngelgroupin2010,2 manytheoreti-
calmodelshavebeenproposedinordertoproperlydescribetheobservedlong-lastingoscillations
in the 2D spectroscopic data. While many methods are based on quantum master equations of
some sort or the other, like Lindblad or Redfield equations,3–8 only recently a more detailed de-
scriptionofdissipativeeffectshasbeenattemptedintheformofPathIntegralMonteCarlo(PIMC)
calculations based on results from atomisticmodelingcombiningmolecular-dynamics (MD) sim-
ulationswithelectronicstructurecalculations.9 Incontrasttoaphenomenologicalmodelingofthe
“environment”(proteinscaffold,watermolecules,etc.), thelatterapproach utilizesBChl-resolved
spectral densities which can be directly incorporated in the PIMC calculations to obtain an exact
account ofthefull excitondynamics.
Most quantum master equation approaches use special analytical forms of spectral densities,
such as Ohmic or Lorentzian forms. Although this allows to obtain solutions for the equations, it
remainstobeaphenomenologicalansatz. Sofartheresultsthathavebeenobtainedbyhierarchical
time-nonlocalmasterequationswithaLorentzianspectraldensityshowreasonableagreementwith
theexperimentalresultsofboththeabsorptionand 2D spectraat 77K.10,11
In contrast, numerically exact methods like PIMC simulations12 or the QUAPI method13–15
2
have proven to be capable to produce the exact quantum dynamics of excitonic energy transfer
overtheexperimentallyrelevant timescales. However,this comes at the priceof ratherlarge com-
putational costs. Here, we combine the respective strengths of numerically exact and approxima-
tive methods while overcoming their respective weaknesses: we use PIMC results for the exact
quantum dynamics of the excitonic population dynamics over intermediate timescales (i.e. 600fs)
whichstillallowforafastproductionoftherespectiveresults,yetaresufficienttoallowforacom-
parisontoatime-convolutionless(TCL)masterequationbasedonaLindbladapproach;wefurther
corroborateourresultsbycomparisonwithPIMCdataforupto1.5ps. Foramodeldimersystems
theauthorshavealreadysuccessfullydemonstratedsuchaconcept.16 Here,ouransatzcapturesthe
long-timebehaviorbyrelatingtoMarcus’theoryofelectrontransport17,18 becauseweassumethat
thelong-timebehaviorisgovernedbyaclassicalhoppingprocessbetweentheindividualsites,18,19
seebelow. Fortheshort-timedephasingbehaviorweintroducenon-Markoviandephasingrates.20
Since, in principle, we then obtain results for arbitrarily long times, we have an efficient yet very
accurate way to calculate arbitrary transfer properties. Furthermore, we use the master equation
to calculate the absorption spectrum of the FMO complex (at 300 K), an observable which can
straightforwardly be obtained experimentally.21 This opens the possibility to estimate the valid-
ity of the underlying microscopic Hamiltonian as well as its respective parametrization based on
recent results from mixed quantum-classical simulations22 by getting into direct contact with ex-
perimentaldata.
Energy transfer on the FMO complex
Microscopic description The dynamics of single excitations on the FMO complex is often de-
scribed by a tight-binding Hamiltonian with 7 localized sites, corresponding to the 7 bacteri-
ochlorophylls(BChls)oftheFMOmonomer.23,24Theinfluenceoftheproteinscaffoldandsolvent
on the excitonicdynamics is treated, in the spirit of the Caldeira-Leggett model,25 as a collection
ofharmonicmodesthatare linearlycoupled toeach BChl. Previousstudiesshowedno significant
3
correlationsbetweenthebathinducedenergyfluctuationsatdifferentsites,26,27 soweassumethat
each BChl is coupled to its own individual environment. The full Hamiltonian of the system can
nowbewrittenas:
H =H +H +H , (1)
S B SB
with
H = (cid:229) e |nihn|+ (cid:229) J |mihn|, (2)
S n mn
n m6=n
HB = n(cid:229) ,k (cid:18)2Pmn2nkk +21mnk w n2k Xn2k (cid:19), (3)
HSB = (cid:229) |nihn| cnk Xnk +L (ncl) . (4)
n
(cid:16) (cid:17)
Thestate|nicorrespondstothesingle-excitationstateofsiten,theparametere denotestheenergy
n
gap between ground and excited state of site n, and J describes the excitonic coupling between
mn
sites m and n. Furthermore, Xnk , Pnk , mnk and w nk denote the position, momentum, mass and
frequency of the bath oscillators, respectively. In the interaction Hamiltonian H , the constants
SB
cnk (in units of eV/m) denote the coupling strength between site n and the bath modes. We have
included the classical reorganization energies L (cl) as a counter-term in H to prevent further
n SB
renormalizationofthesiteenergies bytheenvironment.20,25,28 Thisquantityisdefined as:
h¯ ¥ J (w )
L (cl) = dw n , (5)
n p w
0
Z
whereJ (w )(inunitsof1/s)isthespectraldensityofthebaththatiscoupledtositen. Intermsof
n
thesystemparameters, itis givenby:
Jn(w )= ph¯ (cid:229) k 2mcnk2nkw nk d (w −w nk ). (6)
The precise numerical values of the different parameters entering in the expressions above were
obtained from combined quantum-classical simulations for the full FMO complex including the
4
solvent.22,29
Effective master equation approach
We use now use the microscopic description of the FMO complex to set up a phenomenological
second order time-local quantum master equation. In doing so, we are able to reproduce the dy-
namics obtained from the PIMC simulationsas well as extendingit to, in principle, arbitrary long
times. Additionally, our approach also allows to obtain results for the linear absorption spectrum
whichare incloseaccordance toexperimentalfindings.
The spectral density of the FMO complex22 leads to reorganization energies L (cl) of the order
n
of 0.02−0.09 eV, which is comparable to the differences in the site energies e , while the exci-
n
tonic couplings J are of the order of 1 meV. This implies that we can expect that the protein
mn
environment is relatively strongly coupled to the FMO complex and that it therefore leads to a
strong damping for the population dynamics. This is also reflected by the results of the PIMC
simulations.9
We now assume that in the long-time limit, after most of the coherences (i.e. off-diagonal el-
ements of the reduced density matrix in the site-basis representation) in the system have decayed,
EET can be described by a classical hopping process between the different sites (BChl’s), that is
inducedby theproteinenvironment. Thetransfer rates k havetosatisfydetailed balance, ensur-
mn
ing acorrect equilibriumstate, and are assumed to followfrom Fermi’s golden rule. Furthermore,
the rates should also depend on the reorganization energies L (cl) and L (cl) of the baths that are
n m
coupled to the sites n and m, reflecting the differences in the coupling strengths of the protein en-
vironment to each BChl. Unlike Förster theory, which assumes incoherent hopping between the
energy eigenstatesofH ,30 Marcus’stheoryofelectrontransportsatisfiesalltheseproperties,17,18
s
leadingto transferrates k oftheform:
mn
pb b (e −e +L (cl))2
k = |J |2exp − n m mn , (7)
mn mn
sh¯2L (cl) " 4L (cl) #
mn mn
5
withb =1/k T and L (cl) =L (cl)+L (cl).
B mn m n
Asidefromincoherenttransferbetweenthesites,theenvironmentalsoinducesastrongdephas-
ing on each site. In the framework of the second order TCL master equation,20 these dephasing
rates (in unitsof1/fs)aregivenby:
t ¥
l (t)=2Re ds dw J (w )[coth(b h¯w /2)cos(w s)−isin(w s)]. (8)
n n
0 0
Z Z
Here, we use the spectral densities J (w ) which have been obtained by MD simulations in Ref.22
n
and numericallycalculatethecorrelation function.
TheTCL masterequationthatdescribes theexcitationdynamicscan nowbewrittenas:20
dr (t) i
≡L(t)r (t)=− [H ,r (t)]+D(t)r (t). (9)
s
dt h¯
Our numerical results (not displayed) show that the Lamb shift term that usually appears in this
equation, only leads to a negligible difference in both the population dynamics and the linear
absorptionspectrum (thepositionof thepeak is shifted by approximately-1 meV). The dissipator
D(t)isassumedtotakethefollowingLindbladform, according totheconsiderationsabove:20
D(t)r (t)=(cid:229) g (t) L r L† −1 L† L ,r . (10)
mn mn mn 2 mn mn
mn (cid:18) (cid:19)
n o
TheLindbladoperatorsare defined by L =|mihn|andtherates by g (t)=l (t)and g (t)=
mn mm m mn
k for m 6= n. The operators L model the dephasing process, while the operators L model
mn mm mn
the incoherent transfer between sites m and n. This choice of Lindblad operators will lead - in the
long-time limit - to incoherent hopping transfer between the sites, which is different from, e.g.,
Redfield theory, which requires incoherent transfer between the eigenstates |y i of H , leading to
S
LindbladoperatorsoftheformL∼|y ihy |.18,20
n m
However, we note that the equilibrium state of our master equation (r ) is slightly differ-
eq
ent from the one that follows from Marcus theory r , which is given by detailed balance,
eq,db
6
limt→¥ r nn(t)=(1/Z)exp(−be n) and limt→¥ r mn(t)=0 for m 6=n, where r mn(t)= hm|r (t)|ni.
This can be shown by noting that D(t)r = 0 but [H ,r ] 6= 0. This implies that r is
eq,db S eq,db eq,db
not astationary stateofour masterequation. For thepresent calculation, thedevationsare only of
theorderof1%,so westillexpect ourapproach to givegoodresults.
Path Integral Monte Carlo simulations
PIMCsimulationsallowtoextracttheexactquantumdynamicsinthepresenceofadissipativeen-
vironment,bothforchargetransport12,31,32 aswellasenergytransfer.9 Inshort,thetimeevolution
of the reduced density operator of a dissipative quantum systems is calculated by employing the
pathintegralrepresentationforthepropagatoraccording totheFeynman-Vernontheory33 forfac-
torizingoritsextensiontocorrelatedinitialpreparations.34 Thesepathintegralsarethenevaluated
by a stochastic sampling process based on Markov walks through the configuration space of all
conceivable quantum paths which self-consistently emphasize the physically most relevant ones.
WhilethereisnolimitationwithrespecttothechoiceofsystemparametersforwhichPIMCsimu-
lationare capableof producingnumerically exact results, thisapproach is subjectto the notorious
‘dynamical sign problem’,35 which reflects quantum-mechanical interferences between different
systempathsand resultsinan increase ofthecomputationaleffortnecessary toobtainstatistically
converged results which scales exponentially with the timescale over which the dynamics of the
system is investigated. However, the presence of a dissipativeenvironment substantially weakens
these interference effects and therefore the sign problem. Furthermore, it allows for various effi-
cientoptimizationschemeswhichleadtoafurthersoothingofthiscomputationalbottleneck,thus
significantlyenlargingtheaccessibletimescales.12,36,37
Forthepresentcase, weutilizethePIMCdatapresentedinRef.9 todemonstratethereliability
ofthemasterequationresultsandextendsomeoftheformertolongertimescales. Tothatextend,a
factorizinginitialpreparationhas beenemployed,where, resemblingthesituationpriorto thecre-
ationofan exciton,thebathmodesinitiallyarein thermalequilibriumwithrespect tothemselves,
while the exciton has been modeled to be either initially localized on one of the seven BChl sites
7
orinoneofthesevenexcitoniceigenstates.
Populationdynamics
In Fig. 1 we show the population dynamics that is obtained by solving the TCL master equation,
Eq. (9), forinitialconditionscorrespondingtoalocalizationonthesites|ni,e.g. r (0)=|nihn|,in
comparisontothecorrespondingnumericallyexactPIMCresults. Thedottedcurvesrepresentsthe
latterand the solidlines represent theresults from ourmaster equationapproach. Fig. 2 showsan
extensionoftheresultsupto1200and1500fsforinitialpreparationsinsites1and6,respectively.
In general we observe good quantitative agreement of our approach with the PIMC results
for all localized initial preparations and over all observed timescales. The largest deviations are
observedforaninitialconditionlocalizedonsite4forwhichthebathhasthelowestreorganization
energy (0.025 eV). Since our assumption of a classical hopping process at long-times is based on
having a strong coupling to the environment, we would expect that our approach becomes worse
with decreasing reorganization energy. Also, from Fig. 2 one observes a good agreement in the
approach toequilibrium,althoughthedecay isslightlyslowerthanpredictedby thePIMCresults.
Figure 3 corroborates our results. Here, the excitonic excitation is initially in one of the seven
eigenstates |y i of H . Again we find very good agreement with the PIMC results, where once
n S
more the strongest deviations occur for the initial preparation exhibiting the largest population
on site 4. We note that there is no fitting parameter involved. Introducing a parameter which
interpolates between purely coherent and purely incoherent transfer, as in,38–40 could lead to a
furtherimprovementoftheagreement. Nevertheless,alreadythisrathersimplephenomenological
model captures most of the details which are present in the PIMC calculations. Additionally, it
allowsforacomputationallycheap calculation ofthelinearabsorptionspectrum.
8
1 1
r (0)=|1æÆ 1| r (0)=|2æÆ 2|
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
r (0)=|3æÆ 3| r (0)=|4æÆ 4|
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
r (0)=|5æÆ 5| r (0)=|6æÆ 6|
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
r (0)=|7æÆ 7|0 100 200 300 400 500 600
0.8 t (fs)
0.6
Site 1
0.4 Site 2
Site 3
Site 4
0.2
Site 5
Site 6
Site 7
0
0 100 200 300 400 500 600
t (fs)
Figure1: Comparisonofpopulationdynamicsofthe7differentsitesoftheFMOcomplexobtained
from the numerically exact PIMC results (circles) with the results from the TCL master equation
approach (solid lines) for different localized initial conditions on the sites |ni, n =1,...,7. Note
that the statisticalerror of thePIMC calculationsis typicallysmallerthan thesymbolsize. There-
forewedo notshowtheerrorbars.
9
1 1
Site 1
r (0)=|1æÆ 1| r (0)=|6æÆ 6| Site 2
0.8 0.8 Site 3
Site 4
Site 5
0.6 0.6 Site 6
Site 7
0.4 0.4
0.2 0.2
0 0
0 200 400 600 800 1000 0 200 400 600 800 1000 1200
Figure2: SameasFig.1forlocalizedinitialpreparationsinsites1and6,butnowextendedbeyond
1ps
.
Table 1: The numerical values for the x-, y-, and z- component as well as the absolute value
squared of the transition dipole moments~m =(m ,m ,m ) in units of Debye [D]. The z
m m,x m,y m,z
axisischosenalongtheC -symmetryaxisoftheFMOcomplex,andtheyaxisischosentobe
3
parallelto the N −N axisofBChl 1.
B D
m 1 2 3 4 5 6 7
m 0.0 -6.10 -5.27 0.0 -6.39 5.16 0.0
m,x
m 1.86 1.08 -3.04 2.49 0.0 2.98 -1.14
m,y
m 6.07 1.66 -2.10 5.85 -0.45 2.29 5.85
m,z
|~m |2 40.32 41.09 41.47 40.45 41.09 40.70 35.52
m
10