Modeling Electronic-Nuclear Interactions for Excitation Energy

Jul 29, 2016 - Alexander S. Fokas , Daniel J. Cole , Nicholas D. M. Hine , Stephen A. Wells , Michael C. Payne , and Alex W. Chin. The Journal of Phys...
0 downloads 0 Views 1MB Size
Download PDF
Letter pubs.acs.org/JPCL

Modeling Electronic-Nuclear Interactions for Excitation Energy Transfer Processes in Light-Harvesting Complexes Mi Kyung Lee* and David F. Coker* Department of Chemistry, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, United States S Supporting Information *

ABSTRACT: An accurate approach for computing intermolecular and intrachromophore contributions to spectral densities to describe the electronic−nuclear interactions relevant for modeling excitation energy transfer processes in light harvesting systems is presented. The approach is based on molecular dynamics (MD) calculations of classical correlation functions of long-range contributions to excitation energy fluctuations and a separate harmonic analysis and single-point gradient quantum calculations for electron−intrachromophore vibrational couplings. A simple model is also presented that enables detailed analysis of the shortcomings of standard MD-based excitation energy fluctuation correlation function approaches. The method introduced here avoids these problems, and its reliability is demonstrated in accurate predictions for bacteriochlorophyll molecules in the Fenna−Matthews−Olson pigment−protein complex, where excellent agreement with experimental spectral densities is found. This efficient approach can provide instantaneous spectral densities for treating the influence of fluctuations in environmental dissipation on fast electronic relaxation.

L

ight harvesting pigment−protein complexes (PPC)1−4 of plants, bacteria, and algae boast near-unit efficiencies for transferring absorbed photon energy through molecular excitation energy transfer (EET) to reaction centers where energy transduction, through charge separation, is initiated. Regardless of how the PPC chromophore network is excited, e.g., using a coherent laser, or with incoherent sunlight, and how the energy is transferred, e.g., coherently or incoherently, the interaction between the fast electronic and the slower nuclear degrees of freedom plays a crucial role not only in determining the energy transfer rates, but also the pathways to minimize energy loss. EET dynamics not only depend on the network’s electronic eigenstate (exciton) energies and their extent of delocalization, but it also significantly depends on the coupling between the electronic and phonon-like protein motions for electronic energy dissipation, and on the coupling between the electronic and intramolecular (or intrachromophore) vibrations to describe quantum vibronic states. For systems with closely spaced exciton states, such as the Fenna− Matthews−Olson (FMO) complex, the lower frequency intermolecular environmental modes are resonant with many of the excitonic energy gaps and can assist relaxation and drive EET. For other complexes, such as the phycobiliprotein PE545 complex, the exciton states are more widely separated due to the different chemical nature of the linear-tetrapyrole chromophores, and studies have highlighted the importance of higher frequency internal vibrational motions of the chromophores themselves enhancing EET.5 For example, ultrafast nonlinear spectroscopic studies exploring the early © 2016 American Chemical Society

stages of light harvesting in these biological systems, have measured prominent features arising from vibrational and vibronic coherences in the relaxation dynamics.6−11 Thus, properly accounting for the strength and frequency-dependence of electronic-nuclear interactions is crucial for simulating accurate dissipative open quantum system dynamics, where population transfer, dynamical dephasing, and decoherence processes of the quantum system are highly sensitive to the spectral density models. Although significant advances have been made experimentally to extract the electronic-nuclear couplings, they have only been determined for a very few light-harvesting systems, and the precise microscopic interpretation of many of these results remains ambiguous. Techniques such as hole-burning, fluorescence line narrowing (FLN), and the more recent difference fluorescence line narrowing (ΔFLN) measurements have elucidated the magnitude of the phonon sideband to characterize the coupling to the continuum of protein and solvent modes and the magnitude of the discrete vibrational peaks to estimate electronic couplings to the intramolecular vibrational modes of the chromophores. Unfortunately, the strongly overlapping bands of energy states of the different chromophores make it difficult to analyze and extract high-resolution and detailed sitespecific information. This has motivated the development of Received: June 29, 2016 Accepted: July 29, 2016 Published: July 29, 2016 3171

DOI: 10.1021/acs.jpclett.6b01440 J. Phys. Chem. Lett. 2016, 7, 3171−3178

Letter

The Journal of Physical Chemistry Letters theoretical and computational approaches to compute spectral densities12−16 primarily based on excitation energy correlation functions. These correlation functions are obtained using classical molecular dynamics (MD) trajectories simulated on an approximate ground state molecular mechanics (MM) surface, together with low-level quantum chemistry methods to estimate excitation energies. The level of approximation introduced for computational feasibility, however, can lead to unphysical results that disagree significantly with experimental findings. In this Letter, we describe a simple but practical and accurate alternative to the standard approaches for computing spectral densities. Our approach involves incorporating the influence of lower frequency, intermolecular components of the long-range interactions through analytical electrostatic interactions, together with a separate direct computation of the spectral density associated with the quantum intramolecular vibrational modes. The concepts underlying our theoretical development also enable us to identify the main reasons why standard molecular mechanics correlation function calculations provide spurious spectral density results. We demonstrate the reliability of our alternative approach by computing the site dependent spectral densities of the chromophores in the FMO complex and provide benchmark comparisons with detailed experimental results. For EET processes, the electronic states are often modeled in the site representation, where each site basis function describes the localized excited state of a chromophore with sitedependent transition energy at the Franck−Condon point. Thus, spectral density computation is considerably simplified, since spectral densities can be calculated for each localized excited chromophore, rather than for the excitonic state where the excited state may be delocalized over several sites. Assuming that the ground and excited state energies of the individual chromophores vary quadratically with the nuclear degrees of freedom, each chromophore site, α, is coupled to its own set of independent harmonic modes with frequencies ω(α) i , where subscript i labels the modes of the α chromophore’s environment with normal coordinates Qi(α). To a good approximation in these light harvesting systems, for a given chromophore, the shape of the excited state surface is identical to the ground state, but displaced in these normal mode with respect to the ground state coordinates by amounts d(α) i = 0, for all i. Thus, equilibrium reference geometry where Q(α) i once electronically excited, the nuclear configuration of the chromophore α environment is no longer in equilibrium and must reorganize to the displaced excited state minimum. The reorganization energy for nuclear mode i, is then 1 λi(α) = 2 ωi(α)2di(α)2 ,17 the dimensionless Huang−Rhys (HR) (α) factors are defined as S(α) = λ(α) and the gradient of the i i /ℏωi excited state surface in the direction of this mode at the = ω(α)2 d(α) Franck−Condon point is c(α) i i i . The site excitation (α) energy fluctuation, δϵ , arising from the nuclear motions of the local environment is then simply linear in these coordinates, (α) and δϵ(α) = ∑ic(α) i Qi . The system-environment interactions in these models are then determined entirely by the local bath mode frequencies and linear couplings c(α) ω(α) i i , and the computed correlation function of the excitation energy fluctuations can be expressed as

C(α)(t ) = ⟨δ ϵ(α)(t )δ ϵ(α)(0)⟩ =

∑ ci(α)2⟨Q i(α)(t )Q i(α)(0)⟩ i

= β∑ i

ci(α)2 ωi(α)

cos(ωi(α)t )

(1)

where β = 1/kBT. The spectral density summarizes this information as π J (α)(ω) ≡ ∑ (ci(α)2/ωi(α))δ(ω − ωi(α)) 2 i = π ∑ ωi(α)λi(α)δ(ω − ωi(α)) i

(2)

and is thus readily obtained from the classical excitation energy (α) correlation function as J(α)(ω) = βω ∫ ∞ 0 dtC (t) cos(ωt). 12,13,18−22 Hence, numerous spectral densities have been computed by simulating long MD trajectories governed by a parametrized MM force field and then by performing excited state quantum chemistry calculations in order to compute the time correlation functions. The spectral densities, J(α)(ω), are then extracted from thermally well-averaged classical correlation functions simply through a cosine transform, which does not require any additional factors,18 as discussed in detail in reference 32. Only when applying this spectral density in subsequent calculations of properties such as energy transfer rates or temperature-dependent spectra are quantum correction factors required.23 More recently, a growing number of studies have demonstrated14,21,22,24,25 that the spectral densities obtained from such calculations can be very sensitive to the parametrized MM force field since the MM potentials are not unique. In particular, MM parametrization of the chromophores with generalized force fields may result in different equilibrium (α) structures, {q(α) and iMM}, with different normal modes q different frequencies, Ω(α) from those obtained with ground i and excited state quantum chemistry calculations. The excitation energies computed with the quantum methods would then be inconsistent with the approximate MM ground state surface and lead to spectral densities with an inaccurate description of the electronic-nuclear interaction. In the following, we introduce a simple analysis, with detailed explanations in the first section of the Supporting Information, that relates how these three different potential energy surfacesthe QM excited state surface with the corresponding QM ground state surface, and MM ground state surfaceaffect the spectral density. These J(α) MM(ω) computed from excitation energy gap correlation functions along the MM harmonic ground state surface are shown to have the following approximate form: (α) JMM (ω)

⎛ ω (α)D (α) ⎞2 = πω ∑ ∑ ⎜⎜ k (α)ki ⎟⎟ λk(α)δ(ω − Ω(i α)) ⎠ i k ⎝ Ωi

(3)

where D(α) ki are the elements of the rotation matrix that relates the normal modes of the ground state ab initio quantum surface to the those of the approximate MM surface. The amplitudes of the spectral peaks, in general, involve linear combinations of the actual QM-determined reorganization energies with the combination coefficients determined by both the transformation between the MM surface and the QM surface normal modes, and the ratio of the squares of the frequencies of these modes. Qualitatively, the effect of these linear 3172

DOI: 10.1021/acs.jpclett.6b01440 J. Phys. Chem. Lett. 2016, 7, 3171−3178

Letter

The Journal of Physical Chemistry Letters

Figure 1. (a) Comparisons of (1) current calculations of J(ω) = Jintra(ω) + Jinter(ω) for BChl3 (blue curves) for FMO subunit A (solid line), FMO subunit B (dashed) and FMO subunit C dot-dashed with (2) Gaussian broadened intrachromophore spectral density extracted from ΔFLN experiments36 in light-gray filled curve; and (3) spectral density for BChl3 computed using MD excitation energy correlation function obtained with ZINDO/S-CHARMM protocol21,22 (black curve). (b) Same as (a) except now black curve presents spectral density for BChl3 computed using MD excitation energy correlation function obtained with DFT-Amber protocol.18 Also, panel b includes low-frequency spectral density obtained from FLN experimental data25,37 (dark pink curve).

energies on nuclear motion of the intermolecular environment. The fluctuation of the site energy can then be written as

combinations is to smear out the peaks in the true QM spectral density by adding different amounts of the various QM reorganization energies, rather than having pure peaks associated with each individual mode at its true frequency as in eq 2. In addition, J(α) MM(ω) has peaks at the MM frequencies, Ω(α) , rather than at the true frequencies, ω(α) i i , corresponding to the electronic structure potential surface, as with the true spectral density in eq 2. This dependence will not affect the spectral density as long as the two frequencies are similar. However, as will be later discussed, the mismatch between the QM and MM frequencies results in spurious spectral peaks with unphysical amplitudes. Using ab initio MD simulations employing density functional theory (DFT)26,27 in place of the MM force field, or adapting the force field on-the-fly,28 could alleviate the force field dependence, but running the long simulations necessary to obtain well-converged correlation functions needed to resolve low frequency modes of these PPCs could be prohibitively expensive with these methods. Employing a semiempirical QM/MM dynamics can make these calculations more tractable and enables consistent dynamics and excitation energy calculations.29 However, the reliability of the semiempirical parametrization scheme can be strongly system dependent.30 In the method developed here, the strategy for computing both the electron−phonon and -vibrational couplings borrows significantly from the ideas underlying the description of outersphere and inner-sphere electron transfer,17 which have also been suggested in the context of excitonic transport:13,22,31 the “short-range” inner-sphere interactions describe dependence of the excitation energy of a chromophore on its intramolecular vibrational degrees of freedom and are determined in an environment of static protein configurations. The “long-range” outer-sphere interactions capture the dependence of site

δ ϵ(α)(t ) ∼ δ ϵ(lrα)(t ) + δ ϵ(srα)(t )

(4)

and the excitation energy fluctuation correlation function in eq 1 will then contain direct terms, C(α)(t ) ∼ ⟨δ ϵ(srα)(0)δ ϵ(srα)(t )⟩ + ⟨δ ϵ(lrα)(0)δ ϵ(lrα)(t )⟩

(5)

(α) and cross-terms, ⟨δϵ(α) sr (0)δϵlr (t)⟩. The cross terms will vanish when short and long-range interactions are uncorrelated and is often a reliable approximation18 for light-harvesting complexes. The spectral density then only contains short-range, intramolecular vibrational contributions and long-range intermolecular vibrational (phonon) contributions:

(α) (α) J (α)(ω) ∼ Jinter (ω) + Jintra (ω) 32

(6) (α) Jinter (ω)

Rivera et al. computed site-dependent by implementing the charge density coupling (CDC) method formulated by Renger and co-workers,22,31,33,34 to calculate the changes in the long-range contributions to the site energies along MD trajectories for the FMO complex. The ground and excited state electronic densities of the chromophores in this procedure are replaced with atomic (ESP) point charges that remain fixed when calculating the long-range point-charge Coulomb interaction between the both electronic states of the chromophores and the fluctuating environment. Thus, this approximation avoids running quantum chemistry calculations for each snapshot along the nuclear trajectory. Because the changes in the electronic density arising from intrachromophore nuclear distortions are neglected, essentially only intermolecular effects are taken into account with this approach. Qualitatively, the shape, magnitude and relaxation time of the intermolecular spectral densities computed in this 3173

DOI: 10.1021/acs.jpclett.6b01440 J. Phys. Chem. Lett. 2016, 7, 3171−3178

Letter

The Journal of Physical Chemistry Letters way match reasonably well with the Debye form, JDebye(ω) = 2λωτ/(1 + ω2τ2), employed in many model studies,35 with λ(α) inter = 1/π ∫ dωJ(α) inter(ω)/ω, the typical reorganization energy of magnitude 35 cm−1, and τ and the characteristic relaxation time of 50 fs. Quantitatively, differences from site to site due to structural variations in the local environments of the different chromophores leads to variation from 10 to 40 cm−1 for the total intermolecular reorganization energy component. Constructing the intramolecular contribution, J(α) intra(ω), which complements J(α) inter(ω) in eq 6, only requires calculating the gradient of the excited state surface with respect to the normal mode coordinates of the chromophores since, as mentioned earlier, the electron−nuclear coupling {ci} is simply the gradient of the excited surface at the Franck−Condon point. These vectors can be easily computed by performing ground state geometry optimizations and subsequent normal-mode analysis with reliable quantum chemistry methods. The MD trajectory generates configurations of the full pigment−protein system that sample different locally distorted geometries of the chromophores. Ground state quantum chemistry geometry optimization of each chromophore in the presence of the static electric field of the fixed MM partial charge geometry provides chromophore configurations that sample different local minima of the quantum mechanical potential surface. Assuming Duschkinsky effects are negligible, the excited state forces along the ground state normal mode coordinates can be extracted either numerically or with analytical gradients that are available for different electronic states for many electronic structure methods in most quantum chemistry packages. The Cartesian analytical gradients can then easily be manipulated to give forces in terms of normal mode coordinates. Full details of the implementation of the approach, as well as computational details, can be found in the Supporting Information. An overview of the method appears in reference 39. The performance of the method outlined here has been tested for the well-known FMO trimer antenna complex, where the network of chromophores in each of the β-barrel subunits, A−C, is composed of eight bacteriochlorophyll A chromophores labeled as BChl 1−8. All results reported here were obtained for a single, well-equilibrated MD configuration (see the SI of reference 32 for details), which are validated by comparing with results from ΔFLN experiments. It is important to note that the experimental measurements focus only on the distinct and well-separated lowest energy spectral peak, mostly localized on BChl3, in order to avoid spectral congestion with the many overlapping absorptions at higher energies. Thus, the ΔFLN experiments probe mainly nuclear couplings of the BChl3 chromophore with negligible contributions from other BChls. The total spectral densities for the three BChl3 chromophores of our FMO trimer configuration, obtained as the sum of an intermolecular contribution32 and the instantaneous intramolecular spectral density, according to eq 6, are presented in Figure 1. Our computed site-dependent instantaneous spectral densities for all the eight BChls for the FMO trimer are presented in Figure S2 in the Supporting Information. The intramolecular spectral densities are constructed by artificially broadening each discrete peak of magnitude ci2/ωi with Gaussian functions with standard deviation of σ = 7 cm−1. These broadened site-dependent spectral densities were then normalized to conserve the calculated total reorganization energy (reported in Table 1 of the Supporting Information) of that site. This broadening was applied to mimic the various

line-broadening effects (inhomogeneous broadening, puredephasing, spectral diffusion, etc.) that are averaged over in experimental measurements. The data presented using the different blue line types in Figure 1 are our calculated spectral densities for the BChl3 chromophores in the different subunits of the trimer configuration. These results are compared with other representative work: In panel a, the black curve is the spectral density obtained by computing ZINDO/S excitation energy fluctuations along the CHARMM ground state surface.21,22 In panel b, the black curve gives the BChl3 spectral density obtained using (TD)B3LYP/3-21G excitation energies computed along Amber MM trajectories.18 These results will be referred as ZINDO-CHARMM and DFT-Amber, respectively, in the following discussion. These computed BChl3 spectral densities are also compared with two different experimentally determined spectral densities: the gray filled curve in both panels gives the spectral density obtained from the reported experimental ΔFLN HR factors,36 which have been Gaussian broadened in exactly the same manner as our current calculated results. Finally, the dark pink curve reported in the lowfrequency region of panel b corresponds to the analytic form of the spectral density25 that has been fit to the FLN measurements25,37 and are discussed in reference 18. The similarity in shape and magnitude of the two experimental spectral densities rationalizes the estimated broadening we have employed to model the experimental spectral densities from the discrete experimental HR factors. This approach for constructing the experimental spectral density by simply Gaussian broadening the ΔFLN HR factors in the same way as we do for our calculated results, however, neglects the low frequency phonon response manifested as the broad phonon side bands.36 The removal and neglect of this phonon piece in the construction of the ΔFLN spectral density accounts for the main differences in the low-frequency region between these results and the FLN spectral density:25,37 the positions of the various low frequency peaks agree well between these two sets of experimental results, but the amplitudes differ at the lowest frequencies below about 200 cm−1. Between 250 and 500 cm−1, however, we see that the magnitudes of these two different experimental estimates of the spectral density agree quite closely. Lastly, the fit for FLN spectral density has only been optimized for modes below 500 cm−1; thus, only ΔFLN spectral density is shown for the higher frequency regions where this experimental approach has enhanced resolution. A comprehensive review of FLN and ΔFLN38 spectroscopy is given in reference 38. The ΔFLN spectral density above 500 cm−1 identifies many strongly coupled intrachromophore vibrational modes for the lowest energy electronic eigenstate predominantly localized on BChl3. The correspondence of peak positions and their relative magnitudes between the ΔFLN measurements and our calculated normal mode instantaneous local spectral densities for BChl3 shown in Figure 1 is generally very good. In this same higher frequency region, the calculated peak positions from ZINDO−CHARMM and DFT-Amber for the BChl3 spectral density show essentially no correspondence with those determined experimentally, and the magnitudes of the calculated features are typically an order of magnitude too small, except at the highest frequencies, where these methods give unphysical, excessively large amplitude responses. In the lower frequency region (0−500 cm−1) in Figure 1 we see that, because the ZINDO−CHARMM spectral densities 3174

DOI: 10.1021/acs.jpclett.6b01440 J. Phys. Chem. Lett. 2016, 7, 3171−3178

Letter

The Journal of Physical Chemistry Letters

Figure 2. Comparison of potential surface properties and spectral densities calculated for a gas phase BChl A molecule. The inset shows the correlation of the ground state QM surface vibrational frequency (here computed with B3LYP and estimated using finite difference as outlined in the text) with MM surface normal mode vibrational frequency in wavenumbers. The dark-pink filled circles indicate modes with significant excitonicvibrational coupling, i.e., reorganization energies larger than 1 cm−1. The spectral densities shown in the main part of the figure are calculated using (TD)B3LYP/3-21G with eq 2 (dark pink) and eq 3 (black). The arrows connecting spikes of the same colors under the black and pink spectral densities show examples of how higher frequency modes of the MM surface (black curve) must be substantially shifted to lower frequencies: the ground state QM surface is, in general, significantly softer along these vibrational coordinates and must be shifted to match the peaks in the true QM surface spectral density (pink curve). Results comparing analogous calculations using the ZINDO/S semiempirical electronic structure method are presented in the Supporting Information.

were obtained by first fitting the excitation energy fluctuation correlation functions,21 the results are generally smoothed out and rather featureless, compared to the highly structured experimental spectral densities. For some chromophores (see the data presented in Figure S2) these ZINDO−CHARMM spectral density results on average capture reasonable magnitudes, but they can be significantly larger, e.g., for BChl7. The DFT-Amber averaged spectral densities shown for BChl3 in panel (b) (see also BChl1 comparison in Figure S2) in this low frequency region generally provide a better description of the frequency dependence of the overall magnitudes and resolve some features seen in the experiments, but several peaks are poorly resolved, most likely due to limitation of obtaining well-averaged long-time correlation functions. In marked contrast, our composite method for calculating local spectral densities in the lower frequency region show very good correspondence for mode frequency and spectral density magnitude for BChl3. The independently calculated lower frequency BChl3 peaks for the A, B, and C monomers all show large magnitude peaks just below 200 cm−1 that match well with the FLN experimental results. Other peaks between 250 and 400 cm−1 show more variation from monomer to monomer, but the features are generally unvarying and are in good qualitative agreement with the experiments. The other chromophores (see Figure S2) show slight variations on the basic structural features outlined above. These differences are more clearly observed in the Huang−Rhys (HR) factors reported in Figure S1 of the Supporting Information, where these measures of chromophore electron−nuclear coupling show some local variations between the same chromophores in the three different subunit and also between different chromophores. The differences from chromophore to chromophore and subunit to subunit give an indication of the dispersion estimated in these quantities associated with the different local structures sampled from the room temperature thermal distribution of the protein. They reflect the influence of real differences between the local BChl environments. These differences are also clearly discernible in Figure S2 and in Figure S5 in the Supporting Information, where the sensitivity of the current spectral densities to the various density functionals is explored. Overall, however, the

differences are small, and the results show good agreement with the experimental HR factors. The deviations between the averaged spectral densities computed from MD correlation functions and the experimental results of our normal mode local spectral densities are summarized in the cumulative distribution plots presented Figure S3 of the Supporting Information. We can understand the origin of the discrepancies identified here by comparing the expression for the ab initio spectral density in eq 2 with the simple model of the MD correlation function calculations summarized in eq 3. As discussed earlier, the MD correlation function spectral density will have peaks at the frequencies of the MM surface, and the magnitudes of these features will in general be linear combinations of the true ab initio reorganization energies determined by the square of the ratio of the ab initio and MM frequencies multiplied by the matrix elements of the transformation between the normal modes of these different surfaces. This shortcoming that gives rise to erroneous spectral density peaks obtained with the correlation function method is demonstrated by applying the approach we have introduced earlier in this Letter and in the first section of the Supporting Information. For a gas-phase bacteriochlorophyll A molecule, geometry optimization and normal-mode analysis on the CHARMM potential energy surface were performed to obtain CHARMM normal mode coordinates. Next, the properties of the relevant ab initio surfaces around the minimum of the CHARMM surface are obtained using a finite difference approach to calculate both the local curvatures of the QM ground and excited states and the excited state gradients (obtained using TDDFT and ZINDO/S) by displacing along the CHARMM normal mode eigenvectors. This approach enables analysis of the ground and excited state frequencies and the reorganization energies on the ab initio surfaces, where these surfaces correspond to the those characterized in the correlation function approach that samples geometries from the MM surface minima using MD. The results, presented in Figure 2, are for density functional theory using (TD) B3LYP/3-21G, and in Figure S4 in the Supporting Information we also report ZINDO/S semiempirical electronic structure results. These two methods were chosen in an attempt to model the approaches 3175

DOI: 10.1021/acs.jpclett.6b01440 J. Phys. Chem. Lett. 2016, 7, 3171−3178

Letter

The Journal of Physical Chemistry Letters used in recent MD correlation function studies21,22 and the results presented in Figure 2 show the influence of the nontrivial mismatch between the MM and ab initio potential surfaces on computed spectral densities. The spectral density results reported in Figure 2 do not include the influence of the Duschinsky-like rotations between the MM and QM normal modes (see discussion in the Supporting Information for details) and thus only highlight the effects arising from the differences in local curvature of the MM and QM surfaces. From these calculations, we find that the QM curvatures, {ωi}, can be considerably different from those of the MM local surface curvatures, {Ωi}. While the B3LYP/3-21G and CHARMM frequencies are generally well correlated with each other for frequencies less than 1200 cm−1, consistent deviations are found for higher frequency modes (see inset figures). As presented in Figure S4 of the Supporting Information, the match between the CHARMM and ZINDO surfaces is noticeably worse. The influence of these differences in the surface morphology on the calculated model spectral densities is apparent from the curves in the main panel of Figure 2, which compare the MM spectral densities, JMM(ω), (black curves) computed using eq 3, assuming that Dki = δki (see Supporting Information for details), with the (TD)B3LYP/3-21G electronic structure method. These spectral densities have their peaks at the MM surface frequencies {Ωi}. The spectral density, J(ω), (dark pink curve), obtained using eq 2, on the other hand, has its peaks at the QM surface frequencies {ωi}. These differences in the MM and QM surfaces mean that many of the features in the high frequency region of the model MD correlation function spectral density arise from the inaccurate MM potential energy surface. Thus, many of the higher frequency peaks just below 1500 cm−1 in the J(ω) (dark-pink curve) are artificially shifted to significantly higher frequencies to give the relocated and rescaled features that appear in the JMM(ω) (black curve). The arrows in Figure 2 for example, indicate representative peaks in JMM(ω) that shift to give the true peak positions in J(ω). We find that the high frequency features in the ZINDO− CHARMM and DFT-Amber results21,22 in Figure 1 arise from this problem: analogous with the findings outlined above for our model isolated BChl calculation, the unphysical high frequency features appearing in the MD spectral densities are shifted up by several hundred wavenumbers from where they should appear if the true QM surface frequencies were used, and give consistent results with both our total spectral densities and the experimental results for the FMO complex. In this Letter, we have presented an accurate and efficient approach for calculating spectral densities of electronically excited molecules in complex environments that captures both short-range intrachromophore, and longer range intermolecular electronic-nuclear interactions. Separate estimation of the intermolecular component that does not require expensive electronic structure calculations for correlation functions, and the intramolecular spectral density that requires only a few basic electronic structure calculations, significantly reduces computational cost of our approach compared to other widely used methods. These essential computations include geometry optimizations and normal-mode analysis using standard ground state quantum chemistry calculations, along with the analytical gradients of the excited state calculated consistently with the ground state surface. Ambiguities and artifacts introduced from employing MM potentials are thus removed. The reliability of this method is particularly evident for the high frequency

regions of the spectral density, where we have shown that inconsistency between widely used general purpose MM ground state model potential surfaces and the QM excited state surface leads to unrealistic estimates of reorganization energies. Though these erroneous results may not be so significant for complexes such as FMO, where the small exciton energy gaps are not resonant with high-frequency nuclear modes, realistic description of modes in these regions is crucial for systems such as PE545 and the exciton/charge-transfer states of the reaction center, where the relevant energy gaps are significantly larger and may be resonant with the intramolecular modes that could promote enhanced charge and energy transfer. The method presented here shows significant promise for accurately describing computationally challenging systems such as the reaction center, since the significant reduction of computational time allows the possibility of using higher level and systematically more accurate electronic structure methods necessary to describe chromophores in these and other states such as dark states that are difficult to probe experimentally. Moreover, the method introduced here can easily be applied to not only light-harvesting systems, but for other systems, where, in particular, the ground and excited state potential energy surfaces significantly differ and Duschkinsky effects become important. We conclude by emphasizing that the spectroscopic measurements that are employed in global fitting of parametrized system-bath model Hamiltonians, including spectral densities, generally reflect averages of various configurations sampled from the energy landscape of the complex. The approach introduced here not only gives a convenient way to analyze the local spectral densities of individual chromophores into their intermolecular and intrachromophore components, but also, because of its efficiency, it can in principle be implemented to explore how instantaneous system-bath coupling terms in model Hamiltonians change as an MD trajectory samples different regions of the pigment - protein potential energy landscape. One could thus use the properties of the different local basins sampled along MD trajectories to define these different “instantaneous” EET Hamiltonians and provide meaningful averaged dynamics with the appropriate separation of time scales. Such an approach is more welldefined and reasonable than using an averaged spectral density to describe the fast EET processes as is often done for computational convenience. Our approach for obtaining instantaneous local spectral densities thus produces results that may be combined with information about site excitation energy fluctuations with conformation, and employed in accurate model quantum dynamics methods39 to give a highly versatile approach for analyzing ultrafast single molecule spectroscopy experiments.40

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.6b01440. Simple analytic model for molecular mechanics correlation function calculation of spectral densities; computational details of simulations; detailed results for Huang− Rhys factors and cumulative spectral densities, spectral densities for chlorophyll dimer computed with ZINDO method; and a detailed study of influence of choice of functional on computed spectral densities (PDF). A 3176

DOI: 10.1021/acs.jpclett.6b01440 J. Phys. Chem. Lett. 2016, 7, 3171−3178

Letter

The Journal of Physical Chemistry Letters

■

(12) Damjanovic, A.; Kosztin, I.; Kleinekathöfer, U.; Schulten, K. Excitons in a Photosynthetic Light-Harvesting System: A Combined Molecular Dynamics, Quantum Chemistry, and Polaron Model Study. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2002, 65, 031919−031943. (13) Olbrich, C.; Strümpfer, J.; Schulten, K.; Kleinekathöfer, U. Theory and Simulation of the Environmental Effects on FMO Electronic Transitions. J. Phys. Chem. Lett. 2011, 2, 1771−1776. (14) Renger, T.; Klinger, A.; Steinecker, F.; Schmidt am Busch, M.; Numata, J.; Müh, F. Normal Mode Analysis of the Spectral Density of the Fenna-Matthews-Olson Light-Harvesting Protein: How the Protein Dissipates the Excess Energy of Excitons. J. Phys. Chem. B 2012, 116, 14565−14580. (15) Jurinovich, S.; Viani, L.; Curutchet, C.; Mennucci, B. Limits and Potentials of Quantum Chemical Methods in Modelling Photosynthetic Antennae. Phys. Chem. Chem. Phys. 2015, 17, 30783−30792. (16) Olbrich, C.; Strümpfer, J.; Schulten, K.; Kleinekathöfer, U. Quest for Spatially Correlated Fluctuations in the FMO LightHarvesting Complex. J. Phys. Chem. B 2011, 115, 758−764. (17) Nitzan, A. In Chemical Dynamics in Condensed Phases; Oxford University Press: Oxford, U.K., 2006. (18) Valleau, S.; Eisfeld, A.; Aspuru-Guzik, A. On the Alternatives for Bath Correlators and Spectral Densities from Mixed QuantumClassical Simulations. J. Chem. Phys. 2012, 137, 224103−224116. (19) Aghtar, M.; Strü mpfer, J.; Olbrich, C.; Schulten, K.; Kleinekathöfer, U. Different Types of Vibrations Interacting with Electronic Excitations in Phycoerythrin 545 and Fenna Matthews Olson Antenna Systems. J. Phys. Chem. Lett. 2014, 5, 3131−3137. (20) Viani, L.; Corbella, M.; Curutchet, C.; O’Reilly, E. J.; OlayaCastro, A.; Mennucci, B. Molecular Basis of the Exciton-Phonon Interactions in the PE545 Light-Harvesting Complex. Phys. Chem. Chem. Phys. 2014, 16, 16302−16311. (21) Chandrasekaran, S.; Aghtar, M.; Valleau, S.; Aspuru-Guzik, A.; Kleinekathöfer, U. Comparison between CHARMM & AMBER and ZINDO vs. TDDFT. J. Phys. Chem. B 2015, 119, 9995−10004. (22) Wang, X.; Ritschel, G.; Wuster, S.; Eisfeld, A. Open Quantum System Parameters for Light-Harvesting Complexes from Molecular Dynamics. Phys. Chem. Chem. Phys. 2015, 17, 25629−25641. (23) Skinner, J. L.; Park, K. Calculating Vibrational Energy Relaxation Rates from Classical Molecular Dynamics Simulations: Quantum Correction Factors for Processes Involving Vibration-Vibration Energy Transfer. J. Phys. Chem. B 2001, 105, 6716−6721. (24) Kim, C. W.; Park, J. W.; Rhee, Y. M. Effect of Chromophore Potential Model on the Description of Exciton-Phonon Interactions. J. Phys. Chem. Lett. 2015, 6, 2875−2880. (25) Kreisbeck, C.; Kramer, T. Long-Lived Electronic Coherence in Dissipative Exciton Dynamics of Light-Harvesting Complexes. J. Phys. Chem. Lett. 2012, 3, 2828−2833. (26) Sisto, A.; Glowacki, D. R.; Martinez, T. J. Ab Initio Nonadiabatic Dynamics of Multichromophore Complexes: A Scalable GraphicalProcessing-Unit-Accelerated Exciton Framework. Acc. Chem. Res. 2014, 47, 2857−2866. (27) Shim, S.; Rebentrost, P.; Valleau, S.; Aspuru-Guzik, A. Theory and Simulation of the Environmental Effects on FMO Electronic Transitions. Biophys. J. 2012, 102, 649−660. (28) Higashi, M.; Saito, S. Quantitative Evaluation of Site Energies and Their Fluctuations of Pigments in the Fenna−Matthews−Olson Complex with an Efficient Method for Generating a Potential Energy Surface. J. Chem. Theory Comput. 2016, DOI: 10.1021/ acs.jctc.6b00516. (29) Rosnik, A. M.; Curutchet, C. Theoretical Characterization of the Spectral Density of the Water-Soluble Chlorophyll-Binding Protein from Combined Quantum Mechanics/Molecular Mechanics Molecular Dynamics Simulations. J. Chem. Theory Comput. 2015, 11, 5826− 5837. (30) Gerber, R. B.; Shemesh, D.; Varner, M. E.; Kalinowski, J.; Hirshberg, B. Ab initio and semi-empirical Molecular Dynamics simulations of chemical reactions in isolated molecules and in clusters. Phys. Chem. Chem. Phys. 2014, 16, 9760−9775.

spread sheet containing all computed spectral densities and Huang−Rhys factors as presented in Figures S1 and S2 is also available (XLSX).

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We gratefully acknowledge support for this research from the National Science Foundation under grant number CHE1301157. We acknowledge the computational resources provided by Boston University’s Office of Research Computing, Information Systems and Technology. We are grateful to Eva (Tina) Rivera, Daniel Montemayor, Marco Masia, Alexander Kunitsa, and Pengfei Huo for valuable discussions.

■

REFERENCES

(1) Romero, E.; Augulis, R.; Novoderezhkin, V. I.; Ferretti, M.; Thieme, J.; Zigmantas, D.; van Grondelle, R. Quantum Coherence in Photosynthesis for Efficient Solar-Energy Conversion. Nat. Phys. 2014, 10, 676−682. (2) Turner, D. B.; Dinshaw, R.; Lee, K. K.; Belsley, M. S.; Wilk, K. E.; Curmi, P. M. G.; Scholes, G. D. Quantitative Investigations of Quantum Coherence for a Light-Harvesting Protein at Conditions Simulating Photosynthesis. Phys. Chem. Chem. Phys. 2012, 14, 4857− 4874. (3) Tiwari, V.; Peters, W. K.; Jonas, D. M. Electronic Resonance with Anticorrelated Pigment Vibrations Drives Photosynthetic Energy Transfer Outside the Adiabatic Framework. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 1203−1208. (4) Christensson, N.; Kauffmann, H. F.; Pullerits, T.; Mancal, T. Origin of Long-Lived Coherences in Light-Harvesting Complexes. J. Phys. Chem. B 2012, 116, 7449−7454. (5) Kolli, A.; O’Reilly, E. J.; Scholes, G. D.; Olaya-Castro, A. The Rundamental Role of Quantized Vibrations in Coherent Light Harvesting by Cryptophyte Algae. J. Chem. Phys. 2012, 137, 174109−174124. (6) Engel, G. S.; Calhoun, T. R.; Read, E. L.; Ahn, T. K.; Mancal, T.; Cheng, Y. C.; Blankenship, R.; Fleming, G. Evidence for Wavelike Energy Transfer through Quantum Coherence in Photosynthetic Systems. Nature 2007, 446, 782−786. (7) Panitchayangkoon, G.; Hayes, D.; Fransted, K.; Caram, J. R.; Harel, E.; Wen, J.; Blankenship, R. E.; Engel, G. S. Long-lived Quantum Coherence in Photosynthetic Complexes at Physiological Temperature. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 12766−12770. (8) Collini, E.; Wong, C. Y.; Wilk, K. E.; Curmi, P. M. G.; Brumer, P.; Scholes, G. D. Coherently Wired Light-Harvesting in Photosynthetic Marine Algae at Ambient Temperature. Nature 2010, 463, 644−647. (9) Wong, C. Y.; Alvey, R. M.; Turner, D. B.; Wilk, K. E.; Bryant, D. A.; Curmi, P. M. G.; Silbey, R. J.; Scholes, G. D. Electronic Coherence Lineshapes Reveal Hidden Excitonic Correlations in Photosynthetic Light Harvesting. Nat. Chem. 2012, 4, 396−404. (10) Schlau-Cohen, G. S.; Calhoun, T. R.; Ginsberg, N. S.; Ballottari, M.; Bassi, R.; Fleming, G. R. Spectroscopic Elucidation of Uncoupled Transition Energies in the Major Photosynthetic Light-Harvesting Complex, LHCII. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 13276− 13281. (11) Schlau-Cohen, G. S.; Ishizaki, A.; Calhoun, T. R.; Ginsberg, N. S.; Ballottari, M.; Bassi, R.; Fleming, G. R. Elucidation of the Timescales and Origins of Quantum Electronic Coherence in LHCII. Nat. Chem. 2012, 4, 389−395. 3177

DOI: 10.1021/acs.jpclett.6b01440 J. Phys. Chem. Lett. 2016, 7, 3171−3178

Letter

The Journal of Physical Chemistry Letters (31) Adolphs, J.; Renger, T. How Proteins Trigger Excitation Energy Transfer in the FMO Complex of Green Sulfur Bacteria. Biophys. J. 2006, 91, 2778−2797. (32) Rivera, E.; Montemayor, D.; Masia, M.; Coker, D. F. Influence of Site-Dependent Pigment-Protein Interactions on Excitation Energy Transfer in Photosynthetic Light Harvesting. J. Phys. Chem. B 2013, 117, 5510−5521. (33) Madjet, M. E.; Abdurahman, A.; Renger, T. Intermolecular Coulomb Couplings from Ab Initio Electrostatic Potentials: Application to Optical Transitions of Strongly Coupled Pigments in Photosynthetic Antennae and Reaction Centers. J. Phys. Chem. B 2006, 110, 17268−17281. (34) Jing, Y.; Zheng, R.; Li, H.-X.; Shi, Q. Theoretical Study of the Electronic Vibrational Coupling in the Qy States of the Photosynthetic Reaction Center in Purple Bacteria. J. Phys. Chem. B 2012, 116, 1164− 1171. (35) Ishizaki, A.; Fleming, G. R. Theoretical Examination of Quantum Coherence in a Photosynthetic System at Physiological Temperature. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 17255−17260. (36) Ratsep, M.; Freiberg, A. Electron-Phonon and Vibronic Couplings in the FMO Bacteriochlorophyll A Antenna Complex Studied by Difference Fluorescence Line Narrowing. J. Lumin. 2007, 127, 251−259. (37) Wendling, M.; Pullerits, T.; Przyjalgowski, M. A.; Vulto, S. I. E.; Aartsma, T. J.; van Grondelle, R.; van Amerongen, H. ElectronVibrational Coupling in the Fenna-Matthews-Olson Complex of Prosthecochloris aestuarii Determined by Temperature-Dependent Absorption and Fluorescence Line-Narrowing Measurements. J. Phys. Chem. B 2000, 104, 5825−5831. (38) Pieper, J.; Freiberg, A. In The Biophysics of Photosynthesis; Golbeck, J., van der Est, A., Eds.; Springer: Dordrecht, The Netherlands, 2014; Chapter 2, pp 45−78. (39) Lee, M. K.; Huo, P.; Coker, D. F. Semiclassical Path Integral Dynamics: Photosynthetic Energy Transfer with Realistic Environment Interactions. Annu. Rev. Phys. Chem. 2016, 67, 639−668. (40) Schlau-Cohen, G. S.; Wang, Q.; Southall, J.; Cogdell, R. J.; Moerner, W. E. Single-Molecule Spectroscopy Reveals Photosynthetic LH2 Complexes Switch Between Emissive States. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 10899−10903.

3178

DOI: 10.1021/acs.jpclett.6b01440 J. Phys. Chem. Lett. 2016, 7, 3171−3178