On the alternatives for bath correlators and spectral densities from mixed quantum-classical simulations

We investigate on the procedure of extracting a"spectral density"from mixed QM/MM calculations and employing it in open quantum systems models. In particular, we study the connection between the energy gap correlation function extracted from ground state QM/MM and the bath spectral density used as input in open quantum system approaches. We introduce a simple model which can give intuition on when the ground state QM/MM propagation will give the correct energy gap. We also discuss the role of higher order correlators of the energy-gap fluctuations which can provide useful information on the bath. Further, various semiclassical corrections to the spectral density, are applied and investigated. Finally, we apply our considerations to the photosynthetic Fenna-Matthews-Olson complex. For this system, our results suggest the use of the Harmonic prefactor for the spectral density rather than the Standard one, which was employed in the simulations of the system carried out to date.


I. INTRODUCTION
In the study of the dynamics of large systems such as photosynthetic complexes, reduced models, which provide information on a small set of system degrees of freedom at the price of tracing out the rest of the bath degrees of freedom, have become very popular.Amongst these methods, which are open quantum systems approaches, one can find various Quantum Master Equations  and Stochastic Schrödinger Equations [23,24], which often rely on describing the system-bath interaction through a two-time bath correlation function or a bath spectral density.Therefore, it is of considerable interest to obtain these quantities.While a full quantum mechanical treatment of such large systems is out of reach, one viable approach is to use a mixed quantum-classical approach for the nuclear-electronic degrees of freedom.However, there is no unique way of obtaining a bath correlation function or a bath spectral density when resorting to quasiclassical theories.In this work, we provide criteria that can be helpful to choose an appropriate strategy for this task.
As a case study, we consider the Fenna-Matthews-Olson (FMO) light-harvesting pigment protein complex found in green sulfur bacteria.For this system, recent efforts have been undertaken to extract the bath spectral densities from mixed quantum-classical calculations [25][26][27].The FMO complex has a trimeric structure, where each monomer contains, within the protein scaffolding, eight bacteriochlorophyll (BChl) molecules, which can transport electronic excitation energy.Up to recently, it was thought that only seven of the BChl's actually were present and most of the previous studies have focused on that case.Shim's results [26] which we employ in this work, are indeed based on the case of one monomer with seven BChl molecules.Experimentally, it has been possible to extract a spectral density for the BChl with the lowest transition energy [28,29].However, one can expect that each BChl has a different spectral density, due to its specific protein environment.* aspuru@chemistry.harvard.eduOne theoretical approach to obtain the spectral densities from a microscopic description is a mixed quantum mechanics/classical mechanics (QM/MM) model [30].In this approach, the nuclear degrees of freedom are treated classically and the relevant system quantities are calculated quantum mechanically.Then, from the microscopic description, spectral densities and correlation functions can be extracted and employed in the reduced models.
A specific QM/MM approach, which has become popular in recent years in the context of photosynthetic complexes [25,26,31,32] and has been employed for FMO [25,26], consists in propagating the nuclei in the ground electronic state of the FMO complex, thus the change in the classical forces due to excitation of the BChls is ignored.The bath correlation function and spectral densities are then extracted from the energy gap trajectories, i.e. the electronic transition energies which depend on the time dependent nuclear configuration.This transition energy is calculated using quantum chemistry, for example TDFT [26] or semi-empirical approaches [25].One thus obtains a time dependent energy gap two-time correlation function.Usually, a spectral density (SD) is derived from the time correlation function, to characterize the frequency dependent coupling of the electronic transitions to the environmental degrees of freedom.In the previous investigations on the FMO complex [25][26][27], the spectral densities differ by orders of magnitude respect to each other and also with respect to the SD extracted from experiment [28,29].
In this work we revisit the data of Shim [26].We shed light on the connection between the mixed QM/MM gap correlation function and the open quantum system bath correlation function using a simple model.The mixed QM/MM gap correlation function is real.However, in general, the full quantum correlation function will have an imaginary part.We employ different semiclassical a posteriori corrections to recover this part and compare the resulting spectral densities.Much work has been carried out on these a posteriori semiclassical corrections [33][34][35][36][37][38], but the question of which approximation is best remains open.Towards answering this question, we show that a simple model of shifted harmonic Born-Oppenheimer surfaces leads to two of the semiclassical a posteriori correc-tions, each obtained with a different phase space probability distribution.We thus establish the link to a microscopic picture.This model of shifted harmonic potential surfaces is of particular interest, since the spectral density used in the open quantum system approaches emerges from such a description.
Finally, we will investigate whether the results of the QM/MM fulfill the requirements needed to employ the extracted spectral density in open quantum system methods.These methods rely on the validity of assumptions such as linear coupling of the system to the bath and a bath of harmonic oscillators.We attempt to investigate whether these assumptions are valid in our case by evaluating higher order correlators of the energy gap time traces.In addition, we compare the spectral densities obtained at different temperatures.In most open quantum system approaches, when the harmonic bath approximation is employed, the spectral density is temperature independent.Thus, one can use this invariance as a criteria for choosing which of the applied a posteriori corrections is most reasonable to be employed in these methods.In particular, our findings suggest that the best a posteriori semiclassical approximation for FMO is the Harmonic [39,40] correction rather than the Standard [41][42][43] one, which has so far been employed in the context of the simulation of exciton dynamics in photosynthetic complexes [25][26][27]32].Together, these aspects provide a clearer microscopic picture of the complex approximations involved in combining ground state QM/MM and open quantum system approaches.
The paper is structured as follows: we begin by introducing the general quantum two-time correlation function in Section II.We introduce its time symmetries and its Fourier transform and subsequently we define the spectral density.A brief summary of the general a posteriori semiclassical approximations to the quantum Fourier transform of the correlator from the classical Fourier transform is given in Section II C. In Section III, we introduce the concept of an energy gap correlation function for two-level systems as models for molecules coupled to a bath and show how this leads to a quantum bath correlation function and spectral density which are consistent with the open quantum system approach.In Section IV, we show that one can introduce a microscopic model which leads to some of the same prefactors described in the general case in Section II C. Finally, we investigate the conditions of linear system-bath coupling and harmonic bath in Section V.In particular, we evaluate high-order multi-time correlation functions for the bath.These considerations are applied to our specific QM/MM calculations for FMO in Section VI.We conclude in Section VII by summarizing our findings.

II. THE QUANTUM CORRELATION FUNCTION AND THE SPECTRAL DENSITY
In this section, we introduce the definition of the quantum two-time bath correlation function.The generic Hamiltonian of a system coupled to a bath, in the absence of external fields, can be expressed as Ĥ = ĤS (q,p) + ĤB (Q,P) + ĤSB (q,p, Q, P) , where ĤS is the system Hamiltonian, ĤB is the bath Hamiltonian, ĤSB is the system-bath Hamiltonian.In addition, (q,p) = (q j , p j ) and (Q,P) = (Q k , P k ), indicate the generalized multidimensional conjugated coordinates for the system and the bath respectively.The indexes j = 1, ..., f and k = 1, ..., F run over the system (f ) and bath (F ) degrees of freedom respectively.The system-bath Hamiltonian can be written as a function of the system, Â, and bath, B, operators: The influence of the bath on the system can be described by time-correlation functions.We will mostly focus on the two-time bath correlation function Here, Bm (t, Q, P) = e i ĤBt/ Bm (Q, P)e −i ĤBt/ , and where β = 1/ (k B T ) and T is the temperature.In the following, we will be interested only in the n = m correlators, which we will indicate as C(τ ) with τ = t − t , dropping the subscript notation for simplicity.In section V, we will briefly discuss higher order correlators.The correlator defined above is in general complex and one can show, see e.g.[41,44], that it has the following symmetries with respect to time, A. Fourier transform of the time correlation function and symmetries of the correlator We define G(ω), the Fourier transform of the time correlation function The function G(ω) is in general, temperature-dependent, real and positive.In this work, we will refer to it as the Temperature-Dependent Coupling Density (TDCD).It will be convenient to split G(ω) into a symmetric and antisymmetric component which originate respectively from the real and imaginary parts of C(t), In this definition, we have followed the convention of Ref. [34].Note that in the literature there exist other definitions, e.g. the corresponding equations in Ref. [45], differ by a factor of 2 from the ones used here [46].The detailed-balance condition, which follows directly from the second time symmetry in Eq. 5, implies that the overall TDCD is related to its asymmetric [47] part by = (1 + coth (β ω/2)) G asym (ω).
It will be convenient to abbreviate G asym (ω) by defining Using Eq. 9 and the definition of G(ω), Eq. 6, one can express the correlation function as a function of J(ω),

B. The spectral density
Another quantity which is often of interest is the so-called "spectral density".The spectral density describes the frequency dependent coupling of the system to the bath.There are different definitions of spectral density in the literature (for example J(ω) is sometimes refered to as the spectral density).We follow the convention of defining the spectral density as a positive frequency function Here Θ(ω) is the Heavyside function, which is one for positive arguments and zero for negative ones.The scaling by π has been introduced for later convenience.Note that C. General semiclassical a posteriori approximations For systems of more than a few degrees of freedom, and in general, it is difficult to calculate the exact correlation function, and therefore its Fourier transform, by using a fully quantum mechanical treatment.However, using classical mechanics one can obtain its classical counterpart with much less effort.Therefore, it is common to attempt to construct the quantum spectral density from the classical one.
We define the fully-classical correlation function as the classical → 0 limit of Eq. 3, Here W (Q, P) is the classical bath phase-space density, defined as and the quantum bath operators B in Eq. 3 have been substituted by classical functions of the phase space variables B (t, Q, P).
The classical TDCD is defined as Note that C cl (t) is a real and symmetric function in contrast to its quantum counterpart.This is also the case in the mixed QM/MM simulations employed for FMO [25,26].The QM/MM correlation function obtained is real and no information about the important imaginary part of the quantum correlator is available a priori.
It is now desirable to be able re-construct, at least partially, the exact quantum spectral density from the classical one, through a simple description.Ideally, such a correction should be applied a posteriori and should not require extensive additional computation.Much work has been carried out in this direction, see e.g.[33][34][35][36][37][38].As described in Ref. [34], one can define various semiclassical approximations to the full quantum mechanical G(ω) starting from its classical counterpart G cl (ω).We report each of these approximations in Table I, second column.
These corrections all originate from expansions in and use of the symmetry properties of the two-time correlation function and its Fourier transform.Note that if one expands the quantum correlator C(t) in powers of , the first term is real and symmetric and corresponds to C cl (t).The assumption that C(t) = C cl (t), which leads to the standard approximation, is in general not correct.In fact, since both of the correlation functions are obtained after thermal averaging, we see that they must differ at least by their respective partition functions.
At low frequencies, ωβ ≡ ω b < 1 (i.e.ω < k B T ) all approximations give nearly identical results and give the same value for ω b = 0.
The various approximations for J(ω), and thus for the spectral density, can straightforwardly be derived from those of G(ω) by using Eq. 9.The resulting expressions are reported in column three of Table I and the prefactors follow the same trend as those for G(ω) as a function of frequency.Now, given all the functional forms described above, the question is how to choose the most appropriate one.For the FMO complex, it is unclear at first sight which one would be the best.In Section IV, we will investigate a model to elucidate the origin of these prefactors.This will help to discriminate between these corrections.In Section VI, we will apply all of the corrections listed in Table I to our energy gap traces and discuss the differences between each approach.

III. ENERGY GAP CORRELATION FUNCTION FOR A SIMPLE MODEL
In the mixed QM/MM calculations for photosynthetic systems [25,26,32], the nuclear trajectories are propagated in the electronic ground state using MD with short time steps.For Table I.Column two: Various expressions for obtaining a semiclassical temperature-dependent coupling density TDCD G(ω) from the classical G cl (ω) as discussed in, e.g.[34].Column three: Expressions for obtaining the semiclassical asymmetric TDCD J(ω) from the classical G cl (ω).These follow from the expressions in column two and from detailed balance (Eq.9).

Method
Expression for G(ω) Expression for J(ω) = Gasym(ω) a set of longer times steps within these trajectories, the electronic transition energies of the BChl molecules are computed using an electronic structure calculation method.Because it is computationally costly to calculate the electronic states for the full set of seven/eight coupled BChls simultaneously [25], the system was divided into seven/eight subsystems for which the electronic states were calculated separately.Thus, in these calculations no excited state interactions are included explicitly [50].The Hamiltonian of the coupled BChls is then written as H = N n=1 H n + n<m V nm where H n denotes the Hamiltonian of BChl n and V nm is the Coloumb (transition dipole-dipole) interaction between them.To establish a connection to the open quantum system approach, each BChl is treated as an electronic two level system.These two-level systems and the electronic interaction between them are taken to be the system part.The coupling to internal nuclear degrees of freedom and the surrounding protein will then lead to fluctuations of these quantities in time (for more details see e.g.[32]).From the time dependence of the transition energy between electronic ground and excited state for each BChl, a classical ground-excited state energy-gap correlation function can be obtained.In turn, spectral densities can be extracted from the energy-gap correlation functions.
The gap correlation function, as obtained from the MD simulations, is a quantity which up to the previous section, has not been connected to the open quantum system approach described in Sec.II.In this Section, we will explore a simple model with Born-Oppenheimer (BO) surfaces which can clarify the connection.

A. Quantum correlation function and energy gap correlation function for a molecule
Lets us begin by considering a single molecule (BChl) treated in the Born-Oppenheimer approximation.
The molecule is modeled as a two-level system with an electronic adiabatic ground |g and excited |e state.We can think of the BO-surfaces as having the dependence of the environment (protein and other BChls) already included, ignoring however the resonant dipole-dipole interaction.The approximation of two levels is reasonable in the limit where the next excited state is very far in energy space from the first.Usually, nonadiabatic couplings can be also neglected, as chosen in our calculations.
Given this model, we investigate how the general correlation function, Eq. 3, is related to the energy gap correlation function.
We write the full Hamiltonian formally as where Ĥg (Q) and Ĥe (Q) are the nuclear Hamiltonians for the ground and excited state in the BO approximation.
In mass scaled coordinates Q j = √ m j q j ; P j = p j / √ m j , the Hamiltonians can be expressed as Ĥg (Q,P) = denotes the grounds state potential energy surface.For later purpose, we have expressed the excited state nuclear Hamiltonian with respect to the ground state potential by introducing the energy gap operator, This operator quantifies the energy difference between the excited state and the ground state surface.A coordinate independent constant energy difference ω eg + λ 0 has been explicitly written down, so that the remaining part V e (Q) − V g (Q) does not contain any coordinate independent contributions.This division and the meaning of ω eg and λ 0 will become clear in Section III B.
The total Hamiltonian can be rewritten as where we have defined the reduced gap operator ∆ ≡ ∆eg − ω eg − λ 0 .
To establish a connection to the open quantum system model, as presented in Sec.II, we choose where we have set the energy of the electronic ground state |g to zero.From the form of ĤSB we identify the system operator Âe = |e e| and the bath operator B = ∆eg (Q).We can now define the usual bath correlation function as where we have dropped the dependence on bath coordinates in the notation for simplicity.∆ can be thought of as a "gap" operator, that is, as a measure of the energy difference between the ground and excited state at a given nuclear configuration.
From now on we will indicate reduced gap correlation functions as to distinguish them from the general bath correlation function C(t).Eq. 24 corresponds to the full quantum gap correlation function that one would obtain, e.g., from a quantum simulation on the FMO complex, considering only two electronic levels per molecule and after including the protein environment.

B. Quantum correlation function and energy gap correlation function for harmonic surfaces
While the approach outlined in the previous section is applicable to arbitrary potential surfaces, in most of the open quantum system approaches used to describe the FMO complex, the bath is taken as an (infinite) set of harmonic oscillators for the environment of each BChl.Each oscillator coordinate is then assumed to be linearly coupled to the electronic excitation of the BChls, i.e.H SB = |e e| ⊗ j κj Q j where κj is a coupling constant.
To establish the connection between the reduced gap operator and this system-bath interaction, we now consider identical shifted harmonic potential surfaces, as sketched in Fig. 1.The nuclear Hamiltonians defined in the general case in the previous Section III A become, Ĥg (Q where Ω j is the frequency of the j-th oscillator.This model for a finite small number of vibrational modes of the chromophores, has been successfully employed to describe the optical properties of molecular aggregates [51][52][53][54].These Hamiltonians can be rewritten as function of a † j and a j , the ground state bosonic creation and annihilation operators which are related to the conjugated coordinates by Here, the constant shift λ 0 , previously introduced in Eq. 18, corresponds to the frequently empoyed reorganization energy We have also introduced the so-called Huang-Rhys factor [55]: X j =  j and a (frequency dependent) coupling constant κ j = Ω j X j .Note that the total Hamiltonian is now in the standard form of an open quantum system model, as in Eq. 1, with the relevant quantities given in Table II.In particular the reduced energy gap operator is given by From this expression of the energy gap operator one obtains the quantum two-time bath correlation function (see e.g.[45]) with the temperature independent spectral density Note that from the definition Eq. 13 we have , which is also temperature independent.
To establish a connection to the classical correlator, which is real and symmetric, we note that j(ω) can be obtained from the real part of α(t) via When using this expression to obtain the spectral density from QM/MM simulations one often assumes that C cl (t) ≈ Re{α(t)}, following the Standard approximation.Then, after a Fourier transform and use of symmetry relations for G(ω) one finds the following expression, Table II.Expressions of the system bath quantities for the case of two Born-Oppenheimer harmonic surfaces as sketched in Fig. 1.

System-bath Hamiltonian
This is the expression (up to the constant prefactor 1/ ) used in Refs.[26,32], to obtain spectral densities.

IV. CLASSICAL AND SEMICLASSICAL LIMITS OF THE CORRELATORS AND SPECTRAL DENSITIES FOR HARMONIC SURFACES
As outlined in the previous section, the harmonic model allows for a simple analytic solution in the quantum mechanical case.Now we will show that the system also has a solution in the classical case.In particular, in this section, we will introduce a model to construct exact relations between the classical gap-correlation and the quantum one.To this end, we will consider classical dynamics in the ground state BO potentials within an initial value representation of the initial state which is consistent with the mixed QM/MM approach.For each initial value, we calculate a trajectory and the corresponding reduced classical energy gap between the two surfaces, i.e. ∆(Q(t), P(t)).We then average over many trajectories.

A. Classical equations of motion
The classical equation of motion of the j-th harmonic bath coordinate is Qj + Ω 2 j Q j = 0. Solving this differential equation with the initial condition (Q j0 , P j0 ) = (Q j (t = 0), P j (t = 0)) yields the time dependent coordinate trajectories For each trajectory, the energy gap is then given by where the parametric dependence of Q j and ∆ on the initial conditions (Q j0 , P j0 ) has been explicitly indicated.

B. Energy gap correlator
The evaluation of the reduced gap correlation function, Eq. 24, in the classical limit, results in the following expression where W(Q 0 , P 0 ) is the initial distribution and dP 0 dQ 0 denotes the set of all coordinates, i.e.
For harmonic potential surfaces, Eq 14, is time-evolved following Eq.(32).In this Section, we will investigate two different choices for the initial distribution, namely a Boltzmann distribution, as in Ref. 16, and a Wigner distribution which resembles the quantum thermal state.We will refer to the two cases as the classical limit and the semi-classical limit, respectively.
C. Classical and semiclassical correlation functions

Classical limit
To obtain the classical limit of the correlator, we choose the Boltzmann distribution for the initial coordinates which corresponds to a purely classical thermal state.The distribution is defined as follows with ) , and it is normalized to one, i.e., ´dP j0 dQ j0 W boltz j (Q j0 , P j0 ) = 1.Note that (Ω 2 j δQ j ) 2 = 2 X j Ω 3 j .Using Eq. 33 and the Boltzmann distribution for initial positions and momenta, we obtain, Here we have introduced the abbreviation ζ j ≡ Ω j / (k B T ).

Semiclassical limit
In order to obtain the semiclassical limit, we take the quantum Wigner distribution for the initial coordinates and use it in Eq. 33.The Wigner distribution is given by where we have used the compact notation The normalization of the Wigner distribution is chosen such that ´dPj0dQj0 2π W wig j (Q j0 , P j0 ) = 1.The resulting expression of the energy gap correlation function is

D. Classical and semiclassical spectral densities
After a Fourier transform of the classical correlators in Eq. 35 and 37 we obtain, for the Boltzmann distribution and for the Wigner distribution (39) Here κ j = Ω j X j as in Tab.II.Now, using Eq.28 for the spectral density j(ω) in the quantum case, and using J(ω) = π(j(ω) − j(−ω)) we can write By inverting these equations the exact quantum J(ω) can be expressed in terms of the classical We see that in our harmonic model the semiclassical Wigner distribution yields the same prefactor as for the Standard approximation described in Section II C while the Boltzmann distribution gives the same prefactor as the Harmonic approximation, also described in Section II C.

V. MODELS FOR SYSTEM-BATH COUPLING -HIGHER ORDER CORRELATORS
As discussed in the introduction, there has been a lot of interest in modeling the exciton dynamics of the FMO complex using open quantum system approaches.These usually require as input a bath two-time correlation function or (equivalently) a spectral density and they rely on the assumption of linear coupling to the bath and on a bath described by harmonic oscillators [56].
In the previous Section III B, we have discussed that this model corresponds to shifted adiabatic BO surfaces of identical curvature.We have shown that in this case, the energy gap two-time correlation function for a classical ground-state propagation is directly proportional to the quantum one and we have extracted the appropriate (frequency dependent) proportionality constant.For other shapes of the potential surfaces involved, one will in general obtain different proportionality constants, although the delta-peaks of the spectral densities can be located at the same energies (the positions are determined by the shape of the ground state potential).
It is not clear, a priori, if the approximation of shifted harmonic surfaces (or equivalently linear coupling to a harmonic bath) is a good one for the system under consideration.To gain some insight on this question, from an analysis of QM/MM trajectories, one possibility is to consider higher order correlators.If the approximation of linearly coupled harmonic oscillators is inadequate, one expects that higher order correlators will have a significant relative weight.
We proceed to discuss some properties of correlations of the bath gap operator, Eq. 18.The energy gap operators can be described by a function of the bath coordinates and expanded in terms of these as When only terms up to first order in Q are significant, as in the case of the Harmonic surfaces in the linear system bath coupling limit, Tab.II, we can write the two-time correlation function as Here, we have excluded the zeroth-order term which corresponds, e.g., to a reorganization energy, and is usually renormalized into the system Hamiltonian.The angular brackets ... = tr B {..., ρB } indicate thermal averaging over the bath degrees of freedom.Similarly, the three-time correlation function becomes i ξ (1) In the case of a harmonic bath, the three-time correlation function will vanish, and in general any odd permutation of the harmonic bath coordinates will vanish.
However, if one considers the case where one retains the second order term in Eq. 44, the two-time correlator will become: where we have defined Ξ ij and Q ij (t) as Analogously the three-time correlator becomes If the bath is harmonic, it is straightforward to show that all terms with an odd number of coordinate operators in the averages will vanish.Yet, we see that in general, unless the coupling to the bath coordinates is linear and the bath consists of Harmonic oscillators, the three-point correlator will not vanish.It may therefore be necessary to go beyond the simple description using only the two-time correlator.

VI. APPLICATION TO THE FMO COMPLEX
In this section, we apply the approximations discussed in Section II C, to the energy gap trajectories obtained from the mixed QM/MM simulations for the FMO complex of Prosthecochloris aestuarii as carried out recently by us in Ref. [26].The nuclear trajectories were obtained by classical MD using the AMBER 99 force field.An isothermal-isobaric (NPT) ensemble was employed in the MD simulations.For the calculation of the energy gap, snapshots of the nuclear coordinates were taken at every 4 fs.For each ground state configuration, the gap was obtained by computing the energy corresponding to the Q y transition of the BChl's using time-dependent time-dependent density functional theory with BLYP functional within the Tamm-Dancoff approximation.
The calculations were carried out at 77 and 300K and both temperature were treated on the same footing.We do not expect there to be additional sampling problems for the low temperatures because, up to current knowledge, FMO does not undergo any major conformational changes in this temperature range.More details on the computation can be found in Ref. [26].
The calculation of the SD from the time dependent gap energy is based on the model described in Section III.The actual MD simulation might deviate from this model e.g. because the thermostat could influence the dynamical evolution and thus the correlation function.We plan to investigate this aspect in future work.For now we will assume that the thermostat doesn't influence the dynamics and that the models introduced in Section III provide a reasonable description of a two level molecule treated in the QM/MM approach.
A. TDCD and spectral density from mixed QM/MM with a posteriori semiclassical corrections Using the energy gap trajectories obtained in Ref. [26], we evaluated the different semiclassical approximations as reported in Tab.I. We denote the time-points at which the energy gap is calculated by t i and the corresponding energy gap by X i where i = 0 . . .N − 1 runs over the N the time-points.As in Ref. [26] we evaluated the correlator by using a discrete representation, which implements the k-th element of the twotime correlator as where X is the mean.Here, one assumes that the N − k values X i give a faithful initial distribution which reproduces the Boltzmann distribution.To minimize spurious effects in the Fourier transform, we multiplied the time trace by a Gaussian of variance σ 2 gaussian = 0.09 • t 2 max = 2.304 • 10 5 fs 2 with t max = 1600 fs, the length of the correlation function (as reported in [26]).The Gaussian is normalized to have unitary area in frequency domain [57] following our definition of the Fourier transform in Eq. 6, so that in frequency domain this corresponds to a convolution with a Gaussian with a FWHM of 26 cm −1 .Next, we computed the different semiclassical quantities of Table I using our initial time trace.
In Figure 2, we show the temperature-dependent coupling densities TDCDs (as defined in Eq. 6), for site 1 of the FMO complex (site 1 at 77K and 300K) evaluated using the different approximations listed in Table I column two.We notice how, as expected, there are little differences between the approximations at low frequencies.Only at higher frequencies the TDCD differs significantly for each approach.The Egelstaff approximation incorrectly predicts a negative spectral density for low frequencies in this case and was therefore not shown in the plots.
From the general definition of each semiclassical correction, it isn't clear which one is most accurate.To better reason on which one to choose, we will look at the temperature dependence of the spectral density.Further, we will compare to experimental results and finally we will evaluate the three point correlator (Sec.VI D).

B. Analysis of prefactors in terms of temperature dependence of the spectral density
From our discussion in Section I, we recall that many open quantum system approaches rely on the assumption of lin- ear coupling to a bath of harmonic oscillators.This leads to a temperature-independent spectral density j(ω), as discussed in Section III B. Inspection of Fig. 2 shows that for all but the Harmonic approximation the TDCD (from which one obtains J(ω) which is directly proportional to the spectral density j(ω)) obtained from the QM/MM is not similar at different temperatures.This is more apparent at higher frequencies.To gain further insight into this temperature dependence, in Fig. 3, we compare the asymmetric TDCD (J(ω) = πj(ω) ; ω > 0) obtained using the Standard (panel a) and the Harmonic (panel b) approximations for site 1 of the FMO complex.Results for all sites at both temperatures are reported in the Supplementary Information, Appendix A. and the corresponding data files for the Harmonic J(ω) can be downloaded from Ref. [63].One clearly sees that for the Standard correction there is a huge difference between the 77K and the 300K results.However, in the case of the Harmonic correction the spectral densities obtained at the two temperatures nicely lie on top of each other, as one would require for a temperature-independent spectral density.This result suggests that the Harmonic correction is the appropriate one to employ to obtain spectral densities to be used in open quantum system models which assume linear coupling to a bosonic bath.Note, that the good agreement at both temperatures for the Harmonic correction might be purely accidental or due to the fact that the MD is not fully converged.We would need to run much longer QM/MM trajectories to improve the statistics and check the convergence of the distributions.This lack of statistics could also explain the fact that for the SD aver-aged over all chomophores (panels c) and d)), the agreement between both temperatures is slightly better than for the individual sites.
Finally we would like to remark that a temperature dependence of the reorganization energy has been observed in the context of electron transfer (ET) donor-acceptor energy gap spectral densities [58,59].

C. Comparison to experimental spectral density
In Fig. 4 panel a) and b) we compare the asymmetric TDCD for site 3 (Standard and Harmonic correction), with the asymmetric TDCD obtained from fluorescence line narrowing (FLN) experiments [28].We focus on the low frequency part (up to ~500 cm −1 ), which is relevant for energy transfer in the FMO complex.The FLN results are obtained from the lowest excitonic peak of the FMO absorption spectrum which is believed to be generated almost entirely by BChl 3. Therefore, we compare the experiment to the theoretical spectral density obtained from the QM/MM for BChl 3.
The experimental spectral density shown in Fig. 4 is based on the dotted curve jexp (ω) of Fig. 2 of Ref. [29], which is in good agreement with the one-phonon vibrational profile (OPVP) of Ref. [28], because of the small total Huang-Rhys factor.Note that the extraction of the OPVP uses the same model of shifted harmonic potential surfaces as we did in Section III B. Thus it corresponds to a SD which is suitable as input in the open system approaches.In this harmonic model the profile jexp (ω) is related to our definition of the spectral density by j exp (ω) = ( ω) 2 jexp (ω).The positive frequency part of the asymmetric TDCD, J(ω) is obtained from the spectral density, as defined in Eq. 13, by J exp (ω) = π j exp (ω).
From panel e) and f) of Fig. 4, we see that the magnitude and overall lineshape of both the Standard and the Harmonic correction are in good agreement with the FLN data, in contrast with previous results [26] [60].
A closer inspection of the curves in panels c) and d) of Fig. 4 shows that the width of the peaks obtained from the QM/MM simulation is much broader than that obtained from the FLN data.As described in Section VI A, this broadening is due to the finite length of the numerical correlator, and to the convolution with a gaussian function in frequency domain, which results in a broadening of FWHM 26 cm −1 .Also, the position of the peaks do not perfectly coincide.There might be various reasons for this discrepancy: The trajectories might be too short, the quantum chemical calculations of the transition gap are not accurate enough, or the thermostat leads to some spurious effects.One has also to keep in mind that there are uncertainties in the experimental data as well.The experimental data (in particular at higher frequencies) probably do not represent the actual spectral density of BChl 3 (excitonic effects might play a relevant role, and it was difficult to extract the lineshape from the representation of Refs.[29] and [28]).
Nevertheless, this good agreement in magnitude and overall lineshape makes us confident, that the QM/MM procedure can indeed be useful to extract spectral densities.
Finally, it seems that the Harmonic correction describes the comparison of J(ω) obtained for site 1 with the Harmonic approximation (Tab.I, second line, third column) at 77K and at 300K.We see clearly that the Harmonic prefactor gives a roughly temperature independent J(ω), while large differences are seen using the Standard prefactor.
FLN data slightly better in terms of amplitude, respect to the Standard correction.

D. Higher-order correlation function
From the theory of discrete processes, similarly to Eq. 49, we see that the (k, j)-th element of the three-time correlator is ) with ∆X i = X i − X where X is the mean and N is the number of time points (as defined in Sec.VI A).We compare the two-time and the three-time correlators by dividing them by increasing powers of the standard deviation s ≡ √ m (2) , thus we use Eq.49 for the two-time correlation function and divide it by s 2 and we divide Eq.50 by s 3 .The results for site 1 of the FMO complex at 77 and 300K are reported in Fig. 5.For the two-time correlator, Fig. 5 panels a) and b), we see correlations up to at least 1000 time steps, while for the threetime correlator, panels c), d), e) and f), we see a rather noisy profile with values about one/two orders of magnitude smaller than the largest value of the two-time correlations.This is observed for all sites and temperatures (Results for all sites can be found in the Supplementary information, Appendix A).
This means that since we find a small three-time correlator, the linear coupling to a harmonic bath assumption is probably good.In fact, as described in Sec.V this case corresponds to linear coupling to the bath and Gaussian correlated bath operators.Of course, the statistics of the three-time correlator is not great due to the finite length of the time trajectories, but we think that the general tendency is correct.One should also keep in mind that there may be fortuitous cases in which the three-time correlator is roughly zero and the bath is not harmonic.Further, this comparison is based on the order of magnitude of the correlations, the three-time correlator is only much smaller.It may be that for some modes of the system, certain frequencies, present in the three-time correlator's two dimensional Fourier transform give a more important contribution to the dynamics than other frequencies present in the spectral density.Nonetheless, the above result encourages the idea that the assumption of linear coupling and harmonic bath is valid.This, in turn, implies that one should use the Harmonic semiclassical correction in Sec.II C, which is also Figure 4. Panel a) shows J(ω) for site 3 of the FMO complex calculated with the Standard approximation at 77 and 300K and the green curve corresponds to the experimental spectral density rescaled by π to obtain J(ω) as defined in Eq. 13 [28,29](More details on the experimental spectral density are given in the text).Panel b) shows J(ω) for site 3 calculated with the Harmonic approximation at 77 and 300K and again the green curve corresponds to the experimental spectral density [28,29].The agreement with the experimental (green) spectral density is slightly better for the Harmonic approximation than for the Standard approximation.Panels c) and d) correspond to the same quantities as those of panels a) and b) in the low frequency region, here we note that both approximations are roughly equivalent for ω k B T < 1 (e.g at T = 77 K for ω < 55 cm −1 and at T = 300 K for ω < 200 cm −1 ).Further, the spectral density, as defined in Eq. 13 can be obtained by dividing J(ω) by π.
consistent with the prefactor found in III.
On a final note, to confirm with certainty that the bath is Harmonic, one should evaluate higher order correlators, beyond the three-time correlator.However, to obtain a statistically relevant estimate, much longer time dependent energy gap trajectories, which are expensive in terms of the QM/MM propagation, would be required.Work in this direction is being carried out in our groups.

VII. CONCLUSIONS
In this work, we have investigated the connection between the gap correlation function extracted from ground state QM/MM and the bath spectral density used as input in many open quantum system approaches.
One important point is that the classical bath correlation function is real while the quantum mechanical one is generally complex.There exist several semiclassical a posteriori corrections which aim to fix this and we have employed them on our time traces to recover a part of the imaginary component.
The discussed prefactors originate from general expansions in orders of and do not include information on the specific type of system-bath coupling, etc.We have investigated two simple models and found that the prefactors obtained correspond to two of the general semiclassical expressions.Thus, we have linked the semiclassical limits with a microscopic potential energy surface picture.
We have shown that the gap-correlation function extracted from ground state QM/MM only corresponds to the fully quantum excited state calculations in the case of shifted parabolas.This model for a few vibrational modes of the chromophores has been successfully used to describe the optical properties of molecular aggregates.Including only a finite number of internal vibrations is probably a good approximation for molecules in the gas phase or suprafluid Helium nanodroplets [54].However, for molecules in solution or when a protein environment is present it is no longer a good approximation to include only only a few (undamped) modes.
In particular, one has to take into account the interaction of the vibrations with the environment in addition to the direct interaction of the electronic excitation with the environment.For this general situation, it is no longer clear whether the model of shifted harmonic potential surfaces is indeed a good description of the system.Therefore, we have investigated whether the approximation of harmonic bath and linear coupling is accurate for our QM/MM calculations for the FMO photosynthetic complex by computing the next higher order correlator beyond the twotime correlator.The three-point correlator seems to give a small contribution which, while not being conclusive, suggests to us that the Harmonic/linear coupling model is a good approximation.The evaluation of the four-time correlation function would be useful to bolster this claim.
The analysis of the temperature dependence of prefactors for the spectral density also suggests that the Harmonic approximation is preferred to use for the FMO complex, and perhaps other photosynthetic complexes, rather than the Standard one when employing it in Open Quantum system approaches.
Having made these choices, the theoretical results are in reasonably good agreement with the experimental spectral density.These result in a much better agreement than in our previous work, which underestimated the magnitude of the spectral density [26] and than other QM/MM calculations [25] which overestimate the coupling to the bath by one order of magnitude.
Finally, we have explained the link between bath correlation function and gap correlation function and found models under which the gap correlation function can actually be viewed as a general open quantum system bath correlation function.] site 7 -harmonic -77K site 7 -harmonic -300K Figure 14.Comparison of the asymmetric component of temperature-dependent coupling density Gasym(ω) = J(ω) ; for sites 5-7 of the Fenna-Matthews-Olson complex, obtained with the Harmonic approximation (As described in the text) at 77K and at 300K.Note that, as described in the text, the spectral density can be obtained by dividing J(ω) by π.

Figure 1 .
Figure 1.Shifted identical harmonic Born-Oppenheimer surfaces, Ω is the frequency of each harmonic potential and δQ is the coordinate shift between the minima of the ground and excited state potentials.This model is the one employed in Sec.III B to derive classical and semiclassical expressions for the Fourier transform of the bath correlation function G(ω) and for the spectral density.

Figure 2 .
Figure 2. Positive frequency part of the temperature-dependent coupling densities G(ω) obtained, as described in the Section VI A with each of the Standard, Harmonic, Schofield and Schofield-Harmonic corrections as (see Tab. I, column two).In panel a) results are at 77K and in panel b) 300K.

Figure 3 .
Figure 3. Panel a) comparison of the asymmetric component of temperature-dependent coupling density J(ω) ≡ Gasym(ω) ; for site 1 of the Fenna-Matthews-Olson complex, obtained with the Standard approximation (Tab.I, first line, third column) at 77K and at 300K.Panel b)comparison of J(ω) obtained for site 1 with the Harmonic approximation (Tab.I, second line, third column) at 77K and at 300K.We see clearly that the Harmonic prefactor gives a roughly temperature independent J(ω), while large differences are seen using the Standard prefactor.

Figure 5 .
Figure 5. Panel a): Two-time correlation function of the energy gap fluctuations of site 1 of the FMO complex, normalized by s 2 (the variance) at 77K after evaluating it as in Eq. 49.Panel b): Two-time correlation function for site 1 at 300K.Panel c) Three-time correlation function of the energy gap fluctuations of site 1 of the FMO complex, as defined in Eq.50 normalized by s 3 at 77K.Panel d) Three-time correlation function for site 1 at 300K.

Figure 6 .
Figure 6.Panel a): Two-time correlation function of the energy gap fluctuations of site 1 of the FMO complex, normalized by the variance s 2 at 77K as discussed in the text.Panel b): Two-time correlation function for site 1 at 300K.Panel c) Three-time correlation function of the energy gap fluctuations of site 1 of the FMO complex, as discussed in the text and multiplied by s 3 at 77K.Panel d) Three-time correlation function for site 1 at 300K.

Figure 7 .
Figure 7. Same as for Fig. 6 but for site 2 of the FMO complex.

Figure 8 .
Figure 8. Same as for Fig. 6 but for site 3 of the FMO complex.

Figure 9 .
Figure 9. Same as for Fig. 6 but for site 4 of the FMO complex.

Figure 10 .
Figure 10.Same as for Fig. 6 but for site 5 of the FMO complex.

Figure 11 .
Figure 11.Same as for Fig.6but for site 6 of the FMO complex.

Figure 12 .Figure 13 .
Figure 12.Same as for Fig. 6 but for site 7 of the FMO complex.