Ab initio computation of two-dimensional electronic spectra is an expanding field, whose goal is improving upon simple, few-dimensional models often employed to explain experiments. Here, we propose an accurate and computationally affordable approach, based on the single-trajectory semiclassical thawed Gaussian approximation, to evaluate two-dimensional electronic spectra. Importantly, the method is exact for arbitrary harmonic potentials with mode displacement, changes in the mode frequencies, and inter-mode coupling (Duschinsky effect), but can also account partially for the anharmonicity of the involved potential energy surfaces. We test its accuracy on a set of model Morse potentials and use it to study anharmonicity and Duschinsky effects on the linear and two-dimensional electronic spectra of phenol. We find that in this molecule, the anharmonicity effects are weak, whereas the Duschinsky rotation and the changes in the mode frequencies must be included in accurate simulations. In contrast, the widely used displaced harmonic oscillator model captures only the basic physics of the problem but fails to reproduce the correct vibronic lineshape.

Electronic spectroscopy allows us to study excited electronic states and light-induced nuclear dynamics. To track this ultrafast dynamics on femtosecond time scales, a range of time-resolved and two-dimensional spectroscopic techniques were developed. The complex signals obtained in these experiments are, however, difficult to interpret without the help of theoretical modeling.1 

Most models for describing two-dimensional electronic spectra treat the electronic states as the system and include nuclear dynamics only approximately as bath effects. The simplest approach assumes that nuclear degrees of freedom induce a Gaussian-like or Lorentzian-like broadening, neglecting completely the coherent nuclear dynamics. Such models are appropriate only if the coherent dynamics is strongly suppressed by the surrounding solvent dynamics. To account for inhomogeneous (static) broadening, the energy gap between the electronic states can be averaged over snapshots of different arrangements of the environment.2–4 Alternatively, a swarm of trajectories could be used in Kubo-type calculations, where each trajectory is equipped with a phase, obtained from the time integral of the energy gap between the involved electronic states along the trajectory, and the correlation functions are averaged over the full ensemble;5–10 these approaches are also known as phase averaging,11 Wigner-averaged classical limit,12–14 or dephasing representation.15–17 Such methods are accurate when the curvatures of different potential energy surfaces are similar and in the limit of strong dephasing, i.e., for short times. Even in this case, however, if the on-the-fly dynamics is performed with ab initio electronic structure methods, evolving the full ensemble of classical trajectories can quickly become prohibitively expensive. To account for both intramolecular and intermolecular nuclear dynamics in two-dimensional spectroscopy, one can use the multimode Brownian oscillator model,11,18–20 which considers a few primary harmonic modes coupled to a large number of low-frequency bath oscillators. For this model, the spectra can be computed analytically; moreover, the parameters of the model can be computed efficiently from a single ab initio classical trajectory, as demonstrated in Refs. 21 and 22. This allows one to perform electronic structure computations at a high level of theory using, for example, post-Hartree–Fock multiconfigurational wavefunction methods. Yet, the approach is limited to modeling the molecule as a set of few uncoupled displaced harmonic oscillators (DHOs). Such a simple description is inadequate in systems that exhibit strong mode–mode coupling, changes in the force constants between the ground and excited electronic states, or anharmonicity effects. A generalization to a set of uncoupled harmonic oscillators with both displacement and a possible change in the force constant was proposed by Fidler and Engel, who used the approximate third-order cumulant expansion.23 Very recently, the third-order cumulant expansion was also studied as an approximate way to account for the mode–mode coupling (Duschinsky effect) in both linear24 and nonlinear25 spectra.

There has been little development in the ab initio simulation of two-dimensional electronic spectra beyond the standard semiclassical methods or displaced harmonic models. To account for anharmonicity effects26,27 or more general coupled oscillators, one is forced to employ computationally expensive exact quantum dynamics methods,28–33 such as different flavors of the multiconfigurational time-dependent Hartree (MCTDH) method,34,35 or the hierarchical equations of motion.36,37 These methods require the pre-computation of the full potential energy surfaces and are not suitable for a first-principles on-the-fly implementation. First-principles multi-trajectory semiclassical approaches38–46 and direct quantum dynamics methods, which often use multiple Gaussians,45,47–53 also called coherent or Davydov states,54–56 to represent the evolving wavepacket, are impractical due to the large number of required ab initio evaluations.

Here, we propose an efficient semiclassical method to evaluate vibrationally resolved two-dimensional electronic spectra. The approach, based on Heller’s single-trajectory thawed Gaussian approximation,57 accounts for inter-mode coupling, changes in the force constants, and, at least partially, for the anharmonicities of the ground- and excited-state potential energy surfaces. First, we study how the accuracy of the method depends on the degree of anharmonicity in the one-dimensional Morse system. The results are compared with the exact benchmark and with the harmonic approximation, which neglects the anharmonicity completely. Second, we analyze the effects of Duschinsky coupling and anharmonicity on the linear absorption and two-dimensional spectra of phenol.

The central object in all types of third-order electronic spectroscopy is the third-order polarization11,58

P(3)(t)=0dt30dt20dt1R(t3,t2,t1)E(tt3)×E(tt3t2)E(tt3t2t1),
(1)

where E(t) is the electric field of light (without the polarization vector) and

R(t3,t2,t1)=i3α=14[Rα(t3,t2,t1)Rα(t3,t2,t1)*]
(2)

is the third-order response function, expressed in terms of correlation functions

R1(t3,t2,t1)=C(t2,t3,t1+t2+t3),
(3)
R2(t3,t2,t1)=C(t1+t2,t3,t2+t3),
(4)
R3(t3,t2,t1)=C(t1,t2+t3,t3),
(5)
R4(t3,t2,t1)=C(t3,t2,t1),
(6)

and

C(τa,τb,τc)=Trρ^μ^eiĤ2τa/μ^eiĤ1τb/μ^eiĤ2τc/μ^eiĤ1(τa+τbτc)/.
(7)

In Eq. (7), Ĥi are the vibrational Hamiltonians of the ground (“1”) and excited (“2”) electronic states, ρ^=exp(βĤ1)/Tr[exp(βĤ1)] is the vibrational density operator in the ground electronic state, and μ^ = μ^21=μ^21ϵ is the electronic transition dipole moment projected on the polarization unit vector ϵ of the external electric field. Equations (3)–(6) rely on the following assumptions: (i) due to the large gap between the electronic states, only the ground electronic state is initially populated; (ii) Born–Oppenheimer approximation, i.e., there is no population transfer under field-free evolution; (iii) light pulses are linearly polarized in the same direction ϵ; and (iv) only two electronic states are involved. In the following, we discuss how to evaluate the components Rα of the response function (2).

In the zero-temperature limit, we assume that only the ground (“g”) vibrational state |1, g⟩ of the ground electronic state is populated initially, i.e., ρ^=|1,g1,g|. Then, we may rewrite Eq. (7) in terms of nuclear wavepackets as

C(τa,τb,τc)=1,g|μ^eiĤ2τa/μ^eiĤ1τb/μ^eiĤ2τc/μ^eiĤ1(τa+τbτc)/|1,g
(8)
=1,g|μ^eiĤ2τa/μ^eiĤ1τb/μ^eiĤ2τc/μ^|1,geiω1,g(τa+τbτc)
(9)
=1,g|μ^eiĤ2τa/μ^eiĤ1τb/μ^eiĤ2τc/μ^|1,g
(10)
=ϕτb,τa|ϕ0,τc,
(11)

where ω1,g=1,g|Ĥ1|1,g, Ĥi=Ĥiω1,g, and

ϕτ,t=eiĤ1τ/μ^eiĤ2t/μ^|1,g.

The result (11) has an appealing interpretation in terms of bra and ket wavepackets, which we represent pictorially for R3 [Eq. (5)] in Fig. 1. The bra wavepacket is first evolved in the excited electronic state for a time τa = t1 and then for a time τb = t2 + t3 in the ground state, where it is a non-stationary wavepacket due to the initial t1 dynamics on the excited-state potential energy surface. The ket wavepacket “waits” during the t1 and t2 times and is only evolved for a time τc = t3 in the excited-state. This simple picture has been discussed in the literature in the context of pump–probe59 and two-dimensional28 spectroscopy. In general, during the t1 (coherence) and t3 (detection) times, the bra and ket wavepackets evolve on different potential energy surfaces, i.e., the system is in a state of electronic coherence; during t2, also called population or waiting time, both nuclear wavepackets are in the same electronic state, i.e., the system is in an electronic population state.19 

FIG. 1.

Evolution of the bra [(a), dotted line] and ket [(b), solid line] wavepackets of Eqs. (7) and (11) for τa = t1, τb = t2 + t3, and τc = t3. Their overlap (c) is the R3(t3, t2, t1) component of the third-order response function R(t3, t2, t1) [see Eqs. (5), (7), and (11)].

FIG. 1.

Evolution of the bra [(a), dotted line] and ket [(b), solid line] wavepackets of Eqs. (7) and (11) for τa = t1, τb = t2 + t3, and τc = t3. Their overlap (c) is the R3(t3, t2, t1) component of the third-order response function R(t3, t2, t1) [see Eqs. (5), (7), and (11)].

Close modal

The evaluation of Rα functions requires only one excited-state wavepacket evolution up to time t1 + t2 + t3 and, in addition, wavepackets propagated in the ground electronic state starting from the snapshots along the excited-state trajectory. Since such calculations would be difficult to perform with multiple-trajectory direct dynamics methods, we employ the efficient, single-trajectory thawed Gaussian approximation.

Within the thawed Gaussian approximation, the time-dependent wavepacket takes the form of a Gaussian

ψt(q)=ei[12(qqt)TAt(qqt)+ptT(qqt)+γt],
(12)

where qt and pt are D-dimensional position and momentum vectors, At is a complex and symmetric D × D matrix with a positive-definite imaginary part, and γt is a complex number whose imaginary part ensures the normalization. D is the number of degrees of freedom. The wavepacket (12) solves exactly the time-dependent Schrödinger equation

i|ψ̇t=[T(p^)+VLHAq^]|ψt,
(13)

where T(p)=12pTm1p, m is a D × D symmetric mass matrix, and

VLHA(q)=V(qt)+V(qt)T(qqt)+12(qqt)TV(qt)(qqt)
(14)

is the local harmonic approximation of the true potential V(q) about qt, if the Gaussian’s parameters satisfy the system57 

q̇t=m1pt,
(15)
ṗt=V(qt),
(16)
Ȧt=Atm1AtV(qt),
(17)
γ̇t=Lt+i2Tr(m1At).
(18)

In Eq. (18), Lt = T(pt) − V(qt) is the Lagrangian of the classical trajectory (qt, pt). The above equations are interpreted as follows: the phase-space center (qt, pt) of the Gaussian (12) evolves classically with the exact classical Hamiltonian, the complex matrix At evolves according to the Hessian computed at the current qt, and the complex number γt is updated according to the Lagrangian of the classical trajectory (qt, pt) and the matrix At. Since the only source of error is the local harmonic approximation (14), the thawed Gaussian propagation is exact for arbitrary, multi-dimensional harmonic potentials.

The method was originally proposed for problems involving short-time dynamics, such as photodissociation spectra.60 However, its accuracy appears to be surprisingly satisfactory in molecular systems even at longer times because many molecules are only weakly to moderately anharmonic.61,62 Using a single thawed Gaussian wavepacket, which is the essence of Heller’s thawed Gaussian approximation, is rather restrictive but also very efficient for on-the-fly dynamics coupled to the ab initio electronic structure. The on-the-fly ab initio thawed Gaussian approximation63 proved useful in treating efficiently anharmonicity effects on linear absorption,61,62,64,65 emission,65,66 and photoelectron spectra,64 as well as in understanding nuclei-induced electronic decoherence in attosecond experiments.67 Recently, we extended our on-the-fly implementation of the single-Gaussian approach to simulate frequency- and time-resolved pump–probe spectra,68 similar to the earlier work by Rohrdanz and Cina69 and Braun et al.70 on model potentials. Although ensembles of thawed Gaussians were largely discarded and replaced by frozen Gaussians due to the numerical instabilities that often appear in nonadiabatic dynamics simulations, thawed Gaussians are being reintroduced, e.g., in the semiclassical hybrid dynamics71–73 or Gaussian-based MCTDH,74–76 and especially for spectroscopic applications,77–80 due to their ability to describe couplings between different degrees of freedom.81 

A variety of different third-order experiments can be simulated through the computation of the response function (2).11 For example, the full response function is needed for evaluating transient absorption spectra with finite-duration pulses.82 Here, we focus on the two-dimensional electronic spectra,

Sα(ω3,ω1)=Re0dt30dt1Rα(t3,0,t1)eiω3t3±iω1t1,
(19)

obtained from the individual correlation functions Rα with t2 = 0 (i.e., at zero delay time). Spectra at nonzero delay could be obtained by using t2 > 0 in Eq. (19). Spectra Sα represent the ideal signals obtained in the limit of ultrashort pulses. In a more general setting with finite pulses, the two-dimensional spectra are computed from the time-dependent polarization (1), which involves explicitly the electric fields.19,83 To ensure that all spectra appear at positive frequencies ω1, nonrephasing spectra (α = 1, 4) are computed with the positive sign in the exponent of Eq. (19), while the negative sign is used for the rephasing spectra (α = 2, 3).19 Furthermore, it is easy to see from Eqs. (3)–(6) that R1(t3, 0, t1) = R4(t3, 0, t1) and R2(t3, 0, t1) = R3(t3, 0, t1); hence, we will show only two sets of spectra: S1S4 and S2S3. In general (i.e., for arbitrary t2), correlation functions R1 and R2 are associated with the stimulated emission process because the system evolves in the excited state during the population time t2; functions R3 and R4 correspond to the ground-state bleaching because the system is in the ground electronic state during the t2 delay time. For t2 = 0, one cannot distinguish between these two processes.

To analyze the accuracy of different approximate approaches, we introduce the spectral contrast angle

cosθ=S(ref)SS(ref)S
(20)

between the reference [S(ref)] and approximate (S) spectra, where

S(1)S(2)=dω1dω3S(1)(ω3,ω1)S(2)(ω3,ω1)
(21)

is the inner product of two two-dimensional spectra and S=SS is the associated norm.

An arbitrary one-dimensional harmonic potential,

VHarmonic(q;Veq,qeq,ω)=Veq+12mω2(qqeq)2,
(22)

is described by the equilibrium position qeq, energy minimum Veq, and frequency ω. We set the mass m = 1 in all of our model calculations. Let us also define a one-dimensional Morse potential,

VMorse(q;Veq,qeq,ω,χ)=Veq+ω4χ1e2mωχ(qqeq)2,
(23)

in terms of the anharmonicity parameter χ and the parameters Veq, qeq, and ω, which relate to the harmonic potential (22) fit to the Morse potential at qeq.

We construct a set of one-dimensional systems composed of the ground-state harmonic potential,

V1(q)=VHarmonic(q;V1,eq=0,q1=0,ω1=1),
(24)

and the excited-state Morse potentials,

V2(q)=VMorse(q;V2,eq=0,q2=1.5,ω2=0.9,χ),
(25)

of variable anharmonicity χ ranging from 0.006 to 0.02. The initial vibrational state, i.e., the ground vibrational state of the ground electronic state, is a Gaussian due to the ground-state harmonic potential. The exact two-dimensional electronic spectra are compared to approximate spectra evaluated either with the harmonic approximation or with the thawed Gaussian approximation. Within the harmonic approximation, the excited-state Morse potential is replaced by the harmonic potential

V2(q)VHarmonic(q;V2,eq,q2,ω2).
(26)

Note that the harmonic result does not depend on the anharmonicity parameter χ of the Morse potential.

Next, we compare the harmonic and thawed Gaussian approximations for a one-dimensional system composed of two Morse potentials

V1(q)=VMorse(q;V1,eq,q1,ω1,χ=0.01),
(27)
V2(q)=VMorse(q;V2,eq,q2,ω2,χ=0.01),
(28)

with the same degree of anharmonicity. The exact initial state is no more a Gaussian. However, in the thawed Gaussian simulations, we approximate it by the vibrational ground state of the harmonic potential (24) fitted to the ground-state Morse potential (27) at its minimum. The harmonic approximation replaces both ground-state and excited-state potential energy surfaces by the harmonic potentials; the result is the same as for the harmonic-Morse system described above.

Wavepacket propagation was performed for 150 steps in both t1 and t3 times and with a time step of 0.2. The transition dipole moment was set to 1 (Condon approximation). Correlation functions R1, R2, R3, and R4 were multiplied by a Gaussian damping function exp[a(t12+t32)] with a = 0.014 427, resulting in the Gaussian broadening (half-width at half-maximum of 0.2) of the spectra along both frequency axes. The exact spectra were computed in the eigenstate representation, which is feasible for these one-dimensional systems since both harmonic and Morse eigenfunctions are known;84 the associated Franck–Condon overlaps were computed numerically.

The electronic structure of phenol was modeled using the density functional theory with the PBE0 functional and 6-311G(d, p) basis set, as implemented in the Gaussian 16 quantum chemistry package.85 Excited-state calculations were performed with the time-dependent density functional theory. This choice of electronic structure theory provides ground-state frequencies similar to those computed at the MP2/aug-cc-pVDZ level (see Table II of the supplementary material and Ref. 86) and transition energies along the excited-state trajectory that agree, up to an approximately constant shift (which results only in a shift of the computed spectrum but does not affect its shape), to those evaluated at the EOM-CCSD/6-311G(d, p) level (Fig. 1 of the supplementary material). A single ab initio excited-state trajectory was run for 1000 steps starting from the ground-state optimized geometry; subsequent ground-state classical trajectories were propagated for 500 steps. Overall, the calculations allow the evaluation of the correlation functions with 500 steps in both t1 and t3 delay times; t2 delay was set to zero. All dynamics simulations used a time step of 0.25 fs and a standard second-order Verlet integrator. The ab initio calculations evaluated not only the energies and gradients at each step but also the Hessians of the electronic energy. These potential energy data were transformed to ground-state normal mode coordinates and used to propagate the 33-dimensional wavepacket according to Eqs. (15)–(18). After evolving the ground- and excited-state Gaussian wavepackets, the correlation functions were computed using Eqs. (3)–(6) and (11) and the expression

ψ1|ψ2=(2π)Ddet(iδA)expi12δξT(δA)1δξ+δη
(29)

for the overlap of two thawed Gaussian wavepackets with parameters qi, pi, Ai, and γi (i = 1, 2). In Eq. (29), we defined vectors and scalars

ξipiAiqi,
(30)
ηiγi12(ξi+pi)Tqi,
(31)

as well as the notation δΛΛ2Λ1* for Λ = A, ξ, and η.

To construct the harmonic model, also known as the generalized Brownian oscillator model,24 of phenol, an additional Hessian was computed at the optimized excited-state geometry. This corresponds to the so-called adiabatic Hessian or adiabatic harmonic model.63,87 Two more approximate models were also studied: The uncoupled harmonic model was obtained by neglecting the off-diagonal terms of the excited-state Hessian expressed in terms of the ground-state normal modes. The displaced harmonic oscillator model, also called the Brownian oscillator model, was constructed by replacing the excited-state Hessian in the adiabatic harmonic model by the ground-state Hessian; this specific way of constructing the displaced harmonic oscillator parameters is called the adiabatic shift approach.63,87 When applied to any of these different harmonic potentials, the thawed Gaussian propagation is exact and enables an efficient evaluation of linear and two-dimensional spectra. Although explicit expressions are available for the evaluation of linear absorption and emission spectra of harmonic systems,88 no such analytical approaches have been presented for the two-dimensional spectra of arbitrarily shifted, distorted, and rotated harmonic potentials.

Spectra simulations assumed Condon approximation for the transition dipole moment. Linear absorption spectra were computed from the first 500 steps of the excited-state wavepacket autocorrelation function (see Fig. 2 of the supplementary material, where the convergence is confirmed, and Ref. 63 for more details) and were broadened by a Gaussian with half-width at half-maximum of 120 cm−1; the same broadening was used for the two-dimensional spectra along both ω1 and ω3 frequency axes. This corresponds to a phenomenological inhomogeneous broadening; homogeneous broadening due to direct system-bath interactions is neglected. The system-bath coupling would be needed for spectra at later delay times t2 > 0 as the system would have time to relax and dissipate energy to the environment; we assume that the response functions with t2 = 0 and t1, t3 < 125 fs are only weakly affected by the system-bath coupling.

1. Harmonic-Morse system

Two-dimensional spectra for the harmonic ground-state potential and Morse excited-state potential are shown in Fig. 2. The exact nonrephasing spectrum appears only along the diagonal, whereas the rephasing spectrum exhibits a characteristic checkerboard pattern due to vibronic transitions that involve various ground- and excited-state vibrational states. Already at first sight, it is clear that the spectra evaluated within the thawed Gaussian approximation reproduce the exact spectra well, which is not the case for the harmonic results. The nonrephasing harmonic spectrum flattens out at higher frequencies (spectral region A indicated in the top right panel of Fig. 2), unlike the exact and thawed Gaussian spectra, which exhibit clear vibronic peaks at these frequencies. Similar effects are seen in the rephasing spectra, mostly in the spectral region labeled B (see Fig. 2, bottom right). Again, the exact spectrum is composed of a long vibronic progression up to ω3 = 6, which, in the harmonic approximation, is truncated around ω3 = 4. In region C, the harmonic spectrum is missing negative vibronic peaks, which are reproduced well by the thawed Gaussian approximation. The thawed Gaussian approximation, however, suffers from another form of error: as in linear spectroscopy (see, e.g., Ref. 64), artificial negative peaks may also appear in the two-dimensional spectra, which is most obvious in the nonrephasing spectrum of Fig. 2 around (ω1, ω3) ≈ (−1, −1).

FIG. 2.

Exact, thawed Gaussian, and harmonic two-dimensional electronic spectra of the harmonic-Morse system described in Sec. III A with the anharmonicity of the excited-state Morse potential χ = 0.01. Spectral regions A, B, and C, discussed in the main text, are indicated on the harmonic spectra.

FIG. 2.

Exact, thawed Gaussian, and harmonic two-dimensional electronic spectra of the harmonic-Morse system described in Sec. III A with the anharmonicity of the excited-state Morse potential χ = 0.01. Spectral regions A, B, and C, discussed in the main text, are indicated on the harmonic spectra.

Close modal

We now compare the exact and approximate spectra at different levels of anharmonicity by measuring the error (see Fig. 3) through the spectral contrast angle [Eq. (20)]. The thawed Gaussian approximation exhibits smaller errors in the computed spectra than the harmonic approximation at all levels of anharmonicity and for both rephasing and nonrephasing spectra. As expected, the accuracy of both approximate approaches deteriorates as the anharmonicity of the system increases.

FIG. 3.

Errors of the thawed Gaussian and harmonic spectra of the harmonic-Morse system, measured by the spectral contrast angles [Eq. (20)] at different values of the anharmonicity parameter χ.

FIG. 3.

Errors of the thawed Gaussian and harmonic spectra of the harmonic-Morse system, measured by the spectral contrast angles [Eq. (20)] at different values of the anharmonicity parameter χ.

Close modal

2. Morse–Morse system

Two-dimensional rephasing spectra of the system composed of two Morse potentials, both with the anharmonicity parameter χ = 0.01, are shown in Fig. 4. As in the harmonic-Morse system, the errors of the harmonic spectrum are observed in the spectral regions B and C; the accuracy of the thawed Gaussian spectrum is not much affected by the additional anharmonicity in the ground-state potential surface. To analyze further the two approximate methods, we inspect one-dimensional cuts of the two-dimensional spectra along two different values of ω1 frequency (Fig. 5). We see clearly that the thawed Gaussian approximation recovers the positions and intensities of the vibronic peaks both at low ω1 ≈ 1 and high ω1 ≈ 4 frequencies. Harmonic results qualitatively recover the spectral cut at the lower ω1 frequency (Fig. 5, top) but fail to recover the vibronic peaks at the higher ω1 frequency (Fig. 5, bottom). Notably, the negative peak at ω3 ≈ −1 is missing in the spectrum calculated within the harmonic approximation. Such errors could, in practice, seriously affect the interpretation of the experiments. One of the main challenges in two-dimensional electronic spectroscopy is to assign spectral features to either vibrational or electronic degrees of freedom.89–91 If the simulation, for example, based on a model harmonic potential, cannot reproduce the vibronic peaks found in the experimental spectra, these peaks might end up incorrectly assigned to another electronic state or another excitation process.

FIG. 4.

Exact, thawed Gaussian, and harmonic rephasing spectra of the Morse–Morse system described in Sec. III A (anharmonicity χ = 0.01).

FIG. 4.

Exact, thawed Gaussian, and harmonic rephasing spectra of the Morse–Morse system described in Sec. III A (anharmonicity χ = 0.01).

Close modal
FIG. 5.

One-dimensional cuts of the two-dimensional spectra of Fig. 4 at ω1 ≈ 1 (top) and ω1 ≈ 4 (bottom).

FIG. 5.

One-dimensional cuts of the two-dimensional spectra of Fig. 4 at ω1 ≈ 1 (top) and ω1 ≈ 4 (bottom).

Close modal

Phenol is an ultraviolet chromophore present in proteins as the residue of the naturally occurring amino acid tyrosine. Recently, accurate electronic structure methods were employed to simulate its two-dimensional electronic spectrum94,95 in an attempt to theoretically explore the capabilities of this spectroscopic technique to resolve the features specific to chromophore–chromophore interactions in oligopeptides and, more generally, in proteins.3,96 These recent calculations included multiple electronic states but neglected the vibronic structure of the individual electronic transitions. Here, we present a complementary result: we focus only on the ground and first excited electronic states, i.e., we neglect the excited-state absorption process, but study in detail the vibronic lineshape of the ground-state bleaching/stimulated emission signal. The methods we use neglect the nonadiabatic effects; this is an acceptable approximation for the dynamics in the first excited state of phenol, as demonstrated by the MCTDH simulations performed on a vibronic-coupling Hamiltonian model of phenol.86 

The linear absorption spectrum of phenol was computed with four different approximate methods: the on-the-fly ab initio thawed Gaussian approximation, harmonic approximation, uncoupled harmonic model, and displaced harmonic oscillator model (Fig. 6). Harmonic and on-the-fly thawed Gaussian spectra (Fig. 6, top) are similar in accuracy for this specific system: while the thawed Gaussian propagation results in more accurate intensities of the low-frequency peaks, namely, the 0–0 transition at ≈36 350 cm−1 and the shoulder at ≈36 800 cm−1, the harmonic approximation gives a better estimate of the high-frequency region and the tail of the spectrum. One of the main disadvantages of the thawed Gaussian approximation, the appearance of artificial negative spectral intensities, shows up clearly in the absorption spectrum of phenol. Although the simulated harmonic and thawed Gaussian spectra resemble the experiment, there are remaining differences, most notably in the intensities of the spectral peaks. These errors could be either due to anharmonicity effects not captured by the approximate thawed Gaussian wavepacket propagation or due to the errors in the potential energy data evaluated with an approximate electronic structure method. A fairly small difference between the harmonic and thawed Gaussian spectra suggests that the anharmonicity effects are weak and that the remaining errors in our simulation are due to the inaccuracies of the electronic structure theory. In the bottom panel of Fig. 6, we show the results of two more approximate approaches—the uncoupled harmonic and displaced harmonic oscillator models. These spectra clearly deviate from the experiment, indicating the importance of both mode distortion (changes in mode frequencies) and intermode couplings (Duschinsky effect). When going from the displaced harmonic model, which neglects mode distortion, to the uncoupled harmonic model, which includes mode distortion, the peaks broaden but still exhibit inaccurate intensities. Additional inclusion of the Duschinsky effect, which is achieved by moving to the (coupled) harmonic model, improves the intensities.

FIG. 6.

Experimental linear absorption spectrum of Ref. 92 (data extracted with WebPlotDigitizer93) and the spectra computed with the on-the-fly ab initio thawed Gaussian approximation, harmonic approximation, uncoupled harmonic model (“Uncoupled”), and displaced harmonic oscillator (DHO) model. For ease of comparison, the computed spectra are shifted in frequency and scaled in intensity so that they all match at the maximum of the experimental spectrum (see Sec. IV of the supplementary material for details).

FIG. 6.

Experimental linear absorption spectrum of Ref. 92 (data extracted with WebPlotDigitizer93) and the spectra computed with the on-the-fly ab initio thawed Gaussian approximation, harmonic approximation, uncoupled harmonic model (“Uncoupled”), and displaced harmonic oscillator (DHO) model. For ease of comparison, the computed spectra are shifted in frequency and scaled in intensity so that they all match at the maximum of the experimental spectrum (see Sec. IV of the supplementary material for details).

Close modal

The two-dimensional spectra simulated with different approximate methods are shown in Fig. 7. Again, the uncoupled harmonic and displaced harmonic oscillator models predict the spectra that differ substantially from the harmonic and thawed Gaussian results, which are, in turn, similar to each other. Therefore, based on both linear and two-dimensional spectra simulations, we may conclude that the anharmonicity effects are truly weak in the ground and first excited states of phenol, at least in the region explored by the nuclear wavepacket for short time after the photoexcitation. More precisely, the anharmonicity effects are much weaker than the effects of the Duschinsky rotation and frequency changes, which are, in contrast, significant, as demonstrated by the simulations based on the uncoupled or displaced harmonic models. The anharmonicity could, however, play a role at longer simulation times, needed, for example, to simulate high-resolution spectra. We note that most simulations supporting experimental results are nowadays performed with the simplified displaced harmonic oscillator model, which captures the basic physics of the problem, but is inadequate in certain cases, such as the presented example of phenol.

FIG. 7.

Rephasing two-dimensional electronic spectra of phenol at t2 = 0 computed with the on-the-fly ab initio thawed Gaussian approximation, harmonic approximation, uncoupled harmonic model, and displaced harmonic oscillator (DHO) model. The spectra correspond to the stimulated emission (S2) and ground-state bleaching (S3) processes (S2 = S3 for t2 = 0). Computed spectra were shifted along both frequency axes and scaled in intensity as in Fig. 6.

FIG. 7.

Rephasing two-dimensional electronic spectra of phenol at t2 = 0 computed with the on-the-fly ab initio thawed Gaussian approximation, harmonic approximation, uncoupled harmonic model, and displaced harmonic oscillator (DHO) model. The spectra correspond to the stimulated emission (S2) and ground-state bleaching (S3) processes (S2 = S3 for t2 = 0). Computed spectra were shifted along both frequency axes and scaled in intensity as in Fig. 6.

Close modal

Interestingly, the spectrum spans a broad range of frequencies in both ω1 and ω3, which is in stark contrast with the simulations of Ref. 94. More specifically, the broad vibronic ground-state bleach/stimulated emission spectrum is expected to overlap strongly with the excited-state absorption signals of phenol and even with the signals of other amino acid residues (compare our results with those for a noninteracting benzene-phenol dimer in Fig. 3 of Ref. 94). Hence, an accurate treatment of vibronic effects is needed to simulate realistic spectra and to help explain these overlapping, unresolved spectral features. Our results also indirectly support the concluding part of Ref. 94, where a two-color ultraviolet-visible experiment is proposed to resolve transitions to charge-transfer states (see Fig. 6 of Ref. 94), which appear only when the two chromophores are close to each other. In the visible region of frequency ω3, there are fewer spectroscopic transitions, and these charge transfer states could be easily distinguished from the states of the individual chromophores even with broad vibronic features included.

We have presented a new method for simulating vibrationally resolved two-dimensional electronic spectra that is exact for any shifted, distorted, and coupled harmonic model and, in addition, can approximately account for anharmonicity effects. The method, based on the thawed Gaussian approximation, is shown to be superior to the harmonic approximation for a series of Morse models of varying anharmonicity. On the example of phenol, we show that inter-mode couplings and changes in the mode frequencies, both of which are frequently neglected in simulations, can be crucial for recovering the correct vibronic shape of the two-dimensional electronic spectra. In this specific case, the anharmonicity is shown to be weak, which could allow further studies on the nonlinear spectra of phenol based on the harmonic approximation. For example, our results could be augmented by constructing harmonic models with more accurate electronic structure methods in order to simulate excited-state absorption signals. For systems that do exhibit anharmonicity effects, we propose the on-the-fly ab initio thawed Gaussian approximation as a computationally affordable approach beyond harmonic approximation.

Finally, let us also give a short outlook on how to include features that are missing in the current method. First, as a wavefunction method, the thawed Gaussian approximation is not suitable for treating systems at non-zero temperature. We have shown recently that this limitation can be overcome efficiently with the so-called thermo-field dynamics theory.97 Currently, we are exploring the application of this idea to the computation of nonlinear spectra. Second, the method is originally constructed for isolated systems. An obvious, “ab initio way” to augment the system with an environment would be to include a number of solvent molecules directly into the system. To account for inhomogeneous broadening, the dynamics would have to be repeated for different conformations of the solute–solvent system. Alternatively, the bath effects could be treated through a number of low-frequency harmonic oscillators coupled to the system; the procedures for computing the parameters of the bath oscillators are well-studied in the literature. The extensions that include temperature and environment effects would enable accurate and efficient first-principles simulation of the time-resolved (t2 > 0) two-dimensional electronic spectra in the condensed phase.

See the supplementary material for ground- and excited-state optimized geometries, normal-mode frequencies and displacements, validation of the electronic structure method, wavepacket autocorrelation function, and frequency shifts applied to the computed spectra of phenol and also for Refs. 98–105.

The authors acknowledge the financial support from the Swiss National Science Foundation through the NCCR MUST (Molecular Ultrafast Science and Technology) Network and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 683069-MOLEQULE).

The data that support the findings of this study are openly available in Zenodo at http://doi.org/10.5281/zenodo.4121622.

1.
I.
Conti
,
G.
Cerullo
,
A.
Nenov
, and
M.
Garavelli
,
J. Am. Chem. Soc.
142
,
16117
(
2020
).
2.
I.
Rivalta
,
A.
Nenov
,
G.
Cerullo
,
S.
Mukamel
, and
M.
Garavelli
,
Int. J. Quantum Chem.
114
,
85
(
2014
).
3.
A.
Giussani
,
J.
Marcheselli
,
S.
Mukamel
,
M.
Garavelli
, and
A.
Nenov
,
Photochem. Photobiol.
93
,
1368
(
2017
).
4.
R.
Borrego-Varillas
,
A.
Nenov
,
L.
Ganzer
,
A.
Oriana
,
C.
Manzoni
,
A.
Tolomelli
,
I.
Rivalta
,
S.
Mukamel
,
M.
Garavelli
, and
G.
Cerullo
,
Chem. Sci.
10
,
9907
(
2019
).
5.
S.
Mukamel
,
J. Chem. Phys.
77
,
173
(
1982
).
6.
N. E.
Shemetulskis
and
R. F.
Loring
,
J. Chem. Phys.
97
,
1217
(
1992
).
7.
Z.
Li
,
J. Y.
Fang
, and
C. C.
Martens
,
J. Chem. Phys.
104
,
6919
(
1996
).
8.
C. P.
Van Der Vegte
,
A. G.
Dijkstra
,
J.
Knoester
, and
T. L. C.
Jansen
,
J. Phys. Chem. A
117
,
5970
(
2013
).
9.
R.
Tempelaar
,
C. P.
Van Der Vegte
,
J.
Knoester
, and
T. L. C.
Jansen
,
J. Chem. Phys.
138
,
164106
(
2013
).
10.
A. S.
Petit
and
J. E.
Subotnik
,
J. Chem. Phys.
141
,
154108
(
2014
).
11.
S.
Mukamel
,
Principles of Nonlinear Optical Spectroscopy
, 1st ed. (
Oxford University Press
,
New York
,
1999
).
12.
S. A.
Egorov
,
E.
Rabani
, and
B. J.
Berne
,
J. Chem. Phys.
108
,
1407
(
1998
).
13.
S. A.
Egorov
,
E.
Rabani
, and
B. J.
Berne
,
J. Chem. Phys.
110
,
5238
(
1999
).
14.
Q.
Shi
and
E.
Geva
,
J. Chem. Phys.
121
,
3393
(
2004
).
15.
J.
Vaníček
,
Phys. Rev. E
70
,
055201
(
2004
).
16.
J.
Vaníček
,
Phys. Rev. E
73
,
046204
(
2006
).
17.
T.
Zimmermann
and
J.
Vaníček
,
J. Chem. Phys.
141
,
134102
(
2014
).
18.
A.
Nemeth
,
F.
Milota
,
T.
Mančal
,
V.
Lukeš
,
H. F.
Kauffmann
, and
J.
Sperling
,
Chem. Phys. Lett.
459
,
94
(
2008
).
19.
G. S.
Schlau-Cohen
,
A.
Ishizaki
, and
G. R.
Fleming
,
Chem. Phys.
386
,
1
(
2011
).
20.
J. R.
Caram
,
A. F.
Fidler
, and
G. S.
Engel
,
J. Chem. Phys.
137
,
024507
(
2012
).
21.
A.
Nenov
,
A.
Giussani
,
B. P.
Fingerhut
,
I.
Rivalta
,
E.
Dumont
,
S.
Mukamel
, and
M.
Garavelli
,
Phys. Chem. Chem. Phys.
17
,
30925
(
2015
).
22.
A.
Picchiotti
,
A.
Nenov
,
A.
Giussani
,
V. I.
Prokhorenko
,
R. J. D.
Miller
,
S.
Mukamel
, and
M.
Garavelli
,
J. Phys. Chem. Lett.
10
,
3481
(
2019
).
23.
A. F.
Fidler
and
G. S.
Engel
,
J. Phys. Chem. A
117
,
9444
(
2013
).
24.
T. J.
Zuehlsdorff
,
A.
Montoya-Castillo
,
J. A.
Napoli
,
T. E.
Markland
, and
C. M.
Isborn
,
J. Chem. Phys.
151
,
074111
(
2019
).
25.
T. J.
Zuehlsdorff
,
H.
Hong
,
L.
Shi
, and
C. M.
Isborn
,
J. Chem. Phys.
153
,
044127
(
2020
).
26.
A.
Galestian Pour
,
C. N.
Lincoln
,
V.
Perlík
,
F.
Šanda
, and
J.
Hauer
,
Phys. Chem. Chem. Phys.
19
,
24752
(
2017
).
27.
A.
Anda
,
D.
Abramavičius
, and
T.
Hansen
,
Phys. Chem. Chem. Phys.
20
,
1642
(
2018
).
28.
A.
Schubert
and
V.
Engel
,
J. Chem. Phys.
134
,
104304
(
2011
).
29.
J.
Krčmář
,
M. F.
Gelin
, and
W.
Domcke
,
J. Chem. Phys.
143
,
074308
(
2015
).
30.
M.
Sala
and
D.
Egorova
,
Chem. Phys.
481
,
206
(
2016
).
31.
S.
Choi
and
J.
Vaníček
,
J. Chem. Phys.
150
,
204112
(
2019
).
32.
J.
Roulet
,
S.
Choi
, and
J.
Vaníček
,
J. Chem. Phys.
150
,
204113
(
2019
).
33.
S.
Choi
and
J.
Vaníček
,
J. Chem. Phys.
151
,
234102
(
2019
).
34.
Multidimensional Quantum Dynamics: MCTDH Theory and Applications
, 1st ed., edited by
H.-D.
Meyer
,
F.
Gatti
, and
G. A.
Worth
(
Wiley-VCH
,
Weinheim
,
2009
).
35.
J.
Krčmář
,
M. F.
Gelin
, and
W.
Domcke
,
Chem. Phys.
422
,
53
(
2013
).
36.
J.
Xu
,
R. X.
Xu
,
D.
Abramavicius
,
H. D.
Zhang
, and
Y. J.
Yan
,
Chin. J. Chem. Phys.
24
,
497
(
2011
).
37.
Y.
Tanimura
,
J. Chem. Phys.
153
,
020901
(
2020
).
38.
J.
Tatchen
and
E.
Pollak
,
J. Chem. Phys.
130
,
041103
(
2009
).
39.
M.
Ceotto
,
S.
Atahan
,
S.
Shim
,
G. F.
Tantardini
, and
A.
Aspuru-Guzik
,
Phys. Chem. Chem. Phys.
11
,
3861
(
2009
).
40.
M.
Ceotto
,
S.
Atahan
,
G. F.
Tantardini
, and
A.
Aspuru-Guzik
,
J. Chem. Phys.
130
,
234113
(
2009
).
41.
M.
Buchholz
,
F.
Grossmann
, and
M.
Ceotto
,
J. Chem. Phys.
144
,
094102
(
2016
).
42.
M.
Buchholz
,
F.
Grossmann
, and
M.
Ceotto
,
J. Chem. Phys.
147
,
164110
(
2017
).
43.
F.
Gabas
,
R.
Conte
, and
M.
Ceotto
,
J. Chem. Theory Comput.
13
,
2378
(
2017
).
44.
F.
Gabas
,
G.
Di Liberto
,
R.
Conte
, and
M.
Ceotto
,
Chem. Sci.
9
,
7894
(
2018
).
45.
M.
Bonfanti
,
J.
Petersen
,
P.
Eisenbrandt
,
I.
Burghardt
, and
E.
Pollak
,
J. Chem. Theory Comput.
14
,
5310
(
2018
).
46.
M.
Micciarelli
,
F.
Gabas
,
R.
Conte
, and
M.
Ceotto
,
J. Chem. Phys.
150
,
184113
(
2019
).
47.
T. J.
Martínez
,
M.
Ben-Nun
, and
R. D.
Levine
,
J. Phys. Chem.
100
,
7884
(
1996
).
48.
B. F. E.
Curchod
and
T. J.
Martínez
,
Chem. Rev.
118
,
3305
(
2018
).
49.
D. V.
Makhov
,
C.
Symonds
,
S.
Fernandez-Alberti
, and
D. V.
Shalashilin
,
Chem. Phys.
493
,
200
(
2017
).
50.
M.
Šulc
,
H.
Hernández
,
T. J.
Martínez
, and
J.
Vaníček
,
J. Chem. Phys.
139
,
034112
(
2013
).
51.
G. A.
Worth
,
M. A.
Robb
, and
I.
Burghardt
,
Faraday Discuss.
127
,
307
(
2004
).
52.
G. W.
Richings
,
I.
Polyak
,
K. E.
Spinlove
,
G. A.
Worth
,
I.
Burghardt
, and
B.
Lasorne
,
Int. Rev. Phys. Chem.
34
,
269
(
2015
).
53.
I.
Polyak
,
G. W.
Richings
,
S.
Habershon
, and
P. J.
Knowles
,
J. Chem. Phys.
150
,
041101
(
2019
).
54.
K.-W.
Sun
,
M. F.
Gelin
,
V. Y.
Chernyak
, and
Y.
Zhao
,
J. Chem. Phys.
142
,
212448
(
2015
).
55.
N.
Zhou
,
L.
Chen
,
Z.
Huang
,
K.
Sun
,
Y.
Tanimura
, and
Y.
Zhao
,
J. Phys. Chem. A
120
,
1562
(
2016
).
56.
M.
Werther
and
F.
Großmann
,
Phys. Rev. B
101
,
174315
(
2020
).
57.
E. J.
Heller
,
J. Chem. Phys.
62
,
1544
(
1975
).
58.
M. F.
Gelin
,
D.
Egorova
, and
W.
Domcke
,
J. Chem. Phys.
131
,
124505
(
2009
).
59.
W. T.
Pollard
,
H. L.
Fragnito
,
J.-Y.
Bigot
,
C. V.
Shank
, and
R. A.
Mathies
,
Chem. Phys. Lett.
168
,
239
(
1990
).
60.
S. Y.
Lee
and
E. J.
Heller
,
J. Chem. Phys.
76
,
3035
(
1982
).
61.
A.
Patoz
,
T.
Begušić
, and
J.
Vaníček
,
J. Phys. Chem. Lett.
9
,
2367
(
2018
).
62.
T.
Begušić
,
A.
Patoz
,
M.
Šulc
, and
J.
Vaníček
,
Chem. Phys.
515
,
152
(
2018
).
63.
J.
Vaníček
and
T.
Begušić
, in
Molecular Spectroscopy and Quantum Dynamics
, edited by
R.
Marquardt
and
M.
Quack
(
Elsevier
,
2021
), pp.
199
229
.
64.
M.
Wehrle
,
S.
Oberli
, and
J.
Vaníček
,
J. Phys. Chem. A
119
,
5685
(
2015
).
65.
A.
Prlj
,
T.
Begušić
,
Z. T.
Zhang
,
G. C.
Fish
,
M.
Wehrle
,
T.
Zimmermann
,
S.
Choi
,
J.
Roulet
,
J.-E.
Moser
, and
J.
Vaníček
,
J. Chem. Theory Comput.
16
,
2617
(
2020
).
66.
M.
Wehrle
,
M.
Šulc
, and
J.
Vaníček
,
J. Chem. Phys.
140
,
244114
(
2014
).
67.
N. V.
Golubev
,
T.
Begušić
, and
J.
Vaníček
,
Phys. Rev. Lett.
125
,
083001
(
2020
).
68.
T.
Begušić
,
J.
Roulet
, and
J.
Vaníček
,
J. Chem. Phys.
149
,
244115
(
2018
).
69.
M. A.
Rohrdanz
and
J. A.
Cina
,
Mol. Phys.
104
,
1161
(
2006
).
70.
M.
Braun
,
H.
Metiu
, and
V.
Engel
,
J. Chem. Phys.
108
,
8983
(
1998
).
71.
F.
Grossmann
,
J. Chem. Phys.
125
,
014111
(
2006
).
72.
C.-M.
Goletz
and
F.
Grossmann
,
J. Chem. Phys.
130
,
244107
(
2009
).
73.
M.
Buchholz
,
C.-M.
Goletz
,
F.
Grossmann
,
B.
Schmidt
,
J.
Heyda
, and
P.
Jungwirth
,
J. Phys. Chem. A
116
,
11199
(
2012
).
74.
S.
Römer
,
M.
Ruckenbauer
, and
I.
Burghardt
,
J. Chem. Phys.
138
,
064106
(
2013
).
75.
P.
Eisenbrandt
,
M.
Ruckenbauer
,
S.
Römer
, and
I.
Burghardt
,
J. Chem. Phys.
149
,
174101
(
2018
).
76.
P.
Eisenbrandt
,
M.
Ruckenbauer
, and
I.
Burghardt
,
J. Chem. Phys.
149
,
174102
(
2018
).
77.
P. A.
Kovac
and
J. A.
Cina
,
J. Chem. Phys.
147
,
224112
(
2017
).
78.
D.
Picconi
,
J. A.
Cina
, and
I.
Burghardt
,
J. Chem. Phys.
150
,
064111
(
2019
).
79.
D.
Picconi
,
J. A.
Cina
, and
I.
Burghardt
,
J. Chem. Phys.
150
,
064112
(
2019
).
80.
D.
Picconi
and
I.
Burghardt
,
Faraday Discuss.
221
,
30
(
2019
).
81.
X.
Cheng
and
J. A.
Cina
,
J. Chem. Phys.
141
,
034113
(
2014
).
82.
W. T.
Pollard
and
R. A.
Mathies
,
Annu. Rev. Phys. Chem.
43
,
497
(
1992
).
83.
T. N.
Do
,
M. F.
Gelin
, and
H.-S.
Tan
,
J. Chem. Phys.
147
,
144103
(
2017
).
84.
P. M.
Morse
,
Phys. Rev.
34
,
57
(
1929
).
85.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
G. A.
Petersson
,
H.
Nakatsuji
,
X.
Li
,
M.
Caricato
,
A. V.
Marenich
,
J.
Bloino
,
B. G.
Janesko
,
R.
Gomperts
,
B.
Mennucci
,
H. P.
Hratchian
,
J. V.
Ortiz
,
A. F.
Izmaylov
,
J. L.
Sonnenberg
,
D.
Williams-Young
,
F.
Ding
,
F.
Lipparini
,
F.
Egidi
,
J.
Goings
,
B.
Peng
,
A.
Petrone
,
T.
Henderson
,
D.
Ranasinghe
,
V. G.
Zakrzewski
,
J.
Gao
,
N.
Rega
,
G.
Zheng
,
W.
Liang
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
K.
Throssell
,
J. A.
Montgomery
, Jr.
,
J. E.
Peralta
,
F.
Ogliaro
,
M. J.
Bearpark
,
J. J.
Heyd
,
E. N.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
T. A.
Keith
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A. P.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
J. M.
Millam
,
M.
Klene
,
C.
Adamo
,
R.
Cammi
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
O.
Farkas
,
J. B.
Foresman
, and
D. J.
Fox
, Gaussian 16, Revision C.01,
Gaussian, Inc.
,
Wallingford, CT
,
2016
.
86.
K.
Rajak
,
A.
Ghosh
, and
S.
Mahapatra
,
J. Chem. Phys.
148
,
054301
(
2018
).
87.
F. J.
Avila Ferrer
and
F.
Santoro
,
Phys. Chem. Chem. Phys.
14
,
13549
(
2012
).
88.
A.
Baiardi
,
J.
Bloino
, and
V.
Barone
,
J. Chem. Theory Comput.
9
,
4097
(
2013
).
89.
V.
Butkus
,
L.
Valkunas
, and
D.
Abramavicius
,
J. Chem. Phys.
137
,
044513
(
2012
).
90.
V.
Butkus
,
D.
Zigmantas
,
L.
Valkunas
, and
D.
Abramavicius
,
Chem. Phys. Lett.
545
,
40
(
2012
).
91.
D. B.
Turner
,
K. E.
Wilk
,
P. M. G.
Curmi
, and
G. D.
Scholes
,
J. Phys. Chem. Lett.
2
,
1904
(
2011
).
92.
S.
Karmakar
,
D. P.
Mukhopadhyay
, and
T.
Chakraborty
,
J. Chem. Phys.
142
,
184303
(
2015
).
93.
A.
Rohatgi
, WebPlotDigitizer, https://automeris.io/WebPlotDigitizer; accessed
19 October 2020
.
94.
A.
Nenov
,
S.
Mukamel
,
M.
Garavelli
, and
I.
Rivalta
,
J. Chem. Theory Comput.
11
,
3755
(
2015
).
95.
J.
Segarra-Martí
,
S.
Mukamel
,
M.
Garavelli
,
A.
Nenov
, and
I.
Rivalta
,
Top. Curr. Chem.
376
,
24
(
2018
).
96.
A.
Nenov
,
S. A.
Beccara
,
I.
Rivalta
,
G.
Cerullo
,
S.
Mukamel
, and
M.
Garavelli
,
ChemPhysChem
15
,
3282
(
2014
).
97.
T.
Begušić
and
J.
Vaníček
,
J. Chem. Phys.
153
,
024105
(
2020
).
98.
J.
Gao
,
N.
Li
, and
M.
Freindorf
,
J. Am. Chem. Soc.
118
,
4912
(
1996
).
99.
G.
Granucci
,
J. T.
Hynes
,
P.
Millié
, and
T.-H.
Tran-Thi
,
J. Am. Chem. Soc.
122
,
12243
(
2000
).
100.
D. M.
Rogers
and
J. D.
Hirst
,
J. Phys. Chem. A
107
,
11191
(
2003
).
101.
Z.
Lan
,
W.
Domcke
,
V.
Vallet
,
A. L.
Sobolewski
, and
S.
Mahapatra
,
J. Chem. Phys.
122
,
224315
(
2005
).
102.
L.
Zhang
,
G. H.
Peslherbe
, and
H. M.
Muchall
,
Photochem. Photobiol.
82
,
324
(
2006
).
103.
O. P. J.
Vieuxmaire
,
Z.
Lan
,
A. L.
Sobolewski
, and
W.
Domcke
,
J. Chem. Phys.
129
,
224307
(
2008
).
104.
S.-s.
Kim
,
M.
Kim
, and
H.
Kang
,
Bull. Korean Chem. Soc.
30
,
1481
(
2009
).
105.
K. R.
Yang
,
X.
Xu
,
J.
Zheng
, and
D. G.
Truhlar
,
Chem. Sci.
5
,
4661
(
2014
).

Supplementary Material