Effect of Chromophore Potential Model on the Description of Exciton

Letter pubs.acs.org/JPCL

Effect of Chromophore Potential Model on the Description of Exciton−Phonon Interactions Chang Woo Kim, Jae Woo Park, and Young Min Rhee* Center for Self-Assembly and Complexity, Institute for Basic Science (IBS), Pohang 790-784, Korea Department of Chemistry, Pohang University of Science and Technology (POSTECH), Pohang 790-784, Korea S Supporting Information *

ABSTRACT: System−bath interactions in nonadiabatic simulations are often depicted by first performing molecular dynamics calculations and then by evaluating excitation energies at the trajectory snapshots. Usually, molecular mechanics models and quantum chemical calculations are used in a mixed manner toward a trade-off between efficiency and accuracy. Here we investigate how this mixing scheme affects that depiction by using various potential energy surfaces (PESs) of coumarin-153 chromophore, with the help of interpolated PESs that can closely match the accuracies of quantum chemical calculations. We find that although spectral densities are computed only with second stage data the PES characteristics during the first sampling stage can still prevail in the densities, with limited influences on related reorganization energies. Our results suggest that using the mixed scheme can be acceptable when dynamics is mainly governed by the integrated effect of all phonon modes, but care must be taken for understanding detailed effects from individual modes.

N

snapshots. Through this, the spectral density can be obtained from the time-dependent fluctuations of the excitation energy. Within the regime of linear response,9,10 separately performing two different types of calculations is completely valid; however, a practical consideration often enforces an approximation on the calculations. Because the excitation energy calculations in the second stage can trivially be performed in a parallel fashion, a relatively high level of method can be adopted without too much difficulty; however, the first stage is not as parallelizable and an economically acceptable approach has to be applied. Thus, the natural consideration behind the combination of these two different stages, namely, the MD simulations and the excitation energy calculations, is attaining a balance between efficiency and accuracy. This is the simple reason for separately adopting inexpensive MM models and reliable quantum chemical (QC) approaches in the two stages by convention, while, in principle, both can be performed only with QC calculations or only with MM modeling. Namely, a QC-only approach will be extremely slow for obtaining long trajectories, and an MM-only approach may be questioned regarding its accuracy, especially toward describing excited-state energies. These issues are cleverly avoided by adopting the mixed approach with both MM and QC. Then, for its success, we can easily imagine that the accordance between the two models will be important. Up to this point, however, the effect of the model agreement or disagreement has not been systematically studied yet. The bottleneck toward such verifications is mainly the computational cost associated with adopting higher level approach during MD simulations.

onadiabatic processes involving multiple quantum states take place in many interesting chemical phenomena such as exciton migration, photodynamic relaxation, and electron/ proton transfer. It is known that the interactions of the quantum system with its environment can play important roles in these processes. Because treating large systems in a fully quantum mechanical manner is practically impossible, researchers have developed various approximate methods. For electronic transitions, the nuclei are often considered as classical particles, and only the electronic subsystem is treated quantum mechanically. For this purpose, the collection of the vibrational modes of the entire system is often modeled as harmonic phonon bath from the viewpoint of the electronic subsystem. The spectral density is the effective profile of system−bath or exciton−phonon coupling in the frequency domain, and it provides information about how the system affects and is affected by the nuclear motions.1 Indeed, it serves as a key quantity for describing exciton−phonon interactions in many density matrix-based approaches.2−5 The spectral density can be obtained both by experiments6 and by computations.7,8 In the case of realistic systems, computational characterizations can be achieved by performing molecular dynamics (MD) trajectory simulations. For example, a combination of MD simulations with molecular mechanics (MM) modeling and excitation energy calculations with the time-dependent density functional theory (TDDFT)8 or semiempirical methods7 has been adopted to generate the spectral densities of Fenna− Matthews−Olson (FMO) complex and the reaction center complex. In this combined approach, equilibrium MD simulations are performed first to sample the trajectory snapshots on the ground-state potential energy surface (PES), and the excitation energies are next calculated along the © XXXX American Chemical Society

Received: May 30, 2015 Accepted: July 6, 2015

2875

DOI: 10.1021/acs.jpclett.5b01141 J. Phys. Chem. Lett. 2015, 6, 2875−2880

Letter

The Journal of Physical Chemistry Letters

with the interpolated potential closely matched experimental results, except some overestimation of the Stokes shift (interpolation: ∼24 000 cm−1 (ref 14), experiment: ∼6000 cm−1 (ref 17)) due to limitations of the QC method adopted for the potential energy surface construction (SOS-MP219 for the ground state and SOS-CIS(D0)20 for the excited state) and due to the use of nonpolarizable solvent model.21 Nevertheless, this limitation does not pose a problem in our analysis at all as we can assume that the interpolated potential is the exact surface that the MM model is trying to mimic as closely as possible. Simply, we can elucidate the deviation of the MM + interpolation mixed scheme from the MM-only scheme for describing the characteristics of the system−bath interactions. For all simulations, we have employed a modified version of the GROMACS 4.0 package22 with our in-house library for interpolation. As already explained, the interpolated surface (to be referred as “PES-i”) was adopted as an efficient replacement of any onthe-fly QC calculations. An ultimate accuracy can be achieved by adopting PES-i in both stages, namely, for the first stage of trajectory simulations and for the second stage of excitation energy calculations. When it is necessary to explicitly distinguish the ground and the excited state interpolated surfaces from each other, we will denote them as PES-ig and PES-ie, respectively. When this distinction is not needed, we will just use PES-i. Toward the trajectory simulation stage, we also adopted three different MM models with varying degrees of agreements with the interpolated potential. The first one (PES-1) was created by the Hessian matching process,23 which sets the second derivative of the MM potential to best reproduce the QC reference at the potential energy minimum point. Thus, at the equilibrium geometry, both the interpolated and MM potentials are closely consistent with each other up to the second order. The second MM model (PES-2) was generated by replacing the stretching and bending force constants of PES-1 by Generalized Amber Force Field (GAFF) parameters.24 Thus, PES-2 will be consistent with the interpolated surface only at the equilibrium geometry, namely, to the first order. The last MM model (PES-3) was set up by further replacing the equilibrium bond lengths and angles of PES-2 by GAFF parameters. This surface will, of course, have the least consistency with the reference interpolated surface. Of the three models, PES-1 was taken from a previous work14 and then propagated to generate PES-2 and then PES-3. We adopted the same ground-state nonbonded interaction parameters for PES-ig, PES-1, PES-2, and PES-3. Of course, a separate set of excited-state nonbonded parameters were used for PES-ie. (See ref 14 for the list of these parameters.) For all simulations, one coumarin 153 molecule was solvated with methanol in a cubic box of 33 Å side length, and the system energy was minimized by the steepest descent method. The molecular model for methanol was constructed based on the OPLS all atom force field.25 All simulations were performed with the isothermal−isobaric (NPT) ensemble condition of 300 K and 1 bar. For this, the system was first equilibrated for 500 ps using Berendsen’s weak-coupling thermostat and barostat,26 followed by another equilibration process of 500 ps using velocity rescale thermostat27 and Parrinello−Rahman barostat.28 After this, 10 ns production runs were performed under the same temperature and pressure conditions. A time step of 1 fs was adopted, and the snapshots of the system were recorded at every 5 fs. We measured the energy difference between the ground- and the excited-state PESs along the

Nonetheless, a verification process will still be highly important, especially when one needs to carefully consider electron− phonon coupling. Indeed, molecular complexes with relatively large chromophores, which are heavily studied in recent years,11,12 will likely be in that category, as the chromophore internal vibrations will act as an important part of the phonon bath with potentially stronger coupling to the electronic states generated by the chromophores themselves. We report on the extent to which the model accordance/ discrepancy during the two stages affect the reliability of exciton−phonon coupling aspect for a solvated chromophore system. Instead of adopting a costly QC model in all computations, we will adopt interpolated potentials that can closely match the accuracy of corresponding QC surfaces. By visually inspecting and comparing the spectral densities obtained from simulations with different mixing schemes, we will trace how the model discrepancy affects the results. We will test a series of MM models with varying degrees of discrepancies with the reference interpolated model for more detailed analyses. Of course, results from the interpolation-only scheme will be adopted as references for all comparisons. The interpolation scheme reconstructs the molecular electronic PES from a set of QC data with energies and their derivatives, evaluated at a wide range of dynamically relevant geometrical points.13 (See the Supporting Information for its mathematical formulations.) Conceptually, the interpolated surface can be thought of as a conjoint of locally accurate piecewise potentials. In the mathematical limit of an infinite data set, its accuracy can match the level of the reference QC calculations. In fact, its reliability toward studying condensed phase dynamics of chromophore molecules has been continuously demonstrated.14−16 In this work, where we aim to elucidate the effect of model consistencies for describing a chromophore system, we will adopt the coumarin 153 molecule with an already constructed interpolation data set.14 (See Figure 1 for the system depiction.) We have chosen coumarin 153 here because it has been widely used for studying polar solvation dynamics due to its relatively large reorganization effect17,18 and because its interpolated potential has been verified against experiment.14 Namely, the time scale and other behaviors of excited state relaxation dynamics from simulations

Figure 1. Typical shapshot of the chromophore in methanol solvent. For visual clarity, only one methanol molecule is shown with thick cylindrical representations. 2876

DOI: 10.1021/acs.jpclett.5b01141 J. Phys. Chem. Lett. 2015, 6, 2875−2880

Letter

The Journal of Physical Chemistry Letters trajectory, and its autocorrelation function C(t) was calculated as usual C(tn) =

1 N

the ground-state trajectory simulations (Stage 1) and because the solvent descriptions with PES-ig and PES-1 are the same, the contrast must have been caused by the chromophore vibrations. Indeed, PES-i includes realistic anharmonicity in the surface description. With the anharmonicity, the vibrational normal modes of the chromophore are strongly coupled to each other, and the excess energy of any vibrational mode will quickly dissipate into the chromophore’s own bath of vibrations. On the contrary, PES-1 possesses many simple harmonic functions in the model, and this had led to the overly oscillatory behavior for the scheme with PES-1. Considering that a harmonic oscillator stays on a single orbit in the phase space, it is not surprising to observe similarly regular behaviors in the time correlation with the mostly harmonic PES-1 potential. It also suggests that with PES-1, if a vibrational mode is excited along with electronic excitation, the bath will not efficiently quench that vibrational energy. Similarly, when PES2 and PES-3 instead of PES-1 were adopted during the first MD simulation stage, the time correlations displayed almost identical results except some changes in oscillation frequencies caused by the switches in force constants of bond stretches and angle bendings (Figure S1 in Supporting Information (SI)). After all, we used the same (PES-ig, PES-ie) surface pair for the gap energy evaluations in all of these cases of the time correlation functions, and our comparisons show that the sampling process during the first stage can still affect how the system−bath interactions will be represented during the second stage. Now let us observe the changes in the spectral densities. Of course, the nature of the gap energy autocorrelation will directly affect the overall shapes of the simulated spectral densities, as can be easily expected from the mathematical expression of eq 2. Because of fast decaying oscillations, the most realistic PES-ig/(PES-ig, PES-ie) scheme generates broad peaks in the density (Figure 3a), while PES-1/(PES-ig, PES-ie) with the long-lasting oscillations in its time-correlation produces sharp peaks in J(ω) (Figure 3b). In fact, the sharp peaks from the PES-1/(PES-ig, PES-ie) scheme result from strong vibrational components that are oscillating at the normal-mode frequencies of the PES employed during the MD simulation stage. To confirm this, we also inspected the spectral density from the PES-ig/(PES-1, PES-M) scheme and observed well-overlapping peaks with those from the PES-ig/ (PES-ig, PES-ie) scheme (Figure S2a in the SI). This is simply because the two schemes adopted the same PES for the trajectory sampling during the first stage. Similarly, the peaks in spectral densities from PES-1/(PES-1, PES-M) and PES-1/ (PES-ig, PES-ie) schemes (Figure S2b in the SI) also agree well in their shapes and peak positions. We will now present more detailed comparisons of the results obtained with the schemes with the three MM PESs during the first sampling stage. In our notation, these correspond to PES-1/(PES-ig, PES-ie), PES-2/(PES-ig, PESie), and PES-3/(PES-ig, PES-ie), and they will allow us to investigate how the varying degrees of discrepancies between the PES models during the two separate computational stages will affect the appearance of the system−bath interaction. Figure 3b−d displays the spectral densities obtained for these cases. Let us first compare locations of the sharp and intense peaks in the region above 1000 cm−1. As previously mentioned, their frequencies are closely related to the relevant PES normalmode frequencies. The spectral densities show large similarity between PES-2/(PES-ig, PES-ie) and PES-3/(PES-ig, PES-ie)

N

∑ [ΔE(tk) − ΔE̅ ][ΔE(tk + tn) − ΔE̅ ] k=0

(1)

Here N is the number of windows averaged, ΔE(t) is the excitation energy at time t, and ΔE̅ is the average of the excitation energies over the entire trajectory. The value of N was fixed for any t to make all of the points in the correlation function statistically equivalent. This time correlation can be converted into spectral density by J(ω) =

2 β ℏω πℏ 2

∫0

∞

dt C(t ) cos ωt

(2)

Here, for the prefactor, we adopted βℏω/2 instead of traditional tanh(βℏω/2), following the suggestion by AspuruGuzik and coworkers,29 as the time correlation functions were obtained with classical simulations. It was also shown that the temperature-independent nature of the spectral density can be properly reproduced with this expression.29,30 As previously described, we employed a series of different potential surfaces for MD simulations: PES-ig, PES-1, PES-2, and PES-3. During the second stage of the gap energy calculations, the same set of surfaces can still be adopted. We will adopt a notation of “PESl/(PES-m, PES-n)” when the ground-state MD trajectory was simulated with PES-l, and the gap energy was evaluated between PES-m and PES-n. For example, for a scheme where the trajectory is simulated with PES-1 and the gap energy is calculated with interpolation, it will be denoted as PES-1/(PESig, PES-ie). For comparison purposes, in some situations, we also adopted an MM-style excited-state surface created by Hessian matching process14 during the second stage of calculating the gap energies. We will denote this MM-based gap energy evaluation with PES-M, and then the MM-only approach will be denoted as, for example, PES-1/(PES-1, PESM) in our notation. Figure 2 displays the time correlation functions obtained with PES-ig/(PES-ig, PES-ie) and PES-1/(PES-ig, PES-ie) approaches. One can clearly see that the interpolation-only scheme of PES-ig/(PES-ig, PES-ie) is accompanied by only short-lived oscillations, while the PES-1/(PES-ig, PES-ie) scheme has led to long-lasting oscillations. Because the difference between these two schemes is induced only from

Figure 2. Gap energy autocorrelation functions from PES-ig/(PES-ig, PES-ie) (red) and PES-1/(PES-ig, PES-ie) (black) schemes, obtained from MD simulations of coumarin 153 in methanol. The autocorrelation function from the mixed PES-1/(PES-ig, PES-ie) scheme exhibits strong oscillations due to the excessive harmonicities in PES-1. 2877

DOI: 10.1021/acs.jpclett.5b01141 J. Phys. Chem. Lett. 2015, 6, 2875−2880

Letter

The Journal of Physical Chemistry Letters

by torsional degrees of freedom, which are described identically through all MM models in our schemes. These aspects well explain why the overall spectral shapes in the low-frequency limit appear strikingly similar when the trajectory simulations were performed with PES-1 through PES-3 models (Figure 3). Indeed, when we tried narrowing the chromophore torsional motions by imposing relevant restraints, the spectral densities in this frequency region diminished to a drastic extent (Figure S3 in the SI). Moreover, the spectral densities below 200 cm−1 from all simulations with such restraints were almost identical to each other, and they could be well explained with only the chromophore-solvent nonbonded interactions (Figure S3 in the SI). With this argument and the subsequent observations, one should expect that the spectral density from the interpolationonly scheme (PES-ig/(PES-ig, PES-ie)) with completely different torsional potential may exhibit radical differences, and this indeed is the case for our results shown in Figure 4.

Figure 4. Direct comparisons of the spectral densities from the various schemes in the low-frequency regime.

Phenomenologically with this Figure, we also note that the MM-based models significantly exaggerate the contributions from the low-frequency end region below 100 cm−1. This may be important as coupling in the lower frequency region has a relatively stronger effect on the reorganization event after electronic excitation/de-excitation, as will be explained in the next part. The degree of relaxation upon electronic transition can be represented by the reorganization energy, which is related to the spectral density as

Figure 3. Spectral densities obtained with (a) PES-ig/(PES-ig, PES-ie), (b) PES-1/(PES-ig, PES-ie), (c) PES-2/(PES-ig, PES-ie), and (d) PES3/(PES-ig, PES-ie) schemes. Sharp peaks emerge in the high-frequency region of spectral densities in panels b−d with MM models during the MD stage. For visual clarity, sharp peaks are shown only in the lower part with J < 5 kJ2/mol2.

cases, while the result from PES-1/(PES-ig, PES-ie) shows notable differences from these two. This is due to the accord/ disaccord in the force constant parameters, which directly affect the normal-mode frequencies. Recall that the parameters for PES-1 were generated from chromophore-specific Hessian matching,14 while those for PES-2 and PES-3 were taken from the GAFF set. Although PES-2 and PES-3 differ by equilibrium bond lengths and angles, these geometrical disagreements only resulted in the changes of peak heights, which is seemingly a minor effect. Meanwhile, the low-frequency region (below 1000 cm−1) with rather complex structures is quite similar among these three schemes regardless of the differences in the MM models. We ascribed this similarity to the fact that low-frequency modes tend to involve grouped motions of the molecule. Because the differences in individual stretches and bendings will likely average out, they will be less affected by the model differences. More importantly, the grouped motions will be heavily affected

J(ω) dω (3) ω This expression becomes exact only when the shape of groundand excited-state surfaces are identical, except some shift along the phonon mode coordinates.10 Even in cases where this condition is not met, it is still useful for estimating the magnitude of the reorganization energy.31,32 The reorganization energies computed from the various schemes are listed in Table 1. The coumarin molecule undergoes a distinct relaxation after electronic transition,14 and λ from the reference interpolation-only scheme is relatively large as ∼15 400 cm−1.33 Interestingly, one can notice that the three schemes that adopted MM models during the first MD sampling stage generate very similar results, which are somewhat off from the number from the interpolation-only scheme. It should be stressed again that these three schemes with MM sampling are using the interpolated surfaces during λ≈

2878

∫

DOI: 10.1021/acs.jpclett.5b01141 J. Phys. Chem. Lett. 2015, 6, 2875−2880

Letter

The Journal of Physical Chemistry Letters

manner toward excited-state dynamics simulations will inevitably create some errors due to the discrepancy between the models employed during the MD sampling and the excitation energy calculations. The mixing scheme is for the trade-off between efficiency and accuracy, and investigating the significance of the errors has not been possible in the practical sense. In this work, by employing the potential energy surface interpolation technique, which can closely mimic the reliability of high-level QC calculations at reasonably low expenses,13,16 we have examined the effect of the model discrepancy by focusing on how the chromophore vibrations affect the electronic gap energy correlations with the example of the coumarin 153 chromophore in solution. Specifically, we have evaluated the gap energy time correlation functions by adopting various combinations of potential energy surfaces for the sampling and the gap energy calculation stages and have converted them to spectral densities and reorganization energies. Not surprisingly, when different models were used, detailed features of the system−bath interactions appeared differently. The characteristics of the potential energy surface employed during the MD sampling stage leave their signatures in the gap energy time correlation function and the related spectral density, even when the potential is not used at all during the separate gap energy calculation stage. This should be expected as the chromophore vibrations will be driven by the potential engaged during the sampling stage; however, the effect of the model discrepancy on the overall reorganization energy (λ) was found to be relatively small, suggesting that λdriven processes will not be seriously affected by using the inherently approximate mixed schemes; however, cautions should still be taken as the detailed shapes of the spectral density may turn out quite differently when approximate MM models are used during the sampling stage. Indeed, there have been reports showing that the shape of the spectral density function can actually affect the excited-state dynamics.34 In addition, our results show that even the low-frequency phonon modes are strongly affected by the chromophore model. Thus, when one attempts to handle detailed behaviors of exciton− phonon coupling with chromophores similar to ours or even larger, one will need to utilize a scheme where the potential surfaces throughout all stages of simulations are consistent to each other. Of course, a drawback of achieving consistency will be the potentially unmanageable computational cost. We anticipate that recent advances in the approaches of constructing potential surfaces with both accuracy and efficiency will be of great utility in this area.

Table 1. Reorganization Energies, Numerically Calculated from the Spectral Densities entry

scheme notation

reorganization energy (cm−1)

1 2 3 4 5

PES-ig/(PES-ig, PES-ie) PES-1/(PES-ig, PES-ie) PES-2/(PES-ig, PES-ie) PES-3/(PES-ig, PES-ie) PES-1/(PES-1, PES-M)

15415 (7606)a 17289 (10101)a 17018 (9793)a 17059 (9485)a 8807 (3723)a

a

Numbers in parentheses denote partial reorganization energies obtained by integrating eq 3 up to the 100 cm−1 range.

the second gap energy evaluation stage. It is also interesting to see that the λ values in the table with the mixed MM/ interpolation strategies (∼17 000 cm−1) differ from what we obtained with the interpolation-only scheme (PES-ig/(PES-ig, PES-ie)) by only ∼10%. In contrast, the λ value from the MMonly scheme (PES-1/(PES-1, PES-M)) completely deviates from all others. These aspects suggest that adopting the MM model during the first stage of trajectory sampling simulations will not likely be a too drastic approximation, while adopting it for gap energy evaluations may lead to serious errors. Simply, the precision of MM models is not sufficient for the gap energy evaluations, as we mentioned in a previous part of this paper. In addition, even though the use of an MM model during the sampling stage does not appear to affect the reorganization energy seriously, it still does not provide a full justification of the mixed scheme, as the reorganization energy is in a sense an oversimplified measure. In particular, it does not bear any frequency-dependent information. Nevertheless, physical events that are strongly affected by the magnitude of the reorganization energy but not by the detailed shape of the spectral density will likely be influenced to a limited extent by using MM potentials toward sampling, which has been the only amenable approach in the practical sense until very recently. The trends with which the reorganization energies change with different schemes can actually be easily understood from the shapes of the spectral densities. As can be seen in eq 3, the spectral density affects λ in inverse proportion to ω, and thus the contribution to the reorganization is amplified at a lower frequency. Indeed, the three mixed schemes in Table 1 with similar reorganization energies (entry nos. 2−4 in the Table) have displayed closely overlapping spectral densities in Figure 4. In comparison, with the interpolation-only scheme, the lowfrequency density is significantly smaller, leading to the ∼10% decrease in λ. In fact, when the numerical integration of eq 3 is conducted only in the low-frequency domain under 100 cm−1, the result with the interpolation-only scheme differs more from the ones with other mixed schemes as shown with the numbers in parentheses in Table 1. Considering again that we have adopted the same solvent model for both the interpolation-only and the MM + interpolation mixed schemes, we should note that the phonon modes in the low-frequency limit are actually strongly influenced by how we model the chromophore itself, not just by the way we model the surrounding solvent. In addition, because the phonon modes in the region below 100 cm−1 contribute to ∼50% or even more to the overall reorganization, at least in this case we tested, accurately describing chromophore vibrations will be crucial for quantitatively describing relaxation effects through exciton− phonon coupling. In conclusion, we should mention that a widely adopted protocol of modeling system−bath interaction in a mixed

■

ASSOCIATED CONTENT

S Supporting Information *

Brief description on the interpolation approach; additional gap energy time correlations and spectral densities; and importance of the chromophore torsional motions toward the spectral density in the low-frequency regime. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.5b01141.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. 2879

DOI: 10.1021/acs.jpclett.5b01141 J. Phys. Chem. Lett. 2015, 6, 2875−2880

Letter

The Journal of Physical Chemistry Letters

■

(20) Casanova, D.; Rhee, Y. M.; Head-Gordon, M. Quasidegenerate scaled opposite spin second order perturbation corrections to single excitation configuration interaction. J. Chem. Phys. 2008, 128, 164106. (21) Vladimirov, E.; Ivanova, A.; Rösch, N. Effect of Solvent Polarization on the Reorganization Energy of Electron Transfer from Molecular Dynamics Simulations. J. Chem. Phys. 2008, 129, 194515. (22) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (23) Song, C.-I.; Rhee, Y. M. Development of Force Field Parameters for Oxyluciferin on its Electronic Ground and Excited States. Int. J. Quantum Chem. 2011, 111, 4091−4105. (24) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (25) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (26) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (27) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (28) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (29) 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. (30) 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. (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) Kim, H. W.; Kelly, A.; Park, J. W.; Rhee, Y. M. All-Atom Semiclassical Dynamics Study of Quantum Coherence in Photosynthetic Fenna-Matthews-Olson Complex. J. Am. Chem. Soc. 2012, 134, 11640−11651. (33) The experimental reorganization energy is 2885 cm−1, as reported in Mukerjee, P.; Halder, M.; Hargrove, M. S.; Petrich, J. W. Characterization of the Interactions of Fluorescent Probes with Proteins: Coumarin 153 and 1,8-ANS in Complex with Holo- and Apomyoglobin. Photochem. Photobiol. 2006, 82, 1586−1590. As explained in the introductory part, the large deviation originates from the limitation of the adopted QC method, SOS-CIS(D0), and from the lack of polarization in the model. If we assume that the interpolated potential is the reference surface that the MM model is trying to mimic as closely as possible, this lack of agreement does not pose a problem to our analysis. (34) Mohseni, M.; Shabani, A.; Lloyd, S.; Rabitz, H. Energy-Scales Convergence for Optimal and Robust Quantum Transport in Photosynthetic Complexes. J. Chem. Phys. 2014, 140, 035102.

ACKNOWLEDGMENTS This work has been supported by Institute for Basic Science (IBS) of Korea. C.W.K. was a recipient of Global Ph.D. fellowship awarded by National Research Foundation (NRF) of Korea.

■

REFERENCES

(1) Fleming, G. R.; Cho, M. Chromophore-Solvent Dynamics. Annu. Rev. Phys. Chem. 1996, 47, 109−134. (2) Yang, M.; Fleming, G. R. Influence of Phonons on Exciton Transfer Dynamics: Comparison of the Redfield, Fö rster, and Modified Redfield Equations. Chem. Phys. 2002, 275, 355−372. (3) Heřman, P.; Kleinekathöfer, U.; Barvík, I.; Schreiber, M. Influence of Static and Dynamic Disorder on the Anisotropy of Emission in the Ring Antenna Subunits of Purple Bacteria Photosynthetic Systems. Chem. Phys. 2002, 275, 1−13. (4) Shibata, Y.; Nishi, S.; Kawakami, K.; Shen, J.-R.; Renger, T. Photosystem II Does Not Possess a Simple Excitation Energy Funnel: Time-Resolved Fluorescence Spectroscopy Meets Theory. J. Am. Chem. Soc. 2013, 135, 6903−6914. (5) Aghtar, M.; Liebers, J.; Strümpfer, J.; Schulten, K.; Kleinekathöfer, U. Juxtaposing Density Matrix and Classical Path-Based Wave Packet Dynamics. J. Chem. Phys. 2012, 136, 214101. (6) Fidy, J.; Laberge, M.; Kaposi, A. D.; Vanderkooi, J. M. Fluorescence Line Narrowing Applied to the Study of Proteins. Biochim. Biophys. Acta, Protein Struct. Mol. Enzymol. 1998, 1386, 331− 351. (7) 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. (8) Shim, S.; Rebentrost, P.; Valleau, S.; Aspuru-Guzik, A. Atomistic Study of the Long-Lived Quantum Coherences in the FennaMatthews-Olson Complex. Biophys. J. 2012, 102, 649−660. (9) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (10) Nitzan, A. Chemical Dynamics in Condensed Phases; Oxford University Press: New York, 2006. (11) Cheng, Y. C.; Fleming, G. R. Dynamics of Light Harvesting in Photosynthesis. Annu. Rev. Phys. Chem. 2009, 60, 241−262. (12) Johnson, J. C.; Nozik, A. J.; Michl, J. The Role of Chromophore Coupling in Singlet Fission. Acc. Chem. Res. 2013, 46, 1290−1299. (13) Ischtwan, J.; Collins, M. A. Molecular Potential Energy Surfaces by Interpolation. J. Chem. Phys. 1994, 100, 8080. (14) Park, J. W.; Kim, H. W.; Song, C.-I.; Rhee, Y. M. Condensed Phase Molecular Dynamics Using Interpolated Potential Energy Surfaces with Application to the Resolvation Process of Coumarin 153. J. Chem. Phys. 2011, 135, 014107. (15) Park, J. W.; Rhee, Y. M. Interpolated Mechanics−Molecular Mechanics Study of Internal Rotation Dynamics of the Chromophore Unit in Blue Fluorescent Protein and Its Variants. J. Phys. Chem. B 2012, 116, 11137−11147. (16) Park, J. W.; Rhee, Y. M. Towards the Realization of Ab Initio Dynamics at the Speed of Molecular Mechanics: Simulations with Interpolated Diabatic Hamiltonian. ChemPhysChem 2014, 15, 3183− 3193. (17) Eom, I.; Joo, T. Polar Solvation Dynamics of Coumarin 153 by Ultrafast Time-Resolved Fluorescence. J. Chem. Phys. 2009, 131, 244507. (18) Maroncelli, M.; Fleming, G. R. Picosecond Solvation Dynamics of Coumarin 153: The Importance of Molecular Aspects of Solvation. J. Chem. Phys. 1987, 86, 6221. (19) Jung, Y.; Lochan, R. C.; Dutoi, A. D.; Head-Gordon, M. Scaled opposite-spin second order Møller−Plesset correlation energy: An economical electronic structure method. J. Chem. Phys. 2004, 121, 9793−9802. 2880

DOI: 10.1021/acs.jpclett.5b01141 J. Phys. Chem. Lett. 2015, 6, 2875−2880