Excited state energy fluctuations in the Fenna–Matthews–Olson complex from molecular dynamics simulations with interpolated chromophore potentials

Chang Woo Kim a, Bongsik Choi b and Young Min Rhee *b
aDepartment of Chemistry, Pohang University of Science and Technology (POSTECH), Pohang 37673, Korea
bDepartment of Chemistry, Korea Advanced Institute of Science and Technology (KAIST), Daejeon 34141, Korea. E-mail: ymrhee@kaist.ac.kr

Received 14th September 2017 , Accepted 21st November 2017

First published on 21st November 2017

We analyze the environment-induced fluctuation of pigment excitation energies in the Fenna–Matthews–Olson (FMO) complex from various perspectives, by employing an interpolation-based all-atom 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 non-equilibrium excited state trajectory simulations, we show that the system lies well within the linear response regime after photo-absorption 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.

1. Introduction

Light-harvesting complexes play a key role in capturing photons and transferring their energies efficiently to relevant reaction sites in the early stage of photosynthesis. The underlying mechanisms of the excitation energy transfer (EET) have been extensively studied over the past few decades, and some recent interest has been heavily focused on the role of quantum effects during EET processes. Theoretical studies on the EET processes in photosynthetic systems require appropriate treatments of environmental perturbations exerted on electronic subsystems. Full quantum chemical calculations for entire pigment–protein complexes are generally not feasible, and accordingly one must rely on approximate models. Widely used strategies for modeling environmental effects can be categorized into either abstract or atomistic approaches. In an abstract approach, one employs analytically determined forms of spectral densities for pigment site energies and inter-pigment couplings, and deduces simplified but essential insights into EET dynamics.1–4 On the other hand, in an atomistic approach one usually tries to explicitly include the nuclear degrees of freedom during the simulation. Practically, this is accomplished by propagating the nuclei by molecular dynamics (MD) and then evaluating the subsystem parameters through quantum chemical calculations.5–8 Such an approach also allows one to further investigate how insights obtained from an abstract approach are actually embodied in the realistic system. A major drawback for this approach is of course the heavy computational burden posed by quantum chemical calculations over a large number of frames recorded during the MD simulations. Due to this limitation, attempts of atomistically simulating EET processes have still been reported relatively rarely. A few studies have also pointed out that undesirable effects can arise from the mismatch between different potential energy surface (PES) models employed for propagating the nuclei and for calculating the subsystem parameters.9–11 We have recently reported that these obstacles can be overcome12 by applying an interpolated mechanics/molecular mechanics (IM/MM) method,13 using the case of bacteriochlorophyll a (BChl a) molecules in the Fenna–Matthews–Olson (FMO) complex of green sulfur bacteria Chlorobaculum tepidum.14 Due to its simplicity, the FMO complex has been one of the most widely used prototype systems for studying EET processes.7,11,15–19 We have demonstrated that our BChl a model can properly describe the topology of S0 and S1 electronic states of all chromophores while keeping the computational cost close to that with the conventional molecular mechanics (MM) model.12 In fact, our PES enables us to feasibly perform simulations well over nanosecond timescales, with a consistency established between the stages of the trajectory propagation and the site energy calculations.9 Namely, adopting our PES overcomes the two major drawbacks mentioned above toward performing atomistic simulations.

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 non-equilibrium 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.

2. Methods

All the simulations in this study were conducted by using the GROMACS 4.0 package24 with our in-house library for PES interpolation.13 We constructed our simulation system by placing an FMO monomer of green sulfur bacteria Chlorobaculum Tepidum (PDB ID: 3BSD)25 in a rhombic dodecahedron box of 70.62 Å side length and by solvating the complex with 21053 TIP3P26 water molecules. The AMBER99sb force field27 was employed for describing the surrounding protein. The interpolated PES for S0 and S1 electronic states of BChl a were taken from our previous work,12 while the atomic charges of the pigments were newly fitted by the combination of QM/MM calculations and the restrained electrostatic potential (RESP)28 method. The detailed fitting procedures for atomic charges are described in the ESI with Fig. S1 and S2. As we pointed out in our earlier work,12 large distortions occasionally occurred for nine covalent bonds that were not included in the valence interpolation coordinates. We improved against this unphysical behavior by attaching a stronger MM stabilizer potential13 for the bonds as defined in the force field from the Hessian-matching process29 without any scaling. This reduced the root-mean-square error of the pigment vibrational energy with respect to reference quantum chemical data (density functional theory (DFT) and time-dependent DFT (TDDFT) with the B3LYP30 functional) from 0.15 eV to 0.05 eV for both S0 and S1 electronic states. The accuracy of the S0−S1 energy gap was nearly unaffected by this modification, with the calculated root-mean-square error of 7 meV. With this setting for the potential, the system was first subjected to 1000 steps of energy minimization using the steepest descent method,31 followed by 1 ns equilibration and subsequent 100 ns production MD runs in the S0 electronic state with 1 fs time step. During these MD runs, the system was kept under isothermal–isobaric conditions of 300 K and 1 bar through the control by the velocity-rescale thermostat32 and the Parrinello–Rahman barostat33 with a time constant of 1 ps. Hydrogen atoms were constrained with the LINCS algorithm.34 The electrostatic and dispersion interactions were treated via force shifting with 12 Å cutoffs together with periodic boundary conditions. We saved the atomic positions and velocities every 1 ps during the production run, which resulted in a total of 105 snapshots evenly spread along the 100 ns trajectory for generating histograms toward obtaining excitation energy distributions. These snapshots also served as the initial conditions for additional simulations in later parts, with which we examine the environment-induced fluctuations in detail by recording frames with much finer time resolution.

Describing excited state dynamics requires reliable assessments of site energies and inter-pigment couplings. These are also essential for converting a site-based description to an exciton-based 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 105 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 non-equilibrium 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 inter-pigment 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 Ei and the site-site coupling Vij thus defined, the time-dependent Frenkel exciton Hamiltonian is given as

image file: c7cp06303b-t1.tif(1)
where R is the nuclear coordinates of the entire system. We should note that the time evolution of nuclear position R in this equation relies on the ground state equilibrium simulation in our approach, and that the system will actually evolve differently in the presence of electronic excitation, as the excitation will definitely affect the force exerted on the nuclei.6

With the constructed multistate potential model, the excited state responses can be derived either from equilibrium simulations or from non-equilibrium 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:

image file: c7cp06303b-t2.tif(2)
where N is the number of trajectories and the bracket with t denotes the usual time-averaging toward calculating the time correlation. Of course, ΔEn(t) = En(t) − [E with combining overline]n is the energy fluctuation measured from [E with combining overline]n given as the average energy over the n-th trajectory. We chose N = 100, and the snapshots along the aforementioned 100 ns trajectory were adopted as the initial seeds of these new MD simulations. As was also mentioned, the site energies were recorded much more frequently in these additional simulations at every 5 fs to capture the detailed effects of pigment vibrations with a good time resolution in τ. Meanwhile, with non-equilibrium simulations on the interpolated excited state PES, the Stokes shift function (SSF) of a pigment site was obtained by ensemble-averaging the time series of excited state energies:
image file: c7cp06303b-t3.tif(3)
where the energy En(τ) was calculated by first vertically exciting one pigment site and then performing non-equilibrium MD simulations in time τ on the surface with the excited pigment in its S1 electronic state.13En(∞) is the long-time limit of the S1 → S0 transition energy from the simulation, which reflects the energy gap after the pigment dissipates its excess vibrational energy gained upon electronic excitation. We note that EET is not considered here, and the excitation is kept localized at the initially excited pigment during these SSF simulations. The initial snapshots for the vertical excitations with these non-equilibrium simulations were originally the same 100 conformations adopted for computing C(τ). However, as the s(τ) function based on this trajectory ensemble was excessively noisy, we increased the ensemble size by a factor of 10 by re-collecting 1000 initial shapshots from the 100 ns equilibrium trajectory with a constant time spacing of 100 ps. Each non-equilibrium trajectory was run for 20 ps and the excitation energy of the pigment of interest was recorded at every 5 fs in accordance with C(τ). En(∞) was practically computed by averaging En(τ) using the data over the second half of the simulation, i.e. from 10 ps to 20 ps. Of course, the relaxation of En(τ) was visually complete well ahead of τ = 10 ps.

3. Statistical analyses with trajectory snapshots

3.1 Site energies and couplings

The distributions of the site energies for the snapshots recorded along the trajectory are shown in Fig. 1. The widths of all distributions are similar to each other and are also comparable to the spread of peak positions, leading to substantial mutual overlaps. All the peaks exhibit nearly Gaussian shapes. Table 1 lists the averages and the widths of the site energy and inter-pigment coupling distributions. Our inter-pigment coupling strengths in this table match very well with the corresponding values from point dipole descriptions reported by Adolphs and Renger (data not shown).19 While this is somewhat expected as our transition charges were uniformly scaled to match the experimentally derived transition dipole size, it still signifies that the inter-pigment alignments are quite intact over the relatively long 100 ns simulation time. We also note that having correct coupling strengths is key to meaningful interconversions between site and exciton representations.19
image file: c7cp06303b-f1.tif
Fig. 1 Pigment site energy distributions.
Table 1 The averages and distribution widthsa of site energies and inter-pigment couplings (in cm−1 units)b
Site 1 2 3 4 5 6 7
a Root-mean-square deviation from the average. b Numbers in parentheses are the distribution widths.
1 12[thin space (1/6-em)]445














2 12[thin space (1/6-em)]520












3 12[thin space (1/6-em)]205










4 12[thin space (1/6-em)]335








5 12[thin space (1/6-em)]490






6 12[thin space (1/6-em)]640




7 12[thin space (1/6-em)]450


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 ensemble-averaged 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 two-dimensional electronic spectroscopic experiments41 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 large-scale 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).

image file: c7cp06303b-f2.tif
Fig. 2 Visualization of the site energy fluctuations at different time scales. (a) Site energy fluctuations in a short time scale without any averaging. (b) Site energy drifts in a long time scale after averaging over time windows with T = 1 ns. Note that the horizontal axis scales are very different in the two panels.

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 time-windowed 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, time-averaging 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)
Now, if we further assume that the memory times (τm) of these two components are well separated from T, namely τm(δs) ≪ Tτm(δl), averaging over the time window with a sampling period ΔT will generate a fluctuation size
image file: c7cp06303b-t4.tif(5)

image file: c7cp06303b-f3.tif
Fig. 3 Fluctuation size of time-windowed averages as a function of (a) the window size T and (b) its inverse square root, image file: c7cp06303b-t7.tif. The fluctuation size is defined as the standard deviation of the time-windowed averages. Averages were taken from the 100[thin space (1/6-em)]000 snapshots collected at every ps along the trajectory of 100 ns duration. Thus, when T = 1 ns for example, there are 100 windowed averages and the fluctuation size is the standard deviation of these 100 average values. In (b), the inset shows the data in the long time limit, with the dashed line marking the window size of T = 1 ns.

Thus, the fluctuation size δT of the time-windowed average will change linearly with image file: c7cp06303b-t5.tif, 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 time-windowed averages versusimage file: c7cp06303b-t6.tif 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 image file: c7cp06303b-t8.tif, 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.

3.2 Exciton energies

In the above part, we considered the site energies of individual pigment molecules. Let us now change our viewpoint toward delocalized exciton states. Diagonalization of the Frenkel exciton Hamiltonian given in eqn (1) yields the energy eigenvalues and the eigenvectors for delocalized exciton states. One important caveat in handling exciton states is the fact that indexing the energy eigenstates always has some level of ambiguity. A reasonable criterion for sorting energy eigenstates is to follow the state characters. Namely, the states at a time point can be ordered such that their characters match best with the state characters at the previous time point.51 This is useful when the environment-induced fluctuations of the Hamiltonian matrix elements are small enough. If this is the case for our system, the excitonic states from MD snapshots would closely resemble the eigenstates of the average Hamiltonian in Table 1, making the sorting of the states straightforward. However, as we have seen from Fig. 2a, the ordering of the site energies becomes completely disrupted very frequently at ambient temperature. Under these conditions, unequivocally sorting the states in time by following the state characters is practically impossible. This situation is in fact similar to a recent theoretical study on an LH2 complex, which showed that the eigenstates of the circularly symmetric electronic Hamiltonian lose their meaning as the thermal fluctuation drives the complex into a regime that cannot be accessed by a perturbative approach.52 Therefore, we will simply label the exciton states by indexing them in the order of increasing energies. Namely, “exciton 1” will denote the eigenstate with the lowest energy.

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.

image file: c7cp06303b-f4.tif
Fig. 4 The distribution of adiabatic exciton energies.

4. Characteristics of fluctuations

4.1 Spectral densities at 300 K

In this section, we extract the spectral density from the time-dependent gap energy fluctuations of each site and exciton state and discuss the effect of exciton delocalization on environmental coupling. The spectral density characterizes the environmental coupling in the frequency domain, and there have been a number of attempts to extract it from atomistic models either based on the energy gap ACF5,10,18,21 or the normal mode analysis.11,50 We follow the ACF-based method here, as it can be straightforwardly achieved with our PES model. After evaluating C(τ) by using eqn (2) together with the parameters and the procedure detailed in Section 2 and the ESI, the spectral density was obtained with
image file: c7cp06303b-t9.tif(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.

image file: c7cp06303b-f5.tif
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.

image file: c7cp06303b-f6.tif
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

image file: c7cp06303b-t10.tif(7)
and the Huang–Rhys (HR) factors
image file: c7cp06303b-t11.tif(8)
from the spectral densities of all sites. A frequency cutoff of ωc = 20 cm−1 was used for calculating HR factors for stable integration in the low frequency limit.7 The calculated reorganization energies and HR factors are listed in Table 2.

Table 2 Site-based reorganization energies and HR factors
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 site-based densities depend only weakly on the site identity but the heights in exciton-based 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 environment-induced 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 aggregates56–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

image file: c7cp06303b-t12.tif(9)
where [E with combining overline]i is the electronic energy of the i-th site and sik is the HR factor of the k-th phonon mode represented by the bath coordinate qik with the frequency ωk, coupled to the i-th site. We do not assume any correlation between individual bath modes. Transformation to exciton basis yields
image file: c7cp06303b-t13.tif(10)
where the coefficients cαi and cαj are the elements of the unitary transform matrix that connects site and exciton bases. By equating this to image file: c7cp06303b-t14.tif, the exciton state analog of eqn (9), the bath coordinate coupled with the exciton state can be expressed as linear combination of qik:
image file: c7cp06303b-t15.tif(11)

Using the norm-conserving nature of this bath transform, we get the exciton HR factor sαk for the mode with the frequency ωk as

image file: c7cp06303b-t16.tif(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

image file: c7cp06303b-t17.tif(13)

Because image file: c7cp06303b-t18.tif is independent of the bath mode index, this becomes the scaling factor that converts site-based spectral density into the exciton-based 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.

4.2 Reproduction of experimental spectral density

In this section, we will try to reproduce the experimental spectral density by adjusting the simulation settings to those resembling the experimentally adopted conditions. We target the latest report on the spectral density of BChl a pigment in the FMO complex, which employed the difference fluorescence line narrowing spectroscopy scheme (ΔFLN)59 after selectively exciting the red edge region of the absorption spectrum at 4.5 K.54 As explained in the above section, at the red edge the effect of exciton delocalization becomes minimized and the excitation is mostly localized at site 3. We mimicked the situation in the following manner. We adopted the same 100 snapshots employed toward obtaining C(τ) again as the initial conditions and performed a 1 ns long cooling simulation at 4.5 K. Because the timescale of dynamic disorder is noticeably shorter than 1 ns (Fig. 3), this should be long enough to cool down fast vibrational components. Except this temperature, we adopted the same simulation parameters as described in Section 2. Then, we inspected the last 100 ps time domain from each of these 100 trajectories and picked 10 that displayed the smallest average gap energies in order to mimic the red edge conditions. These were then continued with MD for another 10 ns. During this extended simulation period, the lowest exciton energy showed almost no drift for all 10 trajectories, signifying that they stayed in the red edge region (Fig. S5, ESI). Finally, the spectral density for the lowest exciton was calculated from these 10 trajectories with 100 ns aggregate duration.

The spectral density thus obtained is shown with the experimental one54 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 well-acknowledged 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 log-normal 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 mode-based approach11 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 energies62 did not better reproduce the experimental reorganization energies than other functionals.11

image file: c7cp06303b-f7.tif
Fig. 7 The calculated spectral density of the lowest exciton state based on trajectories near the absorption red edge (upper panel, black). The log-normal 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).

4.3 Comparison with Stokes shift function

In an earlier section, we noted that FMO displayed Gaussian statistics under ambient conditions even with the intrinsic anharmonic nature of the pigment potential. This is in fact closely related to the background of using ground state MD simulations for generating spectral densities. In its relation, let us also test whether the relaxation in the excited state resides in the linear response regime. This can be achieved by comparing the energy gap autocorrelation function with the Stokes shift function, which are guaranteed to be equivalent in the linear regime by the fluctuation–dissipation theorem.23,63 Although it is widely accepted that these conditions are satisfied for most pigment–protein complexes, numerical validation has not been feasible with QM/MM approaches. With IM/MM, indeed, such a task becomes quite feasible. For this, we computed SSF (eqn (3)) of each pigment site by recording the site energy fluctuations from direct non-equilibrium MD simulations.

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 well-resolved 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.

image file: c7cp06303b-f8.tif
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 S0 and S1 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.

Table 3 Comparison between the Stokes shift amounts and the reorganization energies (in cm−1 units)
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

5. Summary and conclusions

In this paper, we characterized the environment-induced energetic fluctuations of BChl a pigments in the FMO complex by using an atomistic model constructed by IM/MM. Our IM/MM model is as reliable as the reference DFT calculations12 and is also very efficient.65 With the model, we could conduct sufficiently long MD simulations with consistent descriptions during trajectory propagations and subsequent energy gap evaluations.9 In addition, based on the energy fluctuations in different timescales, we could directly observe the timescale separation between the dynamic and the static disorder. We also probed the behaviors of delocalized exciton states and compared them to the behaviors of the localized site states, based on their spectral densities and related quantities. Through this, we observed that the exciton state was effectively less perturbed with a larger degree of delocalization. Such a behavior could be readily explained with a formulation regarding chromophore aggregates.56–58 We also attempted to reproduce the experimentally reported spectral density by directly following experimental conditions. Finally, by comparing the gap energy autocorrelation function from the ground state equilibrium simulations and the time-dependent Stokes shift from the excited state non-equilibrium simulations, we demonstrated that the photodynamics of the FMO complex at ambient temperature are well within the linear response regime.

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 “bottom-up” 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 polarizations66–68 will help us to further develop our modeling capability toward reaching that level.

Conflicts of interest

There are no conflicts to declare.


This work was supported by the Institute for Basic Science (IBS), the National Research Foundation of Korea (Grant 2017R1A2B3004946), and Samsung Electronics. We thank Prof. Seogjoo Jang (City University of New York) for discussions on spectral densities, and Dr Hyun Woo Kim (Korea Research Institute of Chemical Technology) for discussions on exciton delocalization effects.


  1. A. Ishizaki, T. R. Calhoun, G. S. Schlau-Cohen and G. R. Fleming, Phys. Chem. Chem. Phys., 2010, 12, 7319–7337 RSC .
  2. C. Kreisbeck and T. Kramer, J. Phys. Chem. Lett., 2012, 3, 2828–2833 CrossRef CAS .
  3. N. Christensson, H. F. Kauffmann, T. Pullerits and T. Mančal, J. Phys. Chem. B, 2012, 116, 7449–7454 CrossRef CAS PubMed .
  4. E. Romero, R. Augulis, V. I. Novoderezhkin, M. Ferretti, J. Thieme, D. Zigmantas and R. van Grondelle, Nat. Phys., 2014, 10, 676–682 CrossRef CAS PubMed .
  5. C. Olbrich, J. Strümpfer, K. Schulten and U. Kleinekathöfer, J. Phys. Chem. Lett., 2011, 2, 1771–1776 CrossRef CAS PubMed .
  6. H. W. Kim, A. Kelly, J. W. Park and Y. M. Rhee, J. Am. Chem. Soc., 2012, 134, 11640–11651 CrossRef CAS PubMed .
  7. X. Wang, G. Ritschel, S. Wuster and A. Eisfeld, Phys. Chem. Chem. Phys., 2015, 17, 25629–25641 RSC .
  8. C. Curutchet and B. Mennucci, Chem. Rev., 2016, 117, 294–343 CrossRef PubMed .
  9. C. W. Kim, J. W. Park and Y. M. Rhee, J. Phys. Chem. Lett., 2015, 6, 2875–2880 CrossRef CAS PubMed .
  10. A. M. Rosnik and C. Curutchet, J. Chem. Theory Comput., 2015, 11, 5826–5837 CrossRef CAS PubMed .
  11. M. K. Lee and D. F. Coker, J. Phys. Chem. Lett., 2016, 7, 3171–3178 CrossRef CAS PubMed .
  12. C. W. Kim and Y. M. Rhee, J. Chem. Theory Comput., 2016, 12, 5235–5246 CrossRef CAS PubMed .
  13. J. W. Park, H. W. Kim, C.-I. Song and Y. M. Rhee, J. Chem. Phys., 2011, 135, 014107 CrossRef PubMed .
  14. R. E. Fenna and B. W. Matthews, Nature, 1975, 258, 573–577 CrossRef CAS .
  15. T. Brixner, J. Stenger, H. M. Vaswani, M. Cho, R. E. Blankenship and G. R. Fleming, Nature, 2005, 434, 625–628 CrossRef CAS PubMed .
  16. J. Adolphs and T. Renger, Biophys. J., 2006, 91, 2778–2797 CrossRef CAS PubMed .
  17. P. Huo and D. F. Coker, J. Chem. Phys., 2011, 135, 201101 CrossRef PubMed .
  18. E. Rivera, D. Montemayor, M. Masia and D. F. Coker, J. Phys. Chem. B, 2013, 117, 5510–5521 CrossRef CAS PubMed .
  19. C. G. Gillis and G. A. Jones, J. Phys. Chem. B, 2015, 119, 4165–4174 CrossRef CAS PubMed .
  20. C. Olbrich, T. L. C. Jansen, J. Liebers, M. Aghtar, J. Strümpfer, K. Schulten, J. Knoester and U. Kleinekathöfer, J. Phys. Chem. B, 2011, 115, 8609–8621 CrossRef CAS PubMed .
  21. S. Shim, P. Rebentrost, S. Valleau and A. Aspuru-Guzik, Biophys. J., 2012, 102, 649–660 CrossRef CAS PubMed .
  22. M. Higashi and S. Saito, J. Chem. Theory Comput., 2016, 12, 4128–4137 CrossRef CAS PubMed .
  23. D. Chandler, Introduction to Modern Statistical Mechanics, Oxford University Press, New York, 1987 Search PubMed .
  24. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447 CrossRef CAS PubMed .
  25. A. Ben-Shem, F. Frolow and N. Nelson, FEBS Lett., 2004, 564, 274–280 CrossRef CAS PubMed .
  26. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935 CrossRef CAS .
  27. V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg and C. Simmerling, Proteins: Struct., Funct., Bioinf., 2006, 65, 712–725 CrossRef CAS PubMed .
  28. C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, J. Phys. Chem., 1993, 97, 10269–10280 CrossRef CAS .
  29. C.-I. Song and Y. M. Rhee, Int. J. Quantum Chem., 2011, 111, 4091–4105 CrossRef CAS .
  30. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789 CrossRef CAS .
  31. P. M. Morse and H. Feshbach, Methods of Theoretical Physics, Part I, McGraw-Hill, New York, 1953 Search PubMed .
  32. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101 CrossRef PubMed .
  33. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190 CrossRef CAS .
  34. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472 CrossRef CAS .
  35. S. Jurinovich, C. Curutchet and B. Mennucci, ChemPhysChem, 2014, 15, 3194–3204 CrossRef CAS PubMed .
  36. S. Jurinovich, L. Viani, I. G. Prandi, T. Renger and B. Mennucci, Phys. Chem. Chem. Phys., 2015, 17, 14405–14416 RSC .
  37. L. Cupellini, S. Jurinovich, M. Campetella, S. Caprasecca, C. A. Guido, S. M. Kelly, A. T. Gardiner, R. Cogdell and B. Mennucci, J. Phys. Chem. B, 2016, 120, 11348–11359 CrossRef CAS PubMed .
  38. M. E. Madjet, A. Abdurahman and T. Renger, J. Phys. Chem. B, 2006, 110, 17268–17281 CrossRef CAS PubMed .
  39. R. S. Knox and B. Q. Spring, Photochem. Photobiol., 2003, 77, 497–501 CrossRef CAS PubMed .
  40. C. Olbrich and U. Kleinekathöfer, J. Phys. Chem. B, 2010, 114, 12427–12437 CrossRef CAS PubMed .
  41. D. M. Jonas, Annu. Rev. Phys. Chem., 2003, 54, 425–463 CrossRef CAS PubMed .
  42. H. P. Lu, L. Xun and X. S. Xie, Science, 1998, 282, 1877–1882 CrossRef CAS PubMed .
  43. I. Barvík, C. Warns, T. Neidlinger and P. Reineker, Chem. Phys., 1999, 240, 173–189 CrossRef .
  44. M. Fuxreiter, Mol. Biosyst., 2012, 8, 168–177 RSC .
  45. K. Henzler-Wildman and D. Kern, Nature, 2007, 450, 964–972 CrossRef CAS PubMed .
  46. H. Lee, Y.-C. Cheng and G. R. Fleming, Science, 2007, 316, 1462–1465 CrossRef CAS PubMed .
  47. V. Tiwari, W. K. Peters and D. M. Jonas, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 1203–1208 CrossRef CAS PubMed .
  48. V. Butkus, D. Zigmantas, D. Abramavicius and L. Valkunas, Chem. Phys. Lett., 2013, 587, 93–98 CrossRef CAS .
  49. V. Butkus, H. Dong, G. R. Fleming, D. Abramavicius and L. Valkunas, J. Phys. Chem. Lett., 2016, 7, 277–282 CrossRef CAS PubMed .
  50. T. Renger, A. Klinger, F. Steinecker, M. Schmidt am Busch, J. Numata and F. Müh, J. Phys. Chem. B, 2012, 116, 14565–14580 CrossRef CAS PubMed .
  51. Y. M. Rhee, D. Casanova and M. Head-Gordon, J. Phys. Chem. A, 2009, 113, 10564–10576 CrossRef CAS PubMed .
  52. C. Stross, M. W. Van der Kamp, T. A. A. Oliver, J. N. Harvey, N. Linden and F. R. Manby, J. Phys. Chem. B, 2016, 120, 11449–11463 CrossRef CAS PubMed .
  53. S. Valleau, A. Eisfeld and A. Aspuru-Guzik, J. Chem. Phys., 2012, 137, 224103 CrossRef PubMed .
  54. M. Rätsep and A. Freiberg, J. Lumin., 2007, 127, 251–259 CrossRef .
  55. D. Abramavicius and L. Valkunas, Photosynth. Res., 2016, 127, 33–47 CrossRef CAS PubMed .
  56. R. M. Hochstrasser and P. N. Prasad, J. Chem. Phys., 1972, 56, 2814–2823 CrossRef CAS .
  57. D. M. Hanson, Chem. Phys. Lett., 1976, 43, 217–220 CrossRef CAS .
  58. J. P. Lemaistre, J. Lumin., 2004, 107, 332–338 CrossRef CAS .
  59. M. Rätsep and A. Freiberg, Chem. Phys. Lett., 2003, 377, 371–376 CrossRef .
  60. A. P. Scott and L. Radom, J. Phys. Chem., 1996, 100, 16502–16513 CrossRef CAS .
  61. A. Kell, X. Feng, M. Reppert and R. Jankowiak, J. Phys. Chem. B, 2013, 117, 7317–7323 CrossRef CAS PubMed .
  62. N. H. List, C. Curutchet, S. Knecht, B. Mennucci and J. Kongsted, J. Chem. Theory Comput., 2013, 9, 4928–4938 CrossRef CAS PubMed .
  63. T. Joo, Y. Jia, J. Y. Yu, M. J. Lang and G. R. Fleming, J. Chem. Phys., 1996, 104, 6089–6108 CrossRef CAS .
  64. If ACF is obtained with equilibrium simulations on the excited state PES, the underdamped vibrational structures will agree better with the results from SSF.
  65. J. W. Park and Y. M. Rhee, ChemPhysChem, 2014, 15, 3183–3193 CrossRef CAS PubMed .
  66. C. Curutchet, A. Muñoz-Losa, S. Monti, J. Kongsted, G. D. Scholes and B. Mennucci, J. Chem. Theory Comput., 2009, 5, 1838–1848 CrossRef CAS PubMed .
  67. J. A. Lemkul, J. Huang, B. Roux and A. D. MacKerell, Chem. Rev., 2016, 116, 4983–5013 CrossRef CAS PubMed .
  68. A. Albaugh, A. M. N. Niklasson and T. Head-Gordon, J. Phys. Chem. Lett., 2017, 8, 1714–1723 CrossRef CAS PubMed .


Electronic supplementary information (ESI) available. See DOI: 10.1039/c7cp06303b

This journal is © the Owner Societies 2018