The intermolecular contribution to the spectral density of the exciton-vibrational coupling of the homotrimeric Fenna–Matthews–Olson (FMO) light-harvesting protein of green sulfur bacteria *P. aestuarii* is analyzed by combining a normal mode analysis of the protein with the charge density coupling method for the calculation of local transition energies of the pigments. Correlations in site energy fluctuations across the whole FMO trimer are found at low vibrational frequencies. Including, additionally, the high-frequency intrapigment part of the spectral density, extracted from line-narrowing spectra, we study intra- and intermonomer exciton transfer. Whereas the intrapigment part of the spectral density is important for fast intramonomer exciton relaxation, the intermolecular contributions (due to pigment-environment coupling) determine the intermonomer exciton transfer. Neither the variations of the local Huang–Rhys factors nor the correlations in site energy fluctuations have a critical influence on energy transfer. At room temperature, the intermonomer transfer in the FMO protein occurs on a 10 ps time scale, whereas intramonomer exciton equilibration is roughly two orders of magnitude faster. At cryogenic temperatures, intermonomer transfer limits the lifetimes of the lowest exciton band. The lifetimes are found to increase between 20 ps in the center of this band up to 100 ps toward lower energies, which is in very good agreement with the estimates from hole burning data. Interestingly, exciton delocalization in the FMO monomers is found to slow down intermonomer energy transfer, at both physiological and cryogenic temperatures.

## I. INTRODUCTION

The Fenna–Matthews–Olson (FMO) protein mediates energy transfer between an outer antenna system in green sulfur bacteria, termed a chlorosome, and the reaction center complex. Since the discovery of its molecular structure in 1975,^{1} it has served as an important model system for the development of structure-based theories and calculation schemes for excitation energy transfer and optical spectra of light-harvesting complexes.^{2–17} The FMO protein is one of the first systems to which the newly developed 2D electronic spectroscopy was applied,^{18} revealing long-lived coherent oscillations^{19–23} that founded the field of quantum biology^{24,25} and brought this protein into the attention of a general audience from many different fields. It is becoming more and more clear that the origin of the long-lived part of the oscillations is due to vibrational rather than electronic coherences.^{21–23,26,27} It is not the protection of inter-exciton state coherences, but the equal strength of nearest neighbor pigment–pigment and local pigment–protein coupling that holds the key for efficient excitation energy transfer.^{28,29} In this way, the excitons become partially localized over just a few pigments, exciton relaxation corresponds to a spatially directed energy transfer toward the low-energy pigments, and the excess energy of excitons can be well dissipated by the vibrations.^{3,4,25} A protection of inter-exciton state coherences would be detrimental for exciton relaxation.^{9}

The FMO protein is a homotrimer, where each monomer binds eight bacteriochlorophyll *a* (BChl *a*) pigments (Fig. 1). Seven of the BChls are bound inside a protein bag made of mostly *β*-sheets, and the eighth is bound at the surface. During biochemical isolation of the FMO protein, the eighth BChl is easily lost. Therefore, in the original crystal structure, only the seven central BChls are present, and the eighth was found later.^{30,32} The eighth BChl is bound between two FMO monomers, and its central magnesium atom has a ligand from both monomers. The assignment of the three respective BChl 8 pigments in the crystal structure^{30} is such that each of them belongs to the protein monomer with the closest atom-to-atom distance between pigment and protein. This assignment has caused some confusion in the literature since in terms of pigment–pigment distances, the assignment would be different. Grouping the pigments according to the strength of their excitonic couplings results in BChls 1–7 and 8″ in the first domain, BChls 8 and 1′–7′ in the second domain, and BChls 8′ and 1″–7″ in the third domain, where a prime (double prime) denotes assignment to monomer 2 (3) in the crystal structure.

The inhomogeneous charge distribution of the protein generates a site energy funnel in the FMO monomer, where those pigments that are closer to the reaction center complex are more redshifted than those facing the baseplate and the outer chlorosome.^{3–5,7,33,34} Therefore, the low-energy exciton states are localized close to the reaction center complex, and directed energy transfer becomes possible. In addition, the inhomogeneous charge distribution of the protein also leads to different coupling constants between the vibrational normal modes of the complex and the electronic excitations of the pigments (site energies).^{9} In this way, correlations between site energy fluctuations are suppressed, and the excess energy of excitons can be dissipated efficiently. In our earlier work^{9} on the spectral density of the monomeric subunit of the FMO protein, we found correlations in site energy fluctuations, but only at vibrational frequencies, which are too small to affect exciton relaxation.

One of the aims of the present study is to investigate whether these correlations in site energy fluctuations extend over the whole FMO trimer and, if so, how they influence the energy transfer between different FMO monomers. Due to the C3 symmetry of the complex, we expect similar local exciton-vibrational coupling constants of equivalent sites and similar exciton energies in the three monomers. Hence, in principle, the former could lead to correlated site energy fluctuations for the delocalized normal modes, and the respective low vibrational frequencies could match energy differences between exciton states in different monomers. In order to investigate this question, we will extend generalized Förster theory to include correlated site energy fluctuations. The intermolecular (low-frequency) part of the spectral density containing these correlations is obtained by combining a normal mode analysis (NMA) of the whole pigment–protein complex with the charge density coupling (CDC) method for the local optical transition energies of the pigments, as described earlier.^{9} Since in this earlier analysis of FMO monomers, it was found that the fluctuations of excitonic couplings are two orders of magnitude smaller than those of the site energies, we do not consider the former in the present work. The high-frequency part of the spectral density that is due to intrapigment vibrations will be assumed site-independent and taken from an analysis of delta-fluorescence line narrowing (ΔFLN) spectra of the FMO protein.^{35} A comparison of our approach to calculate the spectral density with more advanced quantum mechanics/molecular mechanics (QM/MM) methods from the literature^{8,10,15,16,36,37} including specific results on the FMO protein will be given in Sec. V A.

Extensive experimental information exists about intramonomer exciton relaxation and intermonomer energy transfer in the FMO protein from time-resolved pump–probe^{38–43} and 2D electronic spectroscopy^{44} and from frequency-resolved hole burning spectra.^{45–47} Intramonomer exciton relaxation occurs mostly on a sub-picosecond timescale, where the lifetimes of high-energy exciton states are below 100 fs, even at cryogenic temperatures.^{41,44} The intermonomer transfer is orders of magnitude slower and can be detected best at cryogenic temperatures, where the lifetime of the lowest exciton state in two of the three monomers in the FMO protein is determined by the energy transfer to the monomer with the lowest energy exciton state of the trimer. Note that due to static disorder in site energies, for every realization of disorder, the three lowest exciton states are likely to be located in different monomers of FMO having slightly different energies. At higher temperatures, the uphill energy transfer within the FMO monomers most likely out-competes the intermonomer downhill transfer, making the latter hard to detect. To the best of our knowledge, there is only one room temperature study that inferred an intermonomer transfer time in the FMO protein, suggesting a time constant of 1.4 ps–2 ps.^{38}

Intermonomer transfer times in the 20 ps–40 ps range were estimated from low-temperature transient absorption and 2D electronic photon echo experiments.^{42–44} From the homogeneous broadening of the low-energy exciton band, measured in hole burning spectroscopy on FMO complexes of *P. aestuarii* at 4.2 K, a lifetime of 26 ps was obtained at a wavelength of 822.8 nm and one of 99 ps at the longer wavelength of 824.8 nm.^{47} Similar values were determined for the FMO protein from *Chlorobaculum* (formerly *Chlorobium*) *tepidum*.^{46}

Intermonomer excitation energy transfer at cryogenic temperatures was studied theoretically by Jankowiak and co-workers^{48–50} using two models. In the first model, excitation energy transfer between the lowest sites of each monomer was investigated with standard Förster theory.^{48,50} Averaged intermonomer transfer times of 41 ps and 70 ps were obtained. In the second model, exciton delocalization between the pigments in the FMO monomers was taken into account, and the transfer was described by using generalized Förster theory. Averaged times of 17 ps, 16 ps, and 47 ps were reported for downhill transfer between the low-energy exciton states of the different monomers.^{49,50} These results seem to suggest that the intramonomer delocalization of the low-energy exciton states promotes intermonomer exciton transfer. A direct comparison of the transfer times in the two models is, however, difficult since different spectral densities were used in these calculations.

In the present work, besides re-investigating the above effect, we extend these calculations in the following ways: (i) Possible correlations in site energy fluctuations are studied by extending generalized Förster theory and using an intermolecular spectral density, calculated by combining an NMA with the charge density coupling method.^{5,9} (ii) Frequency-resolved time constants of intermonomer transfer are calculated and compared with the results from hole burning spectroscopy.^{47} (iii) In addition to cryogenic temperatures, we investigate intramonomer exciton relaxation and intermonomer transfer at physiological temperatures.

The remainder of this work is organized in the following way. We start with a summary of the Frenkel exciton Hamiltonian and Redfield theory of exciton relaxation. An expression is derived for the generalized Förster rate constant that includes correlations in site energy fluctuations. Certain limiting cases such as a classical description of nuclear motion giving rise to a simple analytical expression for the rate constant, the standard generalized Förster theory, and Förster theory are discussed as well. Afterward, we summarize the calculation of the spectral density and the computational methods. The theories and calculation schemes are finally applied to the FMO protein with particular focus on the intermonomer transfer.

## II. THEORY

For the pigments within a monomeric subunit of the FMO protein, the local reorganization energies of the exciton-vibrational coupling are in the same order of magnitude as their nearest neighbor excitonic couplings. Therefore, in the basis of localized excited states, there is no small parameter that can be treated perturbatively. Numerically exact methods such as the hierarchical equation of motion (HEOM) approach^{51,52} exist, which are, however, computationally expensive, in particular, when arbitrary functional forms of the spectral density and low temperatures are investigated. A numerically efficient alternative is to transform the exciton-vibrational coupling into the basis of exciton states, defined as the eigenstates of the electronic Hamiltonian taken at the equilibrium position of nuclei in the electronic ground state of the complex. For the present system, the off-diagonal elements of the exciton-vibrational coupling in this basis are smaller than the diagonal elements.^{9,25} We take advantage of this difference by treating the diagonal elements exactly and applying a Markov and secular approximation to the off-diagonal elements.^{53} The latter approximations results in a Redfield rate constant for the relaxation between different exciton states. HEOM calculations have found^{54} only small non-secular contributions to the linear absorption spectrum of the FMO monomer, indicating^{55} that the exciton basis is already close to the true eigenbasis that is optically excited.

Due to the large distances between pigments in different FMO monomers (Fig. 1), the excitons dynamically localize inside the FMO monomers. Such a localization of exciton states in the light-harvesting complex II of higher plants was investigated^{56} by comparing the exact HEOM calculations with the combined modified Redfield/Generalized Förster theory.^{57} In the latter approach, dynamic localization is taken into account implicitly by assigning exciton domains of strongly coupled pigments and allowing exciton delocalization only within these domains. A pigment belongs to an exciton domain if it is coupled stronger than a given cut-off value *V*_{c} to at least one other pigment of this domain. In the case of LHC-II, the best agreement between the combined modified Redfield/Generalized Förster theory and the exact HEOM solution was obtained for *V*_{c} = 30 cm^{−1}. In general, *V*_{c} is expected to be in the same order of magnitude as the local reorganization energy of the exciton-vibrational coupling of the pigments. If two pigments are excitonically coupled much weaker than their local reorganization energies, the system can minimize its free energy by localizing the excited states. For the present system, the intermonomer excitonic couplings are at least one order of magnitude smaller than the local reorganization energies of the pigments, and the intramonomer nearest neighbor excitonic couplings are in the same order of magnitude but somewhat larger than these reorganization energies. Therefore, as noted already above, we identify the exciton domains with the monomeric subunits of the FMO protein.

Concerning intradomain (intramonomer) exciton relaxation in the FMO protein, it has been found that modified Redfield theory and Redfield theory give very similar results.^{3} We note in passing that, in general, Redfield theory not necessarily is less accurate than modified Redfield theory, as has been reported in a model study,^{58} where these two theories were compared with the more accurate non-equilibrium modified Redfield theory.^{58,59} For a large parameter range, Redfield theory was found to be even more accurate than modified Redfield theory due to compensation effects between two different approximations. Here, we will describe intramonomer exciton relaxation by Redfield theory and intermonomer energy transfer by generalized Förster theory, where the diagonal elements of the exciton-vibrational coupling in the basis of intramonomer exciton states are treated non-perturbatively. We extend the latter theory by taking into account correlations in site energy fluctuations of pigments in different exciton domains (FMO monomers).

### A. Hamiltonian

The Hamiltonian of the FMO trimer is expressed as

where *H*_{a} contains the Hamiltonian of the exciton domain (FMO monomer) *a* and $V^$ is the excitonic couplings between different domains. We expand the above Hamiltonian in the basis of excited states $Ma$ that are delocalized in the different exciton domains. The exciton state $Ma$ of domain *a* is given as a linear combination of local excited states $ma$ of this domain, $Ma=\u2211macma(Ma)ma$, where the coefficients $cma(Ma)$ are obtained from the eigenvectors of the exciton matrix. This matrix contains in the diagonal the local transition energies $Ema$ (site energies) and in the off-diagonal the excitonic couplings $Vma,na$ between the pigments in this domain. The respective eigenenergy $\u03f5Ma=\u210f\omega Ma$ is the vertical excitation energy of the exciton state $Ma$ at the equilibrium position of nuclei in the electronic ground state.

The Hamiltonian *H*_{a} in Eq. (1) reads

with

where *T*_{nucl} denotes the kinetic energy of nuclei and *Q*_{ξ} is the dimensionless coordinate of normal mode *ξ* with eigenfrequency *ω*_{ξ} that is related to the mass-weighted original normal coordinate *q*_{ξ} by $Q\xi =2\omega \xi /\u210f\u2009q\xi $.^{60} The energy difference $\u210f\omega Ma\u2032$ between the minima of the potential energy surfaces of exciton state $Ma$ and the electronic ground state is given as

with the exciton energy $\u210f\omega Ma$ and the reorganization energy $E\lambda (Ma)=\u2211\xi \u210f\omega \xi g\xi (Ma,Ma)2$. The latter contains the diagonal exciton-vibrational coupling constants *g*_{ξ}(*M*_{a}, *M*_{a}) in the exciton basis. The off-diagonal elements of these couplings are contained in $V^a$ in Eq. (2), reading

The exciton-vibrational coupling constant *g*_{ξ}(*M*_{a}, *N*_{a}) of delocalized states is related to the local coupling constant $g\xi (ma)$ of the modulation of the transition energy of pigment *m*_{a} by

The intermonomer coupling $V^$ in Eq. (1) reads

with

where $Vmanb$ is the excitonic coupling between pigments *m*_{a} in domain *a* and *n*_{b} in domain *b*.

### B. Spectral density

We split the total spectral density $Jmktotal(\omega )$ of the local exciton-vibrational coupling into an intermolecular part obtained from a NMA of the pigment–protein complex^{9} and an intrapigment part that is obtained from difference fluorescence line narrowing experiments,^{35}

The intermolecular part

describes the site energy fluctuations of pigment *m* induced by the intermolecular vibrations for *m* = *k*, and it contains the correlations in site energy fluctuations of pigments *m* and *k* for *m* ≠ *k*. The dimensionless coupling constant $g\xi (m)$ originates from the exciton-vibrational coupling Hamiltonian expanded in the basis of local excited states $m$ of the complex,^{9}

which describes how the local excitation energy of site *m* varies if the mass-weighted normal coordinate *q*_{ξ} is displaced. Note that the sum over *m* runs over all local excited states of the complex and not just over those of a single exciton domain. In this way, we take into account the possibility that the vibrational normal modes of the complex can be delocalized over the whole complex, notwithstanding the fact that the excitons are dynamically localized in certain domains, as discussed above. We will, in fact, find evidence for interdomain delocalization of normal modes. The normal coordinates *q*_{ξ} are related to the Cartesian nuclear coordinates **R**_{j} of the complex by

where $Rj(0)$ denotes the equilibrium position of atom *j* in the electronic ground state of the complex, which is obtained from the crystal structure as described in Sec. III. *M*_{j} is the mass of atom *j*, and $Aj(\xi )$ is part of the eigenvector of normal mode *ξ* with eigenfrequency *ω*_{ξ} obtained from the diagonalization of the Hessian matrix in the NMA. The coupling constants $g\xi (m)$ in Eq. (11) contain the gradients of the local transition energy of pigment *m* with respect to the nuclear coordinates **R**_{j} taken at their equilibrium positions $Rj(0)$. Using our CDC method,^{5} these gradients can be calculated analytically, resulting in the following expression for the coupling constant:^{9}

where the sum over *i* runs over all atoms of the complex, except for those atoms that contribute to the difference charge density of pigment *m*. The respective atomic partial charges *q*_{i} at equilibrium position $Ri(0)$ interact with the difference partial charge $qk(m)(1)\u2212qk(m)(0)$ between the excited and the ground state of the *k*th atom of pigment *m* at equilibrium position $Rk,m(0)$, which are summed over by index *k*. In this way, the change in the difference in Coulomb interaction between the excited and the ground state of pigment *m* and the rest of the complex upon normal mode displacement, giving rise to the fluctuation in the transition energy of the pigment, is determined. The prefactor 1/*ϵ*_{eff} takes into account the screening of the Coulomb interaction due to the electronic polarizability of the protein and solvent environment and uncertainties of the quantum chemical method used to determine the atomic partial charges of the pigments in an effective way. Here, we approximate *ϵ*_{eff} by an average optical dielectric constant *ϵ*_{eff} = *ϵ*_{opt} ≈ 2 that can be estimated, e.g., by evaluating the difference in oscillator strength of pigments between the protein and a solvent environment.^{61} The coupling constants $g\xi (m)$ and $g\xi (k)$ finally enter the intermolecular spectral density *J*_{mk}(*ω*) in Eq. (10) that will be used to calculate rate constants of energy transfer below. Note that if a normal mode *ξ* is delocalized, the spectral density *J*_{mk}(*ω*), in general, will have non-zero off-diagonal elements, reflecting a correlation (*J*_{mk}(*ω*) > 0) or anticorrelation (*J*_{mk}(*ω*) < 0) in site energy fluctuations, that will be analyzed below.

The intrapigment part

describes the modulation of the transition energy of the pigments by their intramolecular vibrations. We assume the same local coupling constants $g\xi intra$ for all pigments. The numerical values for $g\xi intra2$, which equal the Huang–Rhys factors of the different modes *ξ*, and the respective vibrational frequencies *ω*_{ξ} were taken from an analysis of ΔFLN spectra of the FMO protein.^{35} We included the first 18 intrapigment modes with frequencies <500 cm^{−1} from Table I of Ref. 35 in our calculations. The vibrational frequencies of these modes range between *ω*_{1} = 46 cm^{−1} and *ω*_{18} = 481 cm^{−1}. The first eight modes exhibit the largest Huang–Rhys factors ranging between $g3intra2=0.009$ and $g8intra2=0.012$. A graphical representation is given in Fig. 2 (green bars).

Note that there are two different definitions of the spectral density *J*(*ω*) in the literature. The present one is directly related to the spectral shape of the vibrational sideband seen in linear absorption and fluorescence.^{53,60} The alternative definition^{10,36} is obtained from the present one by multiplication with the factor *ω*^{2}*πℏ*.

### C. Energy transfer

#### 1. Intradomain exciton relaxation: Redfield theory

Exciton relaxation within the domains is described by Redfield theory, where the rate constant $kMa\u2192NaRedf$ for relaxation between exciton states $Ma$ and $Na$ is given as

with the Bose–Einstein distribution function

the spectral density $Jmanatotal$ introduced in Eqs. (9)–(14), and the transition frequency $\omega KaLa=(\u03f5Ka\u2212\u03f5La)/\u210f$ between exciton states $Ka$ and $La$.

Note that the site energy correlations are taken into account by the off-diagonal elements of the spectral density. It was found previously^{9} that these correlations have practically no influence on exciton relaxation in the FMO monomers. What remains to be investigated is whether they can influence intermonomer energy transfer in the FMO trimer. For this purpose, we extend the standard generalized Förster theory in the following.

#### 2. Interdomain transfer: Generalized Förster theory with correlations

We assume a fast intradomain relaxation leading to equilibration of the excitation energy within the domains prior to energy transfer between different domains. Hence, the overall rate constant of energy transfer between domain *a* and domain *b* can be written as

with the Boltzmann factor

Here, the generalized Förster rate constant $kMa\u2192NbGF$ for transfer between exciton state $Ma$ of domain *a* and $Nb$ of domain *b*, obtained as described in Appendix A, is given as

The interdomain electronic coupling $VMaNb$ is given in Eq. (8), while the function $FMaNb(t)$ reads

where

contains the fluctuations of the intradomain site energies and

contains the correlations between site energy fluctuations in different exciton domains (*a*$\u2260$*b*), which are described by the intermolecular spectral density $Jmanb(\omega )$ in Eq. (10).

The exciton relaxation-induced inverse dephasing time $\tau Mc\u22121$ is given as

and the shifted transition frequencies are defined as

with

Here, $\omega KC$ is the vertical excitation frequency of exciton state $Kc$, obtained from the diagonalization of the exciton matrix of domain *c*, and $E\lambda (Kc)$ is the reorganization energy introduced above. It is related to the spectral density $Jmcnctotal$ [Eq. (9)] via

Δ*ω*_{c} in Eq. (25) is a small additional shift caused by the off-diagonal elements of the exciton-vibrational coupling,

where *℘* denotes the principal part of the integral. Note that the intermolecular spectral density $Jmcnctotal(\omega )$ in Eqs. (15), (21), and (26) contains, both the site energy fluctuations (*m*_{c} = *n*_{c}) and their correlations (*m*_{c} ≠ *n*_{c}), whereas the interdomain spectral density $Jmanb(\omega )$ in Eq. (22) just describes the interdomain correlations in site energy fluctuations of pigment *m*_{a} in domain *a* and pigment *n*_{b} in domain *b* (*b* ≠ *a*).

##### a. Limit of uncorrelated fluctuation of site energies.

In the limit of uncorrelated fluctuations of site energies in different domains, we can set $GMaNb(t)=0$ and obtain the well-known result for the rate constant^{57,63–66}

that contains the overlap integral between the normalized absorption lineshape function $DNb(\omega )$ of the acceptor and the normalized fluorescence lineshape function $DMa\u2032(\omega )$ of the donor, which are given as

and

respectively.

##### b. Classical limit of nuclear motion.

In the following, we discuss a classical description of nuclear motion. In this high-temperature limit, the Bose–Einstein distribution function in Eq. (16) is approximated as *n*(*ω*) ≈ *k*_{B}*T*/*ℏω*. Neglecting the exciton relaxation-induced dephasing ($\tau Ma\u22121=\tau Nb\u22121=0$) in Eq. (19) and using a short-time approximation [cos(*ωt*) ≈ 1 − *ω*^{2}*t*^{2}/2, sin(*ωt*) ≈ *ωt*] in the function $FMaNb(t)\u2212FMaNb(0)$ in this equation and Eq. (A8) results in

where the reorganization energy is

that contains the spectral density $JMaNb(\omega )$ in Eqs. (A9) and (A10). Performing the time integration in Eq. (19) finally results in the semiclassical generalized Förster rate constant, which is formally identical to the rate constant obtained in Marcus theory of non-adiabatic electron transfer,^{60,67}

The correlations in site energy fluctuations are contained in the spectral density $JMaNb(\omega )$ in Eq. (A10), where the intradomain correlations are described by the off-diagonal elements $Jmcnc(\omega )$, and the correlations between different domains *a* and *b* are given by $Jmanb(\omega )$.

##### c. Limit of localized excitations: Standard Förster theory.

In order to evaluate how the intermonomer excitation energy transfer is influenced by the intramonomer exciton delocalization, we consider the limit of localized excitations. In this case, in analogy to Eq. (17), we obtain a rate constant

with the Boltzmann factor

that takes into account fast equilibration of the population of localized excited states in FMO monomer *a*, where $\u210f\omega ma$ is the site energy of pigment *m*_{a}. The Förster theory rate constant $kma\u2192nbF$ between localized excited states $ma$ in domain *a* and $nb$ in domain *b* reads^{60,68}

where the normalized absorption lineshape function $Dnb(\omega )$ of the acceptor is

and the normalized fluorescence lineshape function of the $Dma\u2032(\omega )$ of the donor is

where *τ*^{pd} is a pure dephasing time. The time-dependent function $Gkc(t)$ of the local excited state $kc$ is given as

where $Jkckctotal(\omega )$ describes the modulation of the site energy of this state. The frequency $\omega \u0303kc$ in Eqs. (37) and (38) is the 0–0 local transition frequency of pigment *k* in domain *c* and is related to the local vertical transition frequency $\omega kc$ by

where $E\lambda (kc)$ is the local reorganization energy of the exciton-vibrational coupling

For the numerical implementation, we use the following expression for the rate constant [obtained by inserting Eqs. (37)–(39) into Eq. (36)]:

An important difference between the above rate constant and the rate constant obtained in the generalized Förster theory in Eq. (28) is given by the prefactor. The excitonic coupling $Vmanb$ between local excited states in Förster theory is replaced by that between delocalized excited states $VMaNb$ in the generalized Förster theory. According to Eq. (8), $VMaNb$ is related to the $Vmanb$ by the coefficients of the intermonomer exciton states. Since the latter can have positive as well as negative signs, exciton delocalization can lead to an enhancement as well as a decrease in the coupling corresponding to constructive and destructive interference effects, respectively.

#### 3. Exciton state lifetime

The lifetime $TMd$ of an exciton state $Md$ is determined by the rate constants of transfer from this state to any other state $Nd$ in the same domain and to any state $Na$ in the other domains *a*≠*d*. We define the frequency-resolved average inverse lifetime *T*(*ω*)^{−1} as

where ⟨$\cdots \u2009$⟩_{dis} denotes an average over static disorder in site energies and

is the density of exciton states.

For analysis purposes, we also introduce functions *T*_{intra}(*ω*) and *T*_{inter}(*ω*) in which only the intramonomer rate constants $kMd\u2192Nd$ or only the intermonomer rate constants $kMd\u2192Na$ (*a* ≠ *d*), respectively, are taken into account on the rhs in Eq. (43).

A variant of the above quantities denoted as *T*′(*ω*) and *DoE*′(*ω*) is used at cryogenic temperatures, where uphill energy transfer is frozen out leading to an infinite lifetime of the lowest exciton state. Therefore, the lowest exciton state of every monomer, in the case of $Tintra\u2032$(*ω*), or of the FMO trimer, in the case of $Tinter\u2032$(*ω*), is excluded from the sum over *M*_{d} on the rhs of Eqs. (43) and (44). In another variant, termed *DoE*′^{r}(*ω*) and *T*′^{r}(*ω*), in addition to leaving out the lowest exciton state(s), instead of the bare exciton frequencies $\omega Md$, the renormalized frequencies $\omega \u0303Md$ [Eq. (25)], which contain a shift by the exciton-vibrational coupling, are used in the rhs of Eqs. (43) and (44).

## III. COMPUTATIONAL DETAILS

Computations are based on the crystal structure of the FMO protein from *P. aestuarii* by Tronrund *et al.*^{30} at a resolution of 1.3 Å that includes all eight BChl pigments per monomer (*holo* form). A starting structure for the geometry optimization was constructed from protein data bank (PDB) file 3EOJ. The protein was included from Thr 8, whose N-terminus was patched with an acetyl group, to the negatively charged C-terminal Gln 366. Hydrogen atoms and missing heavy atoms of the protein (see PDB file 3EOJ) were added by using CHARMM35b3^{69} and the CHARMM22 force field.^{70} All histidine side chains were modeled as neutral except for His 12, which was positively charged. The remaining titratable residues were calculated to be in their standard protonation state (see Appendix B) and modeled accordingly.

BChl *a* pigments were complete in the crystal structure except for BChl 8, where the phytyl chain is replaced by a methyl group.^{30} All water molecules were removed except for the axial ligand of BChl 2 and the hydrogen bonding partner of the 3-acetyl group of BChl 1. The force field of the BChl *a* pigments was adopted from the work of Ceccarelli *et al.*^{71} for the use in CHARMM except for the atomic partial charges (APCs). The latter were obtained from quantum chemical (QC) computations employing QChem.^{72} In QC computations, the geometry of a BChl *a* with the phytyl chain truncated to a methyl group was optimized *in vacuo* with the B3LYP exchange-correlation functional and a 6-31G(d,p) basis set. Excited state calculations were performed on this optimized structure employing time-dependent density functional theory in the random phase approximation with the same basis set and the CAM-B3LYP^{73} exchange-correlation functional, which is superior to B3LYP for excited states.^{74} APCs of the ground and first excited state of BChl *a* were obtained by fitting the electrostatic potential (ESP)^{75,76} of the CAM-B3LYP-derived charge density of the two electronic states. The numerical values of these charges are given in the supplementary material. The ground state APCs were used in all CHARMM computations supplementing the force field by Ceccarelli *et al.*^{71} The ground and excited state APCs were tested in a calculation of site energies employing our Poisson Boltzmann/Quantum chemical (PBQC) method,^{4} revealing a good correlation with the reference values^{7} that describe the optical spectra ( Appendix B). Additionally, the ground state APCs are very similar in magnitude to those obtained in a recent study^{77} by using the CHARMM-compatible CGenFF tool^{78} and, thus, are expected to provide a good description of a polarized BChl *a* molecule in a condensed matter environment according to the CHARMM philosophy.

The geometry of the FMO trimer was optimized with CHARMM using the adopted basis Newton–Raphson (ABNR)^{79} method. CHARMM was also used for the NMA, resulting in 59 373 eigenfrequencies and normal mode vectors. The first six frequencies, corresponding to translation and rotation of the whole FMO trimer, are practically zero (<10^{−3} cm^{−1}). All other eigenfrequencies *ω*_{ξ} are positive, as expected for a stable minimum. The eigenvectors and eigenfrequencies were used in conjunction with the optimized (equilibrium) coordinates of atoms (nuclei), the ground and excited state APCs of BChl *a*, and the ground state APCs of the protein from the CHARMM22 force field to calculate the coupling constants $g\xi (m)$ [Eq. (13)] entering the intermolecular spectral density *J*_{mk}(*ω*) [Eq. (10)] together with the eigenfrequencies *ω*_{ξ} of the normal modes.

## IV. NUMERICAL RESULTS ON THE FMO PROTEIN

### A. Spectral density

The diagonal part *J*_{mm} of the intermolecular spectral density of the eight pigments in the monomeric subunit of the FMO protein, averaged over the three equivalent sites of the trimer, is shown in Fig. 3. It exhibits a strongly asymmetric shape with a sharp rise at small frequencies toward a maximum between 10 cm^{−1} and 20 cm^{−1} and a slower decay for larger frequencies reaching values close to zero at a vibrational frequency of roughly 100 cm^{−1}. The integral over the diagonal part *J*_{mm} of the spectral density, the Huang–Rhys factor

varies between *S*_{5} = 0.15 and *S*_{2} = 0.60.

These Huang–Rhys factors were averaged over the three equivalent sites in the FMO trimer. The individual Huang–Rhys factors of the monomers exhibit a similar trend as the averaged values with some deviations between the equivalent sites in the three FMO-monomers (Fig. 4). Obviously, the geometry optimization of the FMO protein leads to some deviations from the perfect C3 symmetry of the crystal structure.

The Huang–Rhys factors are mostly influenced by the direct protein environments of the pigments, as shown in Fig. 5, where the dependence of *S*_{m} on a cutoff distance *R*_{c} is shown. For a given value of *R*_{c}, only those amino acids (and other components of the complex) were included in the calculation of *S*_{m} that have at least one atom with a smaller distance than *R*_{c} to at least one atom of BChl *m*. For most of the sites *m*, the Huang–Rhys factor *S*_{m} is converged for *R*_{c} ≥ 10 Å–15 Å, demonstrating that electrostatic interactions with the local protein environment are dominant. In contrast, for the two sites *m* = 2 and 6 with the largest *S*_{m}, the latter is influenced by long-range interactions that extend beyond a FMO monomer. *S*_{6} exhibits a sharp initial rise, a decrease for intermediate *R*_{c}, and a slight but steady increase for large *R*_{c} > 20 Å. The largest single amino acid contribution of 0.2 to the initial rise of *S*_{6} up to 0.55 is due to the polar Gln 143. *S*_{2} increases steadily up to *R*_{c} ≈ 40 Å and afterward slightly decreases until convergence at *R*_{c} ≥ 60 Å. The Huang–Rhys factor *S*_{3} = 0.31 of the low-energy site BChl 3, which dominates the lowest exciton state of the FMO monomers, agrees well with the experimental estimates *S* = 0.3 of the intermolecular Huang–Rhys factor from fluorescence line narrowing spectra.^{62} The intermolecular spectral density *J*_{33}(*ω*) obtained for BChl 3 is similar to the average of the diagonal parts of the spectral densities of all pigments in the FMO trimer,

A comparison between *J*_{33}(*ω*) and *J*_{NMA,av}(*ω*) as well as with the experimental spectral density *J*_{FLN}(*ω*) extracted from fluorescence line narrowing experiments^{62} is shown in Fig. 2. *J*_{FLN}(*ω*) is smaller at frequencies *ω* < 25 cm^{−1} and larger at *ω* >25 cm^{−1}. The latter deviations are due to the fact that the experimental spectral density also includes intramolecular contributions that are due to vibrations of the pigments. Huang–Rhys factors of these intrapigment modes have been estimated from experimental ΔFLN data^{35} and are shown as green bars in Fig. 2. The large amplitude of the calculated NMA spectral density at low frequencies could be due to the harmonic approximation that cannot describe soft and anharmonic motion expected at very low frequencies. As discussed later in more detail, the local reorganization energies of the intermolecular spectral density obtained from the present NMA/CDC approach qualitatively agree with those obtained from combined molecular dynamics/CDC calculations,^{10} providing some evidence that, overall, the anharmonic effects are not critical.

An advantage of a microscopic calculation of the spectral density is that correlations in site energy fluctuations can be revealed. These correlations are contained in the off-diagonal parts of *J*_{mk}(*ω*) shown in Fig. 6 for selected pigments located in the same (left half of this figure) and in different (right half) FMO monomers. In contrast to the positive diagonal parts *J*_{mm}(*ω*) (Fig. 3), the off-diagonal elements *J*_{mk}(*ω*) can also become negative at certain frequencies, corresponding to a situation where the normal modes at these frequencies lead to an anticorrelated fluctuation of site energies. Due to the broken C3 symmetry, discussed already above, the correlations are different in the three FMO monomers and also between the monomers. Interestingly, in Fig. 6, the amplitudes of the intramonomer correlations are comparable to those of the intermonomer correlations despite the much larger interpigment distances in the latter case (see Fig. 1). These amplitudes are somewhat smaller but still in the same order of magnitude than those of the diagonal parts *J*_{mm}(*ω*) in Fig. 3. A difference, however, is that the correlations *J*_{mk}(*ω*) with increasing *ω* decrease steeper to zero than the diagonal parts *J*_{mm}(*ω*). Whereas the former are practically zero for *ω* > 40 cm^{−1}, the latter still exhibit non-zero values up to a frequency of 100 cm^{−1}. In order to have a global measure of the correlations, we introduce a generalized Huang–Rhys factor

and a generalized absolute Huang–Rhys factor

representing a hypothetical system, where all anticorrelations were transformed into correlations. As seen in Fig. 7, the latter Huang–Rhys factor shows no dependence on the interpigment distance *R*_{mk}, reflecting the fact of equal amplitudes of intra- and intermonomer correlations (Fig. 6). In contrast, the amplitude of *S*_{mk} decreases with the increasing interpigment distance (Fig. 7). Obviously, the compensations between correlations and anticorrelations in the site energy fluctuations due to different normal modes are more complete for pigments at large distance.

### B. Energy transfer

#### 1. Parameters

In our calculations of energy transfer, presented below, we use the site energies of the FMO protein determined earlier from CDC calculations and a refinement fit of optical spectra.^{7} The excitonic couplings were obtained by using a point-dipole approximation (PDA) and an effective transition dipole moment of 29.8 D^{2}, as determined earlier^{5} from electrostatic calculations. Note that the validity of the PDA for the FMO protein has been demonstrated.^{80} The transition dipole of BChl *a* is assumed to be oriented in the direction of a line connecting the N_{B} and N_{D} atoms (in PDB nomenclature) and placed in the center between these two atoms. The atomic coordinates were taken from the crystal structure of the FMO trimer.^{30} The numerical values of the excitonic couplings are given together with that of the site energies in Appendix C.

Static disorder in site energies is taken into account by assigning a Gaussian distribution function with a full width at half maximum (FWHF) of 100 cm^{−1} to every site,^{3} assuming that there is no correlation between different sites. For every random realization of disorder, generated from the above distribution functions, rate constants of energy transfer are calculated and the distribution and the average values determined. For room temperature calculations, an ensemble of 10 000 sets of site energies is sufficient to obtain convergent results, whereas a 100 times larger ensemble had to be used for cryogenic temperatures. In the calculations of rate constants, we use different spectral densities of the exciton-vibrational coupling, as explained in the following and summarized in Table I. The spectral density *J*_{NMA} corresponds to the intermolecular spectral density *J*_{mk}(*ω*) obtained from the NMA [Eq. (10)]. In *J*_{NMA} (no corr.), only the diagonal part *J*_{mm}(*ω*) of the NMA spectral density is taken into account, and the correlations contained in the off-diagonal parts are neglected. In *J*_{NMA} (pos. corr.), all anticorrelations are turned into correlations of the same amplitude. *J*_{NMA,av} describes the average diagonal part of the NMA spectral density [Eq. (46)], *J*^{intra} is the intrapigment contribution to the spectral density [Eq. (14), green bars in Fig. 2], and *J*_{FLN} is the experimental spectral density extracted from FLN spectra of the FMO protein containing intrapigment as well as intermolecular contributions (black line in Fig. 2).

Spectral density . | Content . |
---|---|

J_{NMA} | J_{mk}(ω), Eq. (10) |

J_{NMA} (no corr.) | J_{mk}(ω) = δ_{m,k}J_{mm}(ω) |

J_{NMA} (pos. corr.) | J_{mk}(ω) → |J_{mk}(ω)| |

J_{NMA,av} | Equation (46) |

J^{intra} | Eq. (14), Fig. 2 |

J_{FLN} | Extracted from FLN spectra,^{62} Fig. 2 |

#### 2. Physiological temperatures

We start with a description of intramonomer exciton relaxation and intermonomer transfer at room temperature (300 K). The density of exciton states *DoE*(*ω*) [Eq. (44)] is shown in Fig. 8 together with the frequency-resolved lifetime *T*(*ω*) [Eq. (43)] calculated with the different spectral densities described above. Including only the intermolecular part of the spectral density results in lifetimes around 200 fs in the central part of the spectrum that increase to more than a ps at the low and high-energy wings. The inclusions of correlations practically has no influence, and also neglecting the site dependence of the spectral density (using *J*_{NMA,av}) leads to very similar results. However, including the high-frequency intrapigment spectral density *J*^{intra}(*ω*) decreased the lifetimes to less than 100 fs in the center of the spectrum and to less than 1 ps at the low and high-energy wings. The lifetimes in Fig. 8 are dominated by intramonomer exciton relaxation, that is, excluding the intermonomer transfer does not change the lifetimes of exciton states.

The distribution of intermonomer transfer times $(ka\u2192b)\u22121$, taking into account the fast intramonomer equilibration [Eq. (17)], are shown in Fig. 9, calculated for different spectral densities. Similar distribution functions are obtained with and without including the intrapigment (high-frequency) part of the spectral density. The maxima of these distribution functions all occur around an intermonomer transfer time of 10 ps that is about two orders of magnitude larger than the typical intramonomer exciton relaxation times (Fig. 8). Including the high-frequency intrapigment contribution to the spectral density leads to a somewhat sharper distribution function but does not change the average dynamics. Hence, the intermonomer transfer is dominated by the low-frequency intermolecular spectral density. Interestingly, practically the same distribution function is obtained with and without including the correlations in site energy fluctuations contained in the off-diagonal elements of the spectral density *J*_{mk}(*ω*), discussed above. In order to see whether the compensation between correlations *J*_{mk}(*ω*) > 0 and anticorrelations *J*_{mk}(*ω*) < 0 can be responsible for this result, we recalculated the distribution function for a hypothetical spectral density, where all anticorrelations were transformed into correlations. For this *J*_{NMA} (pos. corr.), we again find very little deviations from the result obtained for the original spectral density *J*_{NMA}. Therefore, another factor must be responsible for the small effect of the site energy correlations on the transfer times. In addition, the averaged spectral density *J*_{NMA,av} gives very similar results. Even a classical approximation of nuclear motion in the calculation of the intermonomer rate constants $kMa\u2192NbscGF$ [Eq. (33)] provides a qualitatively correct description of the intermonomer transfer, with a distribution function that is slightly shifted toward shorter transfer times (the blue solid line in Fig. 9).

#### 3. Cryogenic temperatures

At cryogenic temperatures, uphill energy transfer between the exciton states in the FMO monomers is frozen out and the lifetime of the lowest exciton state is determined by the intermonomer transfer, which, however, occurs on a much longer time scale than intramonomer exciton relaxation. We, therefore, analyze the lifetimes of the lowest exciton state in the FMO monomers separately. We start with the intramonomer lifetimes of the exciton states, excluding the lowest state in each monomer. The density of these exciton states together with their lifetime $Tintra\u2032$(*ω*) at 77 K is shown in Fig. 10. Similar to the room temperature results (Fig. 8), the high-frequency intrapigment contributions to the spectral density are important for a fast sub-picosecond exciton equilibration. The correlations in the intermolecular part of the spectral density have practically no effect on the lifetimes. Compared to the room temperature lifetimes (Fig. 8), those at 77 K are somewhat larger, in particular, those of the low-energy exciton states.

The lifetime of the lowest exciton state of the FMO trimer is not limited by energy transfer since the latter would be uphill in energy and is frozen out at 77 K. However, the second and third lowest exciton states of the FMO trimer, which are the lowest exciton states in two FMO monomers, can transfer downhill in energy to the lowest exciton state in the trimer. In Fig. 11, we have analyzed the respective lifetimes at 4 K and also provide the density of exciton states of the three lowest exciton states, the second and third lowest exciton states, and the third lowest exciton state of the FMO trimer, calculated with two different spectral densities *J*_{NMA}(*ω*) and *J*_{FLN}(*ω*). We have included the renormalization of the exciton energies by the exciton-vibrational coupling in order to get as close as possible to the excitation energies that enter the absorption lineshape function [Eq. (29)]. Therefore, the density of exciton states depends now also on the spectral density of the exciton-vibrational coupling. The frequency-resolved lifetime *T*′^{r}(*ω*) of exciton states is constant at about 15 ps on the high-energy half of the lowest exciton band and linearly increases up to 130 ps toward lower energies. Such an increase was indeed measured in hole burning spectroscopy,^{47} as discussed above and denoted by two black circles in the upper parts of Fig. 11. An explanation of this behavior is given in terms of the distribution of downhill intermonomer transfer times between the different low-energy exciton states in Fig. 12. The distribution functions of the transfer between the third and the second exciton state (3 → 2) have a larger amplitude at short transfer times than those of the 3 → 1 and 2 → 1 transfers, with average transfer times of 21 ps, 141 ps, and 45 ps, respectively (left upper part of Fig. 12). Therefore, the lifetime of the third exciton state is determined by the fast transfer to the second exciton state, and that of the second exciton state is determined by the slow transfer to the lowest exciton state. When going to lower energies in the lowest exciton band, it becomes less and less likely to find a FMO trimer that has two more exciton states at lower energies, and hence, the fast transfer 3 → 2 will less and less contribute to the lifetime at lower energies, explaining the increase in the lifetime. This effect is lost if the lowest exciton state gets more localized. We have increased the energy gap between the low-energy pigment BChl 3 and the remaining pigments by 300 cm^{−1} and repeated the calculations of the distribution functions of the inverse rate constants and the frequency-resolved exciton state lifetime (middle part in Fig. 12). Now, the distribution functions for the 3 → 2 and the 2 → 1 transfer are almost equal, resulting in very similar average inverse rate constants of 19 ps and 17 ps, respectively. This equality gives rise to a much less frequency-dependent lifetime across the whole low-energy exciton band (middle right part of Fig. 12). For a complete localization of excited states, that is, using standard Förster theory to describe the transfer between localized excited states in different monomers, we find the fastest transfer and practically identical distribution functions for the 3 → 2 and 2 → 1 transfer with an identical average inverse rate constant of 8 ps, resulting in a frequency-independent lifetime.

Finally, we have investigated which pigments essentially contribute to the transfer between the low-energy exciton states of FMO monomers. Besides the lowest-energy pigment BChl 3, obvious candidates are BChls 4, 2, and 7, which are the pigments with the next higher transition energies (see Table II in Appendix C) and relatively large excitonic couplings to pigments in another FMO monomer (Table III in Appendix C). Including just BChls 3 and 4 in the calculation of intermonomer energy transfer lifetimes gives rise to an inverted wavelength dependence as compared to the case, where all pigments are included (blue and black line in Fig. 13, respectively). Taking into account, in addition, BChls 2 and 7 results in a qualitatively correct wavelength dependence (orange line). A more quantitative agreement with the full calculation results, if besides the low-energy BChls 3, 4, 2, and 7 also BChl 5 (which is strongly coupled to BChl 2′ of the neighboring monomer) is included (brown line).

## V. DISCUSSION

### A. Spectral density

The present approach to split the spectral density into an intermolecular part, obtained from a combination of the NMA and the charge density coupling (CDC) method,^{9} and a site-independent intrapigment part, obtained from FLN experiments,^{35} is one of the simplest possible ways to avoid the geometry mismatch problem^{10,15,16,36,37} that one is facing in traditional QM/MM approaches that combine classical molecular dynamics simulations with quantum chemical calculations of transition energy shifts of the pigments. Since the CDC method^{5} just determines the intermolecular pigment–protein Coulomb interactions, it does not suffer from slight distortions of the nuclear coordinates. A workaround, related in spirit, was suggested by Shi and co-workers,^{36} combining the CDC method with a classical MD simulation to obtain the intermolecular part of the spectral density and performing a quantum mechanical normal mode analysis on the isolated, geometry-optimized pigment to obtain a site-independent intrapigment part. An advantage of a MD simulation is that it includes the anharmonic motion of nuclei, whereas our present NMA is able to resolve correlations in site energy calculations at very low vibrational frequencies that are difficult to sample with MD simulations. Concerning the site-independent Huang–Rhys factors of the intrapigment modes that were extracted from FLN spectra of the FMO protein,^{35} we cannot exclude that these factors are biased by the protein environment of the low-energy pigment BChl 3, which dominates the lowest exciton state seen in the low-temperature fluorescence spectrum. On the other hand, a QM NMA on an isolated pigment depends on the limits of the QM method used, as, e.g., on the exchange-correlation functional in (TD)DFT calculations.^{15} Coker and co-workers^{10,15} found a way to take into account the site dependence of the intrapigment part of the spectral density, avoiding the geometry mismatch problem. They performed a local QM geometry optimization of the ground state potential energy surface of the pigments along the classical nuclear trajectories obtained from MD simulations and performed a QM NMA as a basis for the calculation of intramolecular Huang–Rhys factors of the pigments. In their careful analysis, using different XC functionals, they arrive at the conclusion (based on results presented in Fig. S5 of the supplementary material of Ref. 15) that for XC-(hybrid) functionals with small amounts or no exact exchange, the intramolecular spectral density of the different BChls is very similar and does not depend critically on the XC functional. However, with increasing amount of exact exchange, the intrapigment spectral density becomes both site-dependent and also depends on the XC-functional. An alternative solution to the geometry mismatch problem was provided by Rhee^{81,82} and co-workers and by Saito and co-workers^{16,37} using an interpolation between the MM and the QM potential energy surfaces of the BChls. Whereas Saito *et al.*,^{16,37} using the CAM-B3LYP XC-functional, find rather strong variations between sites, the spectral densities calculated by Kim *et al.*^{82} with the B3LYP XC-functional are characterized by very similar reorganization energies and Huang–Rhys factors. There is hope^{15} that the optimally tuned range-separated hybrid functionals^{83,84} will help to settle this issue. At present, it is undecided whether the present assumption of a site-independent intrapigment spectral density is justified.

Concerning the intermolecular contributions to the spectral density, Coker and co-workers,^{10,15} combining the CDC method^{5} with MD simulations, arrived at similar local reorganization energies $E\lambda (k)$ for the BChl pigments in the different binding sites *k*, as seen in Table II ( Appendix C). The largest difference is obtained for BChl 8, which is solvent-exposed. Most likely, interaction with water molecules, not included in the present NMA, but in the MD simulation, is responsible for the three times larger reorganization energy obtained with the latter. Both methods agree on the large intermolecular reorganization energy of BChl 2, which has a water ligand and a couple of additional water molecules in its neighborhood. The analysis of Rivera *et al.*^{10} suggests that about half of the enlarged reorganization energy of BChl 2 is due to the water molecules and the remaining half is due to interaction with polar amino acid side chains of the protein. Since in our analysis only one water molecule, the axial ligand of BChl 2, was taken into account (besides a second water molecule that is hydrogen bonded to BChl 1) the good agreement between the two reorganization energies of BChl 2 seems to be somewhat fortuitous. A remarkable feature of BChl 2 in our calculations is the effect of long-range electrostatic couplings over distances as large as 40 Å (Fig. 5) that increase the Huang–Rhys factor (and the related reorganization energy) of this pigment.

### B. Correlations in site energy fluctuations and energy transfer

A main purpose of the present work was to quantify the correlations in site energy fluctuations in FMO trimers and to investigate their role in intermonomer energy transfer. So far, there are three studies on correlations,^{9,85,86} which investigated the monomeric subunit of the FMO protein. In a QM/MM approach,^{85} no correlations were found, whereas our earlier NMA in combination with the CDC method revealed correlations at very low vibrational frequencies.^{9} Probably, the MD simulations did not resolve these correlations because of the finite simulation times or because the intermolecular site energy correlations were overshadowed by artificially enhanced uncorrelated intrapigment site energy fluctuations caused by the geometry mismatch problem. Cross correlations between site energy fluctuations were reported in another QM/MM study,^{86} despite the geometry mismatch problem and short (40 ps) simulation times. These correlations were found to have practically no influence on intramonomer exciton relaxation, in agreement with our earlier^{9} and our present simulations (Fig. 10).

Interestingly, we find similar amplitudes of the intramonomer and the intermonomer correlations in site energy fluctuations (Fig. 6), despite the much larger distances between pigments in different monomers (Fig. 1). Obviously, the low-frequency normal modes of the FMO protein are delocalized over the entire trimer, causing the correlations in site energy fluctuations. Do they have any functional relevance? No, these correlations have practically no influence on intermonomer energy transfer (Fig. 9).

An important result of our earlier work^{9} is that exciton relaxation were most affected if the correlations in site energy fluctuations, which are contained in the off-diagonal parts of the spectral density *J*_{mk}(*ω*), would be related to the respective diagonal parts *J*_{mm} and *J*_{kk}, describing the site energy fluctuations, in the following way:

In this hypothetical case, using the equality $\u2211kacka(Ma)cka(Na)=\delta Ma,Na$, it is seen that the Redfield relaxation rate constant $kMa\u2192NaRedf$ in Eq. (15) becomes zero. In addition, it was found that coherences between different exciton states would live very long.^{9,25} These results led us to conclude that nature does not protect inter-exciton coherences since excitons would not relax.

Introducing Eq. (49) into the expression for the generalized Förster rate constant $kMa\u2192NbGF$ [Eq. (19)], noting that the dephasing times $\tau Ma$ and $\tau Nb$ are infinite, because of the missing intermonomer exciton relaxation, and the function $FMaNb(t)$ in Eqs. (20)–(22) is zero, we obtain

This equation is simply Fermi’s Golden rule, where the only effect of the exciton-vibrational coupling left is the small renormalization of the exciton energies [Eq. (25)]. For the present site energies and excitonic couplings, an inverse average intermonomer rate constant of 24 ps is obtained, where the average was taken with respect to disorder in site energies and the different intramonomer exciton states. This time constant is more than twice as large as the ≈10 ps transfer time obtained for the realistic intermolecular spectral density at room temperature (Fig. 9). Hence, the effect of the hypothetical correlations in Eq. (49) on the intermonomer transfer is not as dramatic as on the intramonomer exciton relaxation but still present. In the hypothetical case that intramonomer exciton relaxation is still fast, the disorder averaged intermonomer rate constant *k*_{a→b} [Eq. (17)], using the $kMa\u2192NbGF$ for perfectly correlated site energies [Eq. (50)], corresponds to a time constant of 5 ps, which is even faster than that of the uncorrelated transfer. However, nature cannot take advantage of this effect because excitons would not relax in the case of perfect correlations.

What still needs to be understood is why the real correlations in Fig. 6 are obviously so far away from the hypothetical case in Eq. (49). One difference is that the hypothetical *J*_{mk}(*ω*) is positive at all frequencies, whereas the real *J*_{mk}(*ω*) can be positive at some frequencies and negative at others. However, inverting the negative values in the latter has practically no influence on the calculated rate constants for the intermonomer transfer (Fig. 9). Hence, the explanation must be related to the different amplitudes and spectral shapes of the diagonal (Fig. 3) and off-diagonal (Fig. 6) parts of the spectral density *J*_{mk}(*ω*). The off-diagonal parts with increasing frequency decrease faster to zero than the diagonal parts, reflecting the fact that correlations require delocalized vibrational normal modes, whereas the diagonal parts *J*_{mm}(*ω*) are influenced strongest by the intermediate environments of the pigments (Fig. 5). Hence, high-frequency normal modes that are localized in the close environment of the pigments have an influence on the diagonal parts *J*_{mm}(*ω*) but do not affect the off-diagonal parts *J*_{mk}(*ω*). Therefore, Eq. (49) cannot be fulfilled, in particular, at large intermolecular frequencies. Since anharmonic vibrational effects, neglected in the present work, are expected to contribute mainly at low frequencies, it is likely that these effects will not work in favor of Eq. (49) either. Nevertheless, it is of general interest to clarify the role of anharmonic nuclear motion in the correlation of site energy fluctuations.

### C. Intramonomer exciton relaxation

Intramonomer exciton relaxation in the FMO protein occurs on an ultrafast subpicosecond timescale both at room temperature (Fig. 8) and at cryogenic temperatures (Fig. 10). In both cases, the exciton relaxation is accelerated by the dissipation of high-frequency intrapigment vibrational quanta, reflecting the fact that the differences in exciton energies are in the 0 cm^{−1}–500 cm^{−1} range.^{3,9} Good agreement between calculated lifetimes with those estimated in 2D experiments on *C. tepidum* at 77 K is obtained except for the bottleneck around 810 nm, which is much more pronounced in the experiment (Fig. 10). The deviations could be due to the fact that the present calculations were performed with the Hamiltonian of *P. aestuarii*, whereas the experiments were done on *C. tepidum*. Although both structures are very similar, their optical spectra show some differences^{3} that still need to be understood. Another explanation for the long experimental lifetime around 810 nm could be that the pigments contributing to the underlying exciton state (BChls 1 and/or 2)^{3,25} exhibit a weaker intrapigment exciton-vibrational coupling and therefore the respective exciton state a longer lifetime (Fig. 10). Recent QM/MM calculations by Saito and co-workers,^{16,37} which report a particularly strong exciton-vibrational coupling of BChl 2, seem to suggest that a candidate for weak intrapigment coupling could be BChl 1. On the other hand, as discussed above, the QM/MM calculations of Kim *et al.*^{82} obtain only very small variations in Huang–Rhys factors of the different pigments. Saito *et al.*^{16} demonstrated that their variation in local reorganization energies optimizes exciton relaxation. It is, however, difficult to imagine that nature actively worked on this optimization, considering that the transfer from the baseplate into the FMO protein and the transfer from the FMO protein to the reaction center complex occur on orders of magnitudes slower timescales.^{87}

### D. Intermonomer energy transfer

Intermonomer energy transfer is calculated to occur with a time constant of about 10 ps at room temperature (Fig. 9). As noted earlier, this transfer is hard, if not impossible, to detect experimentally since the lifetime of the exciton states including the lowest one is determined by intramonomer exciton transfer. Even the uphill transfer in the FMO monomer is more than an order of magnitude faster (Fig. 8) than the intermonomer transfer. Most likely, the time constant of 1.4 ps–2 ps suggested in an early room temperature pump–probe study^{38} to reflect intermonomer transfer has a different origin. According to the present calculations, this transfer is about five times slower.

In contrast to the intramonomer transfer, the intermonomer transfer is governed by the couplings of excitons to low-frequency intermolecular vibrations (Fig. 9), suggesting that the most important intermonomer channels involve exciton states at close energies in different FMO monomers. Because of the small vibrational frequencies involved in the transfer, the classical description of nuclear motion works so well at room temperature (Fig. 9).

At cryogenic temperatures, the lifetimes of two of the three exciton states in the lowest exciton band of the trimer are limited by intermonomer transfer since intramonomer transfer would be uphill in energy and is, therefore, frozen out. In good agreement with estimates from hole burning data,^{47} we obtain a lifetime of 15 ps–25 ps for the high-energy half of the low-energy exciton band and an increase to more than 100 ps at the low-energy edge. The origin of the increase lies in the intramonomer delocalization of the lowest exciton state in the FMO monomers. Interestingly, the minimal model to describe this effect includes not just the two low-energy pigments BChls 3 and 4 but also BChls 2, 5, and 7 (Fig. 13).

Due to this delocalization, the square of the excitonic coupling entering the rate constant [Eqs. (19) and (8)] between the third and the second exciton state of the FMO trimer on average is about a factor of 2.4 larger than that between the second and the first exciton state, giving rise to the different rate constants that explain the increase in lifetime (Fig. 12). Upon localizing the lowest states of the FMO monomers, this effect is lost, and the lifetime becomes frequency-independent (Figs. 12 and 13). Interestingly, the intermonomer transfer is faster in this hypothetical limit. So, exciton delocalization in the FMO monomers slows down intermonomer transfer, in contrast to the results of Jankowiak and co-workers,^{49,50} who, however, did not focus on the present effect and used a different spectral density in their Förster and in their generalized Förster calculations.

We note that an opposite effect was found in the LH2 complex of purple bacteria, where standard Förster theory gave an about five times too small rate constant for the transfer between the B800 and the B850 states of the antenna.^{88} Generalized Förster theory taking into account the delocalization of the B850 states explained^{64–66} the experimental result. Optically dark exciton states of the B850 manifold accept excitation energy from the B800 states, an effect that cannot be described by Förster theory, where only bright states contribute to the transfer.

The FMO trimer represents an example, where quantum mechanic delocalization of excited states acts in the opposite direction, slowing down the transfer. Obviously, the light-harvesting apparatus of green sulfur bacteria can tolerate this effect. At room temperature, we find that Förster theory would result in an average intermonomer transfer time of 6 ps, which is roughly half the time we obtain by using generalized Förster theory. Any quantitative estimate of the delocalization effect on the overall energy transfer in green sulfur bacteria has to await the molecular structure of FMO-reaction center supercomplexes. From 2D experiments on whole cells of *C. tepidum* at 77 K, a time constant of 17 ps has been inferred for the transfer between the FMO protein and the reaction center complex.^{87} At physiological temperatures, the transfer is most likely somewhat faster, but still in the same order of magnitude as the intermonomer transfer occurring on a 10 ps time scale according to the present calculations. Hence, the system might still take advantage of the intermonomer transfer in the FMO protein in order to transfer excitation energy efficiently to the reaction center complex.

## VI. CONCLUSIONS

The main conclusions of the present work are as follows: (i) Intermonomer transfer in the FMO protein occurs with a time constant of about 10 ps at room temperature and may well compete with the transfer between the FMO protein and the reaction center complex. (ii) The low-frequency normal modes of the FMO protein are delocalized over the whole trimer and cause correlations and anticorrelations in the fluctuations of site energies of different pigments in the same and in different FMO monomers. (iii) The effect of these correlations on intra- as well as intermonomer energy transfer is negligible mainly because of the different spectral shapes of the diagonal and the off-diagonal parts of the intermolecular spectral density. (iv) Exciton delocalization in the FMO monomer slows down the intermonomer transfer by about a factor of 2 at physiological temperatures. (v) A direct signature of this exciton delocalization is seen in the wavelength dependence of the lifetime of low-energy exciton states at cryogenic temperatures detected with hole burning spectroscopy and explained in the present work. (vi) Intramonomer exciton relaxation involves the dissipation of intrapigment vibrational quanta, while intermonomer transfer is dominated by the coupling to low-frequency intermolecular vibrations. Any fine details of the spectral density are not so critical for intramonomer exciton relaxation and intermonomer energy transfer. A classical description of nuclear motion works well for the intermonomer transfer at room temperature.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the atomic partial charges of the ground and excited state of BChl *a*.

## ACKNOWLEDGMENTS

Financial support by the Austrian Science Fund (FWF) (Grant No. P 33155-NBL) and by the Austrian Federal Ministry of Science, Research and Economy and the State of Upper Austria (Grant No. LIT-2017-4-SEE-009) is gratefully acknowledged.

## DATA AVAILABILITY

The atomic partial charges of the ground and excited state of BChl *a*, as well as the parameters of the exciton Hamiltonian, the optimized structural coordinates of the FMO trimer, the eigenfrequencies *ω*_{ξ} of the NMA, and the coupling constants *g*_{ξ}(*m*) of the exciton-vibrational coupling can be downloaded from the Zenodo repository via https://doi.org/10.5281/zenodo.4110066.

### APPENDIX A: DERIVATION OF RATE CONSTANT

The rate constant $kMa\u2192NbGF$ is obtained in second-order perturbation theory in the interdomain coupling $V^$,

where $R$ denotes the real part, Tr_{vib} denotes a trace over the vibrational degrees of freedom (DOF), $Weq(Ma)$ is the equilibrium statistical operator of the vibrational DOF in the potential energy surface (PES) of exciton state $Ma$, and *U*_{c}(*τ*) is the time-evolution operator of domain *c*,

with the time evolution operator

containing $Hc(0)$ in Eq. (3) and the *S*-operator of domain *c*,

where $Vc^(\tau )=(Uc(0)(\tau ))\u2020Vc^Uc(0)(\tau )$ contains the time-ordering operator $T$ and *V*_{c} in Eq. (5). In secular approximation, the trace over the vibrational degrees of freedom in Eq. (A1) can be factorized according to Ref. 55,

where $UMc(t)$ is the time evolution operator of the vibrational degrees of freedom in the PES of exciton state $Mc$,

and $Weq(g)$ is the equilibrium statistical operator of the vibrational degrees of freedom in the electronic ground state. The trace in the second line of Eq. (A5) is obtained as

with $\omega MaNb\u2032=\omega Ma\u2032\u2212\omega Nb\u2032$, where these frequencies are defined in Eq. (4), and the time-dependent function

that contains the spectral density

With the help of Eq. (6), the above spectral density can be related to the spectral density of the local exciton-vibrational coupling

In this way, the function $FMaNb(t)$ in Eq. (A8) can be decomposed into

with the function $GKc(t)$ in Eq. (39) that contains the site energy fluctuations and their correlations in domain *c* and the function $GMaNb(t)$, defined in Eq. (22), that contains the correlations in site energy fluctuations of pigments in different domains *a* and *b*. The traces in the first and third line of Eq. (A5) are obtained as^{55}

and

### APPENDIX B: PROTONATION STATES AND SITE ENERGIES

Based on the geometry-optimized structure of the FMO trimer (see Sec. III), site energies of the eight BChl *a* pigments were computed by employing the PBQC method^{4} to test the APCs derived from CAM-B3LYP. The PBQC method was implemented as described in Ref. 89 using dielectric constants of *ϵ*_{p} = 4.0, *ϵ*_{solv} = 80.0 for protonation patterns and $\u03f5\u0303p=1.5$, *ϵ*_{solv} = 80.0 for pigment–protein interactions as well as an ionic strength in the outer medium corresponding to a 0.1M NaCl solution with solvent radius 1.4 Å and ion radius 2.0 Å. In the Monte Carlo titration of protonation states, it was found that His 12, the only histidine considered as titratable, was largely positively charged in the pH range from 8 to 11 (i.e., 100% at pH 8.0; 60% at pH 11.0). In the same pH range, all other titratable groups were found to be in their respective standard protonation states. Consequently, no re-optimization of the FMO trimer (net charge: −6*e*) was necessary to account for non-standard protonation states, while His 12 was modeled as charged in all CHARMM computations (see part III). The resulting site energy shifts (in cm^{−1}: 1, 59; 2, 62; 3, −255; 4,-184; 5, 101; 6, 41; 7, 31; and 8, 174; averaged over the three monomers) are plotted against fitted site energies (“*holo*” from Ref. 7) in Fig. 14, revealing a good correlation and confirming the consistency of our methods employing the CAM-B3LYP functional.

### APPENDIX C: INTERMOLECULAR REORGANIZATION ENERGIES AND EXCITON MATRIX

#### 1. Intermolecular reorganization energies

Table II shows reorganization energies of the intermolecular exciton-vibrational coupling of the BChl pigments in units of cm^{−1}. Table III shows site energies (diagonal) and excitonic couplings (off-diagonal) of BChl pigments in the monomeric subunit (exciton domain) of the FMO protein in units of cm^{-1}.

#### 2. FMO monomer

1 | 2 | 3 | 4 | 5 | 6 | 7 | 8″(24) | |

1 | 12 505 | −94.8 | 5.5 | −5.9 | 7.1 | −15.1 | −12.2 | 39.5 |

2 | 12 425 | 29.8 | 7.6 | 1.6 | 13.1 | 5.7 | 7.9 | |

3 | 12 195 | −58.9 | −1.2 | −9.3 | 3.4 | 1.4 | ||

4 | 12 375 | −64.1 | −17.4 | −62.3 | −1.6 | |||

5 | 12 600 | 89.5 | −4.6 | 4.4 | ||||

6 | 12 515 | 35.1 | −9.1 | |||||

7 | 12 465 | −11.1 | ||||||

8″(24) | 12 700 |

1 | 2 | 3 | 4 | 5 | 6 | 7 | 8″(24) | |

1 | 12 505 | −94.8 | 5.5 | −5.9 | 7.1 | −15.1 | −12.2 | 39.5 |

2 | 12 425 | 29.8 | 7.6 | 1.6 | 13.1 | 5.7 | 7.9 | |

3 | 12 195 | −58.9 | −1.2 | −9.3 | 3.4 | 1.4 | ||

4 | 12 375 | −64.1 | −17.4 | −62.3 | −1.6 | |||

5 | 12 600 | 89.5 | −4.6 | 4.4 | ||||

6 | 12 515 | 35.1 | −9.1 | |||||

7 | 12 465 | −11.1 | ||||||

8″(24) | 12 700 |

#### 3. Intermonomer couplings

Table IV shows excitonic couplings between BChl pigments in two monomeric subunits of the FMO protein in units of cm^{−1}.

1′(9) | 2′(10) | 3′(11) | 4′(12) | 5′(13) | 6′(14) | 7′(15) | 8 (8) | |

1 | 1.7 | 2.4 | 2.1 | 0.4 | 1.0 | 0.0 | 0.5 | 0.2 |

2 | 0.7 | −0.1 | 0.6 | 0.6 | 1.5 | 1.0 | 0.5 | 0.9 |

3 | −0.8 | −3.5 | −3.6 | 0.6 | 1.7 | 1.0 | −0.8 | 2.3 |

4 | 0.3 | −3.3 | 7.6 | 3.0 | −0.5 | 2.2 | 7.2 | −2.3 |

5 | 4.0 | 11.3 | 6.4 | −1.2 | 2.6 | −2.4 | −2.4 | 6.6 |

6 | 2.4 | 7.1 | 2.9 | −0.5 | −0.0 | −2.2 | 0.0 | −3.5 |

7 | 1.3 | 2.8 | 7.3 | 2.9 | −0.9 | 2.7 | 9.6 | −8.6 |

8″(24) | 0.2 | 1.1 | 1.2 | −1.4 | 2.0 | −1.4 | −3.3 | 5.3 |

1′(9) | 2′(10) | 3′(11) | 4′(12) | 5′(13) | 6′(14) | 7′(15) | 8 (8) | |

1 | 1.7 | 2.4 | 2.1 | 0.4 | 1.0 | 0.0 | 0.5 | 0.2 |

2 | 0.7 | −0.1 | 0.6 | 0.6 | 1.5 | 1.0 | 0.5 | 0.9 |

3 | −0.8 | −3.5 | −3.6 | 0.6 | 1.7 | 1.0 | −0.8 | 2.3 |

4 | 0.3 | −3.3 | 7.6 | 3.0 | −0.5 | 2.2 | 7.2 | −2.3 |

5 | 4.0 | 11.3 | 6.4 | −1.2 | 2.6 | −2.4 | −2.4 | 6.6 |

6 | 2.4 | 7.1 | 2.9 | −0.5 | −0.0 | −2.2 | 0.0 | −3.5 |

7 | 1.3 | 2.8 | 7.3 | 2.9 | −0.9 | 2.7 | 9.6 | −8.6 |

8″(24) | 0.2 | 1.1 | 1.2 | −1.4 | 2.0 | −1.4 | −3.3 | 5.3 |