Chang Woo
Kim
^{a},
Bongsik
Choi
^{b} and
Young Min
Rhee
*^{b}
^{a}Department of Chemistry, Pohang University of Science and Technology (POSTECH), Pohang 37673, Korea
^{b}Department of Chemistry, Korea Advanced Institute of Science and Technology (KAIST), Daejeon 34141, Korea. Email: ymrhee@kaist.ac.kr
First published on 21st November 2017
We analyze the environmentinduced fluctuation of pigment excitation energies in the Fenna–Matthews–Olson (FMO) complex from various perspectives, by employing an interpolationbased allatom potential energy model for describing realistic pigment vibrations. We conduct molecular dynamics simulations on a 100 ns timescale, which is an extent that can enclose the effect of static disorder, and demonstrate its timescale separation from fast dynamic disorder. We extract the spectral densities of the complex by considering both the site and the exciton bases. We show that exciton delocalization reduces the effective environmental fluctuation and rationalize this aspect based on a model of fluctuating molecular aggregates. We also obtained the spectral density of the lowest exciton state under low temperature conditions and show that it reasonably well reproduces the experimental result. Finally, by additionally performing nonequilibrium excited state trajectory simulations, we show that the system lies well within the linear response regime after photoabsorption and that the pigments do not visit anharmonic regions of the potential surface to a significant extent. This indicates that methodologies based on harmonic bath models are indeed reasonable approaches for describing the excited state dynamics of the FMO complex.
With such a reliable and efficient PES available, it will be worthwhile characterizing and visualizing the nature of environmental perturbations engraved on the surface. This study features the results from a series of molecular dynamics simulations based on our FMO model, with a primary focus on observing how environmental fluctuations affect the energies of the involved states. For such observations, a spectral density has a good utility as it shows the distribution of the fluctuating modes and their coupling strengths. In addition to the individual site based consideration, which has been already reported in a number of earlier studies,^{5,18,20–22} we extend to analyzing environmental perturbations on delocalized exciton states and discuss how the basis conversion affects the analysis results. We then attempt to reproduce the experimental spectral density in a situation that resembles the experimental conditions at low temperature. We also show that the results from nonequilibrium simulations on the excited state surface correspond quite well with the ones from time correlations obtained from ground state equilibrium simulations. Namely, we show that the system is well within the linear response regime.^{23} This suggests that the spectral density obtained from equilibrium simulations can be safely used for characterizing excited state dynamics of the FMO complex.
This article is organized as follows: Section 2 presents the details of our model construction and the specifics of our trajectory simulations. In Section 3, statistical analyses are performed for the snapshots along the trajectories by adopting both site and exciton bases. Section 4 is dedicated to characterizing environmental perturbations through the use of spectral densities. These are followed by some concluding remarks in Section 5.
Describing excited state dynamics requires reliable assessments of site energies and interpigment couplings. These are also essential for converting a sitebased description to an excitonbased one. For site energies, accurately reproducing experimental values is not trivial and often requires QM/MM calculations with a polarizable MM environment with large QM regions.^{22,35–37} As we were again adopting a rather simple point charge model from our earlier work,^{12} we decided to apply a posteriori calibrations to the pigment site energies. This was practically attained by applying a constant energy shift to each individual pigment such that the average site energy of each pigment over all 10^{5} snapshots matched the corresponding experimental value. We note that the experimental averages themselves were also obtained by similarly adjusting the site energies simultaneously so that good fits were attained against experimentally determined absorption and linear/circular dichroism spectra.^{16} The calibrated site energies were employed for both generating gap energy autocorrelation and performing excited state nonequilibrium simulations as detailed in the next paragraph. Philosophically, our nonpolarizable approach is based on an assumption that the pigment vibration is insensitive to the pigment polarization as it affects the average of the site energy more than its fluctuation. Namely, polarization will mostly shift the average site energy but will affect the excited state energy fluctuations to a limited extent. For interpigment couplings, the method of transition charges from electrostatic potential (TrESP)^{38} was used. This evaluates the electronic coupling between pigments as interaction energy between point transition charges placed at atomic positions of the pigments. The TrESP point charges were obtained as described in the ESI,† and were uniformly scaled by a factor of 0.86 to match the experimentally determined transition dipole moment of 6.1 Debye.^{39} With the pigment site energy E_{i} and the sitesite coupling V_{ij} thus defined, the timedependent Frenkel exciton Hamiltonian is given as
(1) 
With the constructed multistate potential model, the excited state responses can be derived either from equilibrium simulations or from nonequilibrium simulations. In the equilibrium simulation case, we included the possible variations arising from static disorder by averaging the autocorrelation functions (ACFs) from a batch of trajectories. This aspect is in contrast to earlier studies based on single trajectories,^{5,18,40} and exemplifies the speedy nature of our interpolation scheme. More specifically, we calculated the ACF from the time series of excitation energy by employing multiple trajectories:
(2) 
(3) 
Site  1  2  3  4  5  6  7 

a Rootmeansquare deviation from the average. b Numbers in parentheses are the distribution widths.  
1 
12445
(260) 
−112.1
(12.7) 
5.0
(1.3) 
−9.9
(1.6) 
10.3
(2.0) 
−28.9
(8.4) 
−13.7
(3.5) 
2 
12520
(253) 
48.2
(2.5) 
7.4
(1.3) 
5.3
(2.1) 
9.5
(2.3) 
3.5
(4.4) 

3 
12205
(278) 
−81.4
(11.1) 
−2.5
(1.6) 
−10.5
(0.5) 
−5.8
(9.2) 

4 
12335
(240) 
−86.0
(10.8) 
−18.9
(2.1) 
−72.7
(9.5) 

5 
12490
(243) 
86.9
(15.9) 
0.7
(7.5) 

6 
12640
(251) 
52.5
(10.4) 

7 
12450
(257) 
The widths of the site energy distributions reflect the extent of site energy fluctuations induced by environmental disorder. In the physical process involving macroscopic systems, the overall disorder is often disentangled into fast dynamic and slow static components. Each type of disorder plays a different role in ensembleaveraged dynamics observed by experiments. Namely, dynamic disorder governs the behavior of a single ensemble element while static disorder dictates the differences between different ensemble elements. Strictly speaking, therefore, the distinction relies on the type of involved experiments. The timescale of EET processes in the FMO complex is on the order of a few picoseconds,^{15} and the twodimensional electronic spectroscopic experiments^{41} on this system have mainly targeted up to around this time regime. For protein systems, especially with simulations, a somewhat loose criterion based on atomic motions is also often used for convenience. With this, dynamic disorder refers to fast fluctuations induced by localized thermal motions such as bond vibrations, while static disorder refers to slow energy drifts induced by the largescale motions of the surrounding protein.^{42–44} According to this classification, the timescale of static disorder is considered to be in the order of nanoseconds.^{45} At any rate, the timescale of static disorder is much longer than the timescale of exciton transfer within a single complex, and is also beyond the reach of many sophisticated nonadiabatic dynamics simulation tools. For this reason, studies based on abstract Hamiltonian models have sometimes taken the static disorder into account by assuming an ensemble of Hamiltonians with artificial random distributions of site energies.^{3,46–49} In our case with the atomistic model, such an assumption is not necessary as simulating for a sufficiently long period is quite feasible and can naturally include both dynamic and static disorder. Indeed, it will be interesting to examine how the two types of disorder manifest themselves during the simulation.
We can trivially infer the magnitude of the dynamic disorder of any given pigment by observing the fluctuating pattern of its site energy. Fig. 2a presents such patterns for the seven pigments in the FMO complex. One can easily notice that the dynamic disorder is large enough to completely disrupt the energy ordering besides site 3, which spends most of the time at the bottom of the energy manifold. Although strictly separating the static disorder components from the fluctuating patterns is in principle impossible, one can still apply a heuristic approach of estimating them. Let us suppose plotting the site energies after averaging over time windows of duration T. For example, Fig. 2b displays the average site energies in every ns time window (namely, T = 1 ns). Now, this figure only displays fluctuations that occur beyond the 1 ns timescale. One can see that the energy drift in this case is much smaller than the one shown in Fig. 2a, and the averaged site energies are generally well ordered in the long timescale (t > 1 ns).
What is actually more intriguing is the fact that we could observe clear timescale separation by scanning T. When we measured the fluctuation size of the timewindowed averages as a function of T for each pigment, we obtained a sharply changing pattern as displayed in Fig. 3a. For all pigment sites, up to ∼0.2 ns, increasing T feasibly averaged out fast fluctuations (dynamic disorder). Beyond that timescale, however, timeaveraging does not wash out site energy fluctuations that well even at T = 10 ns. This behavior strongly indicates that there are at least two different types of fluctuations that have markedly different temporal behaviors. For a quantitative analysis, let us assume that the site energy changes in time with two sources: long time fluctuation (δ_{l}) and short time fluctuation (δ_{s}). Then, the total fluctuation is simply given as
δ = δ_{l} + δ_{s}  (4) 
(5) 
Thus, the fluctuation size δ_{T} of the timewindowed average will change linearly with , and its extrapolation to the T → ∞ limit, namely δ_{l}, will suggest the size of the static disorder. In addition, because the assumption T ≪ τ_{m} (δ_{l}) will eventually break down as T approaches the timescale of static disorder, δ_{T} will also start to deviate from linearity represented by eqn (5) at large enough T. Thus, by observing the deviating point, we can also estimate the timescale of static disorder. When we plotted the fluctuation size of the timewindowed averages versus for each pigment (Fig. 3b), we indeed observed good linearity in a wide range. This suggests that our distinction between dynamic and static disorder is quite physical.
In addition, by extrapolating the linear portions toward , we can infer that the sizes of static disorder for all pigments are in the 20–50 cm^{−1} range. This corresponds favorably to the typically adopted parameters for simulating exciton dynamics in earlier studies.^{2,50} Also, from the point where the linearity breaks down (T ∼ 1 ns), we can have a direct sense of the timescale of the static disorder.
The distributions of exciton energies thus obtained are plotted in Fig. 4. Again, they form a series of heavily overlapping peaks. By comparing with the site energy distributions shown in Fig. 1, one can notice that the exciton energy distributions in the middle of the manifold become narrower than the others. It can be inferred that the degree of fluctuation for exciton states depends on the position in the manifold, being small at the middle and large at the edge. We postpone the detailed analysis of such a behavior until the next section, where we discuss the effect of exciton delocalization on environmental coupling.
(6) 
Practically, this was achieved with a Fourier transform of ACF in the 0 ≤ τ ≤ 20 ps range. Also, the harmonic prefactor βℏω/2 was used instead of tanh(βℏω/2) to account for the fact that classical MD simulations were used to propagate the system.^{53} According to eqn (6), the relation between ACF and the spectral density is linear. Thus, obtaining J(ω) with the average ACF based on multiple trajectories described in eqn (2) is equivalent to obtaining multiple J(ω) first from each trajectory and then averaging later on. Namely, the same result will be obtained even if we exchange the order of the sum and the integral in eqn (2) and (6). The spirit of employing multiple trajectories toward obtaining an average spectral density is indeed similar to what was adopted in ref. 7. In cases where multiple trajectories at different temperatures should be used, J(ω) must first be calculated from individual trajectories and then averaged as was done in the earlier work.
Fig. 5 displays the complete set of site spectral densities calculated from ACF, with comparisons to the results from some earlier studies. Because the spectral density in an experiment is obtained in the J(ω)/ω^{2} form rather than J(ω) itself,^{55} we also displayed the density in the scaled form in Fig. S3 (ESI†). The spectral densities of the exciton states could also be evaluated in the same manner and are displayed in Fig. S4 (ESI†). Among these, the densities of site 3, exciton 1, and exciton 4 were picked as representative cases and displayed in Fig. 6. Excitons 1 and 4 were selected in order to compare the characteristics of exciton bands at the edge and in the middle, while site 3 with the lowest site energy was selected as its contribution to exciton 1 was the largest.
Fig. 5 The calculated spectral densities J(ω) for all seven pigment sites from 300 K simulations (in colors). Spectral densities calculated by averaging monomer spectral densities from ref. 11 (in gray) are also shown for comparison. For site 3, we also display the experimentally obtained spectral density (in black). The vibrational part of the experimental spectral density was reconstructed up to 1550 cm^{−1} by Gaussian broadening of the Huang–Rhys factors from ref. 54, with a standard deviation of 7 cm^{−1}. 
Fig. 6 The calculated scaled spectral densities J(ω)/ω^{2} for site 3 (black), exciton 1 (red) and exciton 4 (blue) from 300 K simulations. Baselines were vertically shifted for visual clarity. 
Let us first examine the spectral density of site 3 and discuss its general feature. The spectral density can be roughly divided into two parts: a continuous band in the low frequency region and discrete, sharp peaks at high frequencies. The former arises from the combination of intermolecular interactions and valence vibrations involving large scale motions of the pigment, while the latter consists of the remaining vibrational modes that mostly originate from local stretch and bending motions of the pigment.^{9} The heterogeneity of a local protein environment usually induces only slight differences in pigment vibrations, and the spectral densities of other pigments are highly similar in shape to the one from site 3 (Fig. S3, ESI†). We also evaluated the pigment reorganization energies
(7) 
(8) 
Site  1  2  3  4  5  6  7 

λ (cm^{−1})  147  145  168  131  134  142  147 
S  0.46  0.46  0.64  0.31  0.35  0.41  0.45 
When we consider the three spectral densities in Fig. 6 altogether, we can easily see that the spectral peak positions are well overlapping but their relative heights are quite different. In fact, the peak heights are generally largest with site 3 and smallest with exciton 4. In fact, if we examine all site and exciton spectral densities (Fig. S3 and S4, ESI†), the peak heights in sitebased densities depend only weakly on the site identity but the heights in excitonbased densities vary noticeably depending on the exciton level. Interestingly, the overall shape of an exciton spectral density appears to be a reduction of a site spectral density with a constant scaling factor over the entire region. The magnitude of spectral density is proportional to the strength of environmentinduced perturbation, and in Section 3.2 we already demonstrated that the fluctuation of the exciton energy was larger at the edge of the manifold than at the middle (Fig. 4). This explains the general trend observed for relative intensities of exciton spectral densities, yet only in a phenomenological way. For more quantitative analysis, we adapt the theoretical treatment of periodic molecular aggregates^{56–58} which can also be applied to a coupled chromophore system considered here.
For this, let us write the site energy with respect to environmental bath coordinates as
(9) 
(10) 
(11) 
Using the normconserving nature of this bath transform, we get the exciton HR factor s_{αk} for the mode with the frequency ω_{k} as
(12) 
If we assume that the site spectral densities are nearly identical regardless of the pigment site index i as was indeed observed in Fig. S3 (ESI†), we have
(13) 
Because is independent of the bath mode index, this becomes the scaling factor that converts sitebased spectral density into the excitonbased one. Moreover, this factor is less than unity and becomes smaller as the extent of exciton delocalization becomes larger. Then, Fig. 4 and 6 imply that the exciton states in the middle of the energy ladder are more delocalized than the ones at the edge. The lack of delocalization in the edge regions can be explained by the fact that sites 3 and 6, which are respectively the lowest and the highest in energy, form their own exciton states with reduced contributions by other sites.
The spectral density thus obtained is shown with the experimental one^{54} in Fig. 7. The overall shapes are quite similar to each other, although some disagreements in the details are also visible. In particular, a group of experimental peaks around ∼1200 cm^{−1} appear on the bluer side in our result and the peak heights are generally smaller. The discrepancy of the peak positions arises from the wellacknowledged overestimation of experimental vibrational frequency using quantum chemical methods.^{60} The discrepancy in the heights is somewhat more substantial and is reflected in the reorganization energy (λ) of 118 cm^{−1}, in comparison to an experimental value of 220 cm^{−1}. To further analyze the source of this discrepancy, we split the spectral density into phonon and vibrational parts by extracting the phonon band from the calculated spectral density. The continuous, low frequency phonon band was fitted with a lognormal distribution for this purpose and is additionally displayed in Fig. 7. This band gave λ = 21.2 cm^{−1} and S = 0.33, which are slightly larger, not smaller, than the experimental values of λ = 14.6 cm^{−1} and S = 0.29.^{54,61} Therefore, we can infer that the underestimation of environmental perturbation by our model originates from valence vibrations. A recent work which extracted spectral densities of the FMO complex with a normal modebased approach^{11} showed that the B3LYP functional, from which our model is constructed, yields reorganization energies of 140–160 cm^{−1} for the seven pigment sites embedded in the protein. As this agrees very well with our site reorganization energies presented in Table 2, we believe that the underestimation of the intramolecular spectral density components arises from the choice of DFT functional that we engaged toward generating the IM/MM PES. In a sense, this is a limitation that we have to take when applying DFT, because the accuracy level of 100 cm^{−1} is difficult to attain with presently available functionals. It has also been shown that the functional that yielded the most reliable site excitation energies^{62} did not better reproduce the experimental reorganization energies than other functionals.^{11}
Fig. 7 The calculated spectral density of the lowest exciton state based on trajectories near the absorption red edge (upper panel, black). The lognormal fit of the continuous low frequency phonon band is overlapped (upper panel, red) with the density. For comparison, the experimental spectral density from ref. 54 is also displayed (lower panel). 
Fig. 8a displays the comparison for site 1. The detailed shapes of the oscillatory components also appear quite similar until 0.6 ps, with some deviations beyond that time point. The oscillations mainly reflect underdamped intramolecular pigment vibrations and are strongly affected by the PES on which the nuclear motions evolve. Therefore, the difference between the shapes of the ground and excited state PES (e.g. Duschinsky rotation) may introduce some degree of inconsistency in the oscillatory structures of the two curves. The oscillatory behaviors will be better represented with Fourier transforms, shown in Fig. 8b. From this figure, we can see that the discrepancies are likely originating from the inherent lack of sampling in SSF and the overall features generally agree quite well with matching peak positions. As the lengths of the MD trajectories used to obtain C(τ) and s(τ) in Fig. 8a are on the same order, this highlights the superior statistical nature of ACF and the improvement is simply associated with the fact that a snapshot on the equilibrium trajectory is used for evaluating multiple data points in ACF. Nonetheless for this overall agreement, at least one change in the normal mode behaviors can be observed. The sharp and wellresolved peak at 1745 cm^{−1} from ACF is shifted to 1730 cm^{−1} with SSF as marked with a dashed line in Fig. 8b. The oscillations in ACF and SSF will respectively reflect the topologies of the ground and the excited state PES's.^{64} Upon close inspection, we deduced that the mode involves the methyl ester carbonyl stretch. Therefore, the smaller carbonyl stretch frequency in the SSF can be interpreted as slight weakening of this carbonyl bond upon electronic excitation.
Fig. 8 (a) Comparison between the gap energy autocorrelation function (ACF) and the Stokes shift function (SSF) for site 1. (b) Fourier transforms of the same data. 
In the above, we have suggested that the normal modes of pigment vibrations are generally agreeing in the ground and excited states. To further check the consistency of the harmonicity of the system, we can also compare the half of the total Stokes shift, s(0)/2, with the reorganization energy calculated from the spectral density (namely, eqn (3)versuseqn (7)). The two quantities will become identical when the ground and the excited state PES's are both perfectly harmonic and overlapping in their shapes. Actually, this is an assumption that is usually made for analyzing many photosynthetic systems with simplified models. Table 3 presents the reorganization energies calculated from the two separate methods for all seven pigment sites. The agreements are generally quite good for all pigments, and we can infer that the harmonic assumption is indeed well satisfied. Of course, due to the inherent anharmonicities and morphological differences between S_{0} and S_{1} surfaces, the s(0)/2 values are somewhat larger than the λ values derived from the spectral densities. In general, however, we can infer that the BChl a molecules in FMO do not undergo large changes in their charge distributions or structures upon electronic excitation.
Site  1  2  3  4  5  6  7 

s(0)/2  153  166  191  141  153  150  163 
λ  147  145  168  131  134  142  147 
We stress that access to an efficient and reliable model such as our IM/MM potential can be crucial in characterizing and understanding diverse features of an energy transfer system. At the same time, we note that the use of empirical site energy shifts and transition dipole moment scaling prevents our study from being a complete “bottomup” approach that is free from any experimental parameters. Indeed, accurate assessments of such quantities are crucial for truly ab initio description of excited state dynamics, as they directly affect the system Hamiltonian that governs the collective behaviors from all ensemble elements. We anticipate that recent advancements in more reliable quantum chemical and molecular mechanical methods with appropriately accounted mutual polarizations^{66–68} will help us to further develop our modeling capability toward reaching that level.
Footnote 
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7cp06303b 
This journal is © the Owner Societies 2018 