A widely used strategy for simulating the charge transfer between donor and acceptor electronic states in an all-atom anharmonic condensed-phase system is based on invoking linear response theory to describe the system in terms of an effective spin-boson model Hamiltonian. Extending this strategy to photoinduced charge transfer processes requires also taking into consideration the ground electronic state in addition to the excited donor and acceptor electronic states. In this paper, we revisit the problem of describing such nonequilibrium processes in terms of an effective three-state harmonic model. We do so within the framework of nonequilibrium Fermi’s golden rule (NE-FGR) in the context of photoinduced charge transfer in the carotenoid–porphyrin–C60 (CPC60) molecular triad dissolved in explicit tetrahydrofuran (THF). To this end, we consider different ways for obtaining a three-state harmonic model from the equilibrium autocorrelation functions of the donor–acceptor, donor–ground, and acceptor–ground energy gaps, as obtained from all-atom molecular dynamics simulations of the CPC60/THF system. The quantum-mechanically exact time-dependent NE-FGR rate coefficients for two different charge transfer processes in two different triad conformations are then calculated using the effective three-state model Hamiltonians as well as a hierarchy of more approximate expressions that lead to the instantaneous Marcus theory limit. Our results show that the photoinduced charge transfer in CPC60/THF can be described accurately by the effective harmonic three-state models and that nuclear quantum effects are small in this system.

Photoinduced charge transfer (CT) plays an important role in many systems of practical interest, for example, solar energy conversion in photosynthesis and organic photovoltaic (OPV) devices.1–11 The ability to accurately simulate photoinduced CT dynamics could therefore offer invaluable insights toward the discovery of more efficient artificial energy-conversion materials. Many quantum dynamical methods have been proposed to study electronic transitions in the condensed phase, such as the semiclassical initial value representation,12–14 mean-field Ehrenfest,15–18 fewest switches surface hopping,19–21 mixed quantum-classical Liouville,22–25 generalized quantum master equation,26–35 and hierarchical equations of motion,36–41 to name a few. Despite these advances, simulating CT dynamics in realistic complex systems rigorously and accurately remains highly challenging.42–49 

Importantly, Marcus theory50–54 cannot account for the nonequilibrium nature of photoinduced CT, which is due to the fact that the initial state after photoexcitation from the ground state to the excited donor state corresponds to equilibrium on the ground state potential energy surface (PES) rather than on the donor state PES.55–60 The nonequilibrium nature of the initial state can have a significant effect on the CT dynamics when the timescale of nuclear relaxation to thermal equilibrium on the donor PES is comparable to or longer than the timescale of the CT.

Accounting for the above-mentioned nonequilibrium effects is made possible by the nonequilibrium Fermi’s golden rule (NE-FGR).61–65 Combined with the linearized semiclassical (LSC) approximation, the NE-FGR can also be applied to complex condensed-phase systems described by all-atom anharmonic Hamiltonians.66 Furthermore, starting from the LSC NE-FGR, a hierarchy of more approximate methods can be derived, which ultimately leads to the instantaneous Marcus theory (IMT), where the time-dependent CT rate coefficient is given by a Marcus-like expression with explicitly time-dependent reorganization energy and reaction free energy.66 The time-dependent IMT rate coefficient is expected to be accurate when (1) the nuclear degrees of freedom (DOF) can be treated as classical, (2) the nuclear motion occurs on a timescale slower than the electronic dephasing time, and (3) the time-dependent distribution of the donor–acceptor energy gap maintains its Gaussian nature.

An alternative and widely used approach relies on mapping the all-atom anharmonic Hamiltonian onto an effective harmonic model Hamiltonian, such as the popular spin-boson model for two-state systems.44–47,67–70 The mapping onto the spin-boson model typically relies on invoking linear response theory and is based on the donor–acceptor energy gap time correlation function (TCF) as obtained from a classical molecular dynamics (MD) simulation of the all-atom anharmonic system.71,72 Such mapping onto an effective spin-boson model Hamiltonian72 has been recently found to lead to remarkably accurate equilibrium FGR (E-FGR) CT rate constants in the case of the carotenoid–porphyrin–C60 (CPC60) molecular triad11,57,73–76 dissolved in an explicit tetrahydrofuran (THF) solvent.77–81 The main advantage of using the effective harmonic model is that once the mapping is established, the quantum-mechanically exact expression for both E-FGR and NE-FGR and the corresponding hierarchy of approximations are known in closed form.42,64

In this paper, we go beyond the two-state spin-boson Hamiltonian to construct effective three-state harmonic Hamiltonians for simulating NE-FGR CT rate coefficients as well as a progression of increasingly more approximate expressions, which ultimately lead to IMT. We benchmark the results obtained for the harmonic model against the all-atom predictions in the context of the following two photoinduced CT processes in the bent and linear CPC60 triad conformations in explicit THF (see Fig. 1). The triad is initially equilibrated on the ground (G) state, CPC60, before it is photoexcited impulsively to the P-localized ππ* state, CP*C60. Following the impulsive photoexcitation at time 0, there could be a nonradiative transition to the excited P-to-C60 CT state, CP+C60, which is denoted as CT1,

CPC60(G)hνCP*C60(ππ*)CP+C60(CT1),
(1)

or to the excited C-to-C60 charge separated state, C+PC60, which is denoted as CT2,

CPC60(G)hνCP*C60(ππ*)C+PC60(CT2).
(2)
FIG. 1.

Two characteristic conformations of the carotenoid–porphyrin–C60 (CPC60) molecular triad in explicit THF, (a) the bent triad conformation and (b) the linear triad conformation, and schematic representations of the potential energy surfaces involved in the photoinduced charge transfer processes in (c) the bent triad and (d) the linear triad, where vertical photoexcitation brings it from the ground state to the ππ* state (black arrow), and then, during the nuclear relaxation on the ππ* state (orange arrow), electronic transition to CT1 or CT2 states occurs. The molecular structure visualization was generated using the Visual Molecular Dynamics (VMD) package.82 

FIG. 1.

Two characteristic conformations of the carotenoid–porphyrin–C60 (CPC60) molecular triad in explicit THF, (a) the bent triad conformation and (b) the linear triad conformation, and schematic representations of the potential energy surfaces involved in the photoinduced charge transfer processes in (c) the bent triad and (d) the linear triad, where vertical photoexcitation brings it from the ground state to the ππ* state (black arrow), and then, during the nuclear relaxation on the ππ* state (orange arrow), electronic transition to CT1 or CT2 states occurs. The molecular structure visualization was generated using the Visual Molecular Dynamics (VMD) package.82 

Close modal

In what follows, we construct unique three-state harmonic Hamiltonians for each of the two conformations and photoinduced CT processes. It should be noted that mapping a three-state all-atom anharmonic system onto an effective three-state harmonic system is not as straightforward as the two-state case.83,84 This is because, as pointed out by Cho and Silbey in Ref. 62, the ground, donor, and acceptor PESs are not independent, and therefore, one needs to consider the relations between each pair of the PESs consistently. In particular, three different energy gap TCFs (donor–acceptor, donor–ground, and acceptor–ground) are needed in order to construct the harmonic model. In this paper, we propose and test several routes for achieving the consistency requirement by satisfying constraints imposed by various reorganization energies. Our previous LSC NE-FGR and IMT calculations showed a significant nonequilibrium transient effect in the bent triad conformation that enhanced the porphyrin-to-C60 [Eq. (1)] CT rate coefficient by a factor of ∼40.66 Therefore, testing the ability of the harmonic three-state Hamiltonians to accurately reproduce this effect serves as an important benchmark. The fact that the effective three-state harmonic model makes it easy to evaluate the fully quantum mechanical NE-FGR also allows us to assess the importance of nuclear quantum effects in this system.72 

The remainder of this paper is organized as follows. Section II summarizes the NE-FGR CT rate coefficient and the corresponding hierarchy of LSC-based approximations that leads to the IMT. Section III presents three different ways to construct three-state harmonic models for the nonequilibrium photoinduced CT process (models 1–3) as well as a strategy for calculating the IMT CT rate coefficients directly from the three spectral densities associated with the three above-mentioned TCFs. The all-atom anharmonic model for the CPC60 molecular triad dissolved in liquid THF and MD simulation techniques are described in Sec. IV. The results and discussion are reported in Sec. V. The conclusions and outlook are provided in Sec. VI. The derivations of the IMT rate coefficient for the three-state model and additional supporting tables and figures are provided in the supplementary material.

Consider a charge transfer system with the Hamiltonian

Ĥ=ĤD|DD|+ĤA|AA|+ΓDA(|DA|+|AD|).
(3)

Here, |D⟩ and |A⟩ represent the diabatic donor and acceptor electronic states, respectively, and ĤD/A are the corresponding nuclear Hamiltonians, ĤD/A=P^2/2+VD/A(R^), where R^=R^1,,R^N and P^=P^1,,P^N are the mass-weighted nuclear coordinates and momenta, VD/A(R^) are the donor/acceptor PESs, and ΓDA is the electronic coupling coefficient that is assumed to be constant within the Condon approximation.

We assume that the system starts out at the donor electronic state with the initial nuclear DOF described by the ground state equilibrium nuclear density operator ρ^G=eβĤG/Trn[eβĤG]. Here, ĤG=P^2/2+VG(R^) is the ground state nuclear Hamiltonian, VG(R^) is the ground state PES, β = 1/kBT is the inverse temperature, Trn(·) denotes the trace over the nuclear Hilbert space, and Trn Tre denotes the trace over both the nuclear and electronic Hilbert spaces. The donor state population, PD(t), at a later time t is given by

PD(t)=TrnTreeiĤt/ρ^GeiĤt/|DD|.
(4)

Applying second-order perturbation theory then leads to the following approximate NE-FGR expression for the donor state population:61 

PD(t)exp0tdtkt,
(5)

where the time-dependent rate coefficient is defined as

kt=22Re0tdτCt,τ.
(6)

Here,

Ct,τ=ΓDA2TrneiĤDt/ρ^GeiĤDt/eiĤAτ/eiĤDτ/.
(7)

It should be noted that the NE-FGR expression in Eq. (7) is fully quantum-mechanical, which can still be difficult to evaluate in complex systems.

We also note in passing that assuming that the initial state of the nuclear DOF corresponds to equilibrium on the donor state PES, ρ^D=eβĤD/Trn[eβĤD], C(t, τ) reduces to

Cτ=ΓDA2Trnρ^DeiĤAτ/eiĤDτ/.
(8)

Also assuming that the CT happens on a timescale longer than the electronic dephasing time, which is defined by the lifetime of Cτ, then leads to the time-dependent rate coefficient k(t) turning into the time-independent E-FGR CT rate constant,42 

kDA=22Re0dτC(τ).
(9)

Under those conditions, the donor state population could be estimated by the exponential decay PD(t) = exp(−kDAt). It should be noted that the rate constant kDA is the long-time asymptotic limit of the time-dependent rate coefficient, k(t).64,66 Accounting for the nuclear initial preparation effect is important when k(t) and kDA are significantly different and the timescale for reaching thermal equilibrium on the donor PES is comparable to or longer than the timescale of CT, kDA1.

The linearized semiclassical (LSC) approximation provides a computationally feasible route for calculating NE-FGR rates. The main idea behind LSC is to express Ct,τ in terms of the real-time path integral and then apply the linearization approximation to the differences between the forward and backward paths, which filters out the highly oscillatory contributions and leads us to a classical-like procedure for calculating the time-dependent CT rate coefficient. The hierarchy of LSC-based approximations for C(t, τ) is summarized as follows:64 

CW-AV(t,τ)=|ΓDA|2dR0dP0ρG,W(R0,P0)×expittτdtU(Rtav),
(10)
CW-0(t,τ)=|ΓDA|2dR0dP0ρG,W(R0,P0)expiU(Rt)τ,
(11)
CC-AV(t,τ)=|ΓDA|2dR0dP0ρG,Cl(R0,P0)×expittτdtU(Rtav),
(12)
CC-D(t,τ)=|ΓDA|2dR0dP0ρG,Cl(R0,P0)×expittτdtU(RtD),
(13)
CC-0(t,τ)=|ΓDA|2dR0dP0ρG,Cl(R0,P0)expiU(Rt)τ.
(14)

Here, “W” denotes the semiclassical Wigner distribution for initial sampling of the nuclear DOF,

ρG,W(R,P)=12πNdZexpiZPRZ2ρ^GR+Z2,
(15)

and “C” denotes the classical nuclear sampling using the corresponding classical distribution ρG,Cl(R, P). Following the initial sampling of (R0, P0) via Wigner or classical nuclear densities, the system is propagated on the donor PES forward in time over time interval t, arriving at the configuration Rt, and then backward in time from time t to tτ (τ = 0 → t), integrating the donor–acceptor energy gap U(R) = UDA(R) = VD(R) − VA(R) in the phase factor for each (t, τ). “AV” denotes that the dynamics over τ is propagated on the average PES, Vav(R) = (VD(R) + VA(R))/2, leading to trajectories of Rτav, “D” denotes that the τ-dynamics is on the donor surface, leading to trajectories of RτD, and “0” denotes no τ-dynamics beyond the t-relaxation on the donor PES until Rt. The W-AV level is the original LSC approximation for the NE-FGR.

From the C-0 approximation for the NE-FGR, one can derive the instantaneous Marcus theory (IMT) expression of k(t),66 

kIMT(t)=|ΓDA|22πσt2expUt¯22σt2,
(16)

where we assume that the instantaneous distribution of the energy gap Ut is Gaussian with a time-dependent mean Ut¯ and the corresponding variance σt2=Ut2¯Ut¯2 at time t after the photoexcitation so that Prob(Ut)=1/2πσt2exp(UtUt¯)2/(2σt2). Here, Ut¯ and σt2 can be calculated using energies obtained from all-atom molecular dynamics simulations.66 The IMT CT rate coefficient kIMT(t) can also be expressed equivalently in terms of time-dependent reaction free energy ΔE(t) and reorganization energy Er(t),

Er(t)=σt22kBT=ΔE(t)Ut¯.
(17)

In addition, the hierarchy of LSC-based approximations for the E-FGR rate constants corresponds to replacing the t-relaxation on the donor PES with equilibrium sampling on the donor PES, so C(t, τ) in Eqs. (10)(14) would reduce to W-AV, W-0, C-AV, C-D, and C-0 levels of approximation for E-FGR C(τ) in Eq. (8), respectively. Assuming the second-order cumulant approximation for U at the C-0 level of approximation, the classical Marcus rate constant is thus given by42 

kDAM=|ΓDA|2πkBTErexp(ΔE+Er)24kBTEr.
(18)

The corresponding reorganization energy and reaction free energy satisfy the relation Er=σeq2/(2kBT)=ΔEU and the activation energy Ea=kBTU2/(2σeq2), where ⟨·⟩ is the equilibrium average on the donor state PES and the equilibrium variance of the energy gap σeq2=U2U2.

The spin-boson model has been widely used for simulating CT dynamics in the condensed phase.44–47,67–70 The prototypical spin-boson Hamiltonian is given by47 

Ĥ=ΓDAσ^x+ωDA2σ^z+j=1NP^j22+12ωj2R^j2cjR^jσ^z.
(19)

Here, σ^x=|DA|+|AD| and σ^z=|DD||AA| are the Pauli operators; ΓDA is the electronic coupling coefficient; ΔE = −ℏωDA is the donor-to-acceptor reaction free energy; {R^j,P^j,ωj|j=1,,N} are the mass-weighted coordinates, momenta, and frequencies associated with the nuclear normal modes, respectively; and {cj} are the coupling coefficients between the electronic DOF and nuclear normal modes. The frequencies and the coupling coefficients, {ωj, cj}, are often given in terms of spectral density, which is defined by

J(ω)=π2j=1Ncj2ωjδ(ωωj).
(20)

Next, we outline the procedure typically followed for mapping an all-atom anharmonic Hamiltonian such as that of the triad in liquid THF solution onto the spin-boson Hamiltonian. It starts out by obtaining the donor–acceptor energy gap time correlation function, CUU(t), from an equilibrium MD simulation on the donor PES,72 

CUU(t)=U(t)U(0)U2,
(21)

and using the following relation to obtain the spectral density (note that in some literature,71 the TCF of half energy gap is used, so we see a factor of 1/4 difference):

J(ω)=βω40dtCUU(t)cos(ωt).
(22)

It is noted that within this framework, the reorganization energy can be identified as being given in terms of the donor–acceptor energy gap variance, CUU(0)=σeq2,

Er=CUU(0)2kBT=4π0dωJ(ω)ω.
(23)

Obtaining a discrete representation of the nuclear DOF in terms of N normal modes (N ≫ 1) from the continuous spectral density function in Eq. (22) is based on solving the following equation for the frequency ωj of the jth mode (j = 1, 2, …, N):71 

2NωjπCUU(0)0dtCUU(t)ωjtsin(ωjt)=j12.
(24)

It is worth noting that following this procedure, the obtained discrete frequency intervals correspond to equal fractions of the reorganization energy. Once the frequency of the jth mode is determined, the corresponding coupling coefficients, cj, are obtained via the relation

cj=Er2Nωj(j=1,2,,N).
(25)

Next, we consider generalizing this procedure to the case of a nonequilibrium CT process that takes place in a three-state system. To this end, we define the following effective harmonic Hamiltonian:

Ĥ=ĤG|GG|+ĤD|DD|+ĤA|AA|+ΓDA(|DA|+|AD|),
(26)

where the nuclear Hamiltonians for the donor (D), acceptor (A), and ground (G) electronic states are given by

ĤD=j=1NP^j22+j=1N12ωj2R^j2+ωDA,ĤA=j=1NP^j22+j=1N12ωj2R^jRjeq2,ĤG=j=1NP^j22+j=1N12ωj2R^j+Sj2+εG.
(27)

Here, the energy difference between the donor and acceptor states, ΔE = −ℏωDA, corresponds to the CT reaction free energy, while the ground state energy, εG, is defined relative to the energy of the acceptor state at its equilibrium geometry. It should be noted that the photoinduced CT rate is independent of the value of εG since we assume that the photoexcitation corresponds to a vertical excitation from the ground to the donor state. The displacement between the donor and acceptor equilibrium geometries along the jth mode (j = 1, 2, …, N) is given by [see Eqs. (25) and (19) for the definition of cj]

Rjeq=2cjωj2=2ErN1ωj.
(28)

Thus, the reorganization energy between donor and acceptor states is given by

Er=j=1N2cj2ωj2=j=1N12ωj2(Rjeq)2.
(29)

Expressed in terms of {Rjeq} instead of {cj}, the spectral density can be written in the following form:

J(ω)=π8j=1Nωj3(Rjeq)2δ(ωωj).
(30)

Finally, {Sj} correspond to the displacements between the ground and donor equilibrium geometries along the normal mode coordinates, which is responsible for the nonequilibrium initial preparation of the CT process (see below).

For the three-state harmonic model defined in Eqs. (26) and (27), we can express the NE-FGR hierarchy of approximations in the closed form [it should be noted that the LSC approximation (W-AV) coincides with the exact quantum-mechanical expression in this case]64 

Cexact/W-AV(t,τ)=ΓDA2expiωDAτj=1Nωj(Rjeq)22×cothβωj21cos(ωjτ)+isin(ωjτ)j=1NiωjRjeqSj[sin(ωjt)+sin(ωjτωjt)],
(31)
CW-0(t,τ)=ΓDA2expiωDAτj=1Nωj(Rjeq)22×cothβωj2ωj2τ22+iωjτj=1NiωjRjeqSjcos(ωjt)ωjτ,
(32)
CC-AV(t,τ)=ΓDA2expiωDAτj=1Nωj(Rjeq)22×2βωj1cos(ωjτ)+isin(ωjτ)j=1NiωjRjeqSj[sin(ωjt)+sin(ωjτωjt)],
(33)
CC-D(t,τ)=ΓDA2expiωDAτj=1Nωj(Rjeq)22×2βωj1cos(ωjτ)+iωjτj=1NiωjRjeqSj[sin(ωjt)+sin(ωjτωjt)],
(34)
CC-0(t,τ)=ΓDA2expiωDAτj=1Nωj(Rjeq)22×ωjτ2β+iωjτj=1NiωjRjeqSjcos(ωjt)ωjτ.
(35)

The time-dependent IMT parameters Ut¯ and σt2 in Eq. (16) for the three-state harmonic model can be expressed in the closed form (the derivation can be found in the supplementary material)

Ut¯=ωDAj=1N12ωj2(Rjeq)2+2RjeqSjcos(ωjt)=U+βC1L(t),
(36)
σt2=σ02=1βj=1Nωj2(Rjeq)2.
(37)

Here, C1L(t) is the cross correlation function of the donor–acceptor and donor–ground energy gaps,

C1L(t)=UDA(t)UDG(0)UDAUDG=1βj=1Nωj2RjeqSjcos(ωjt).
(38)

It should be noted that σt2=σ02 for the harmonic model, which implies that any significant time-dependence of σt2 is a signature for deviations from this model. It should also be noted that the same expression can be derived from linear response theory,85 which is exact in the case of the three-state harmonic model (see the supplementary material).

We now consider the construction of three-state harmonic models of the form of Eqs. (26) and (27) from inputs obtained from MD simulations. To this end, and in analogy to the spin-boson case, we use as inputs the correlation functions of the three energy gaps (donor–acceptor, donor–ground, and acceptor–ground), which are calculated from classical equilibrium MD simulations on the donor PES,

CUUDA(t)=UDA(t)UDA(0)UDA2CUU(t),CUUDG(t)=UDG(t)UDG(0)UDG2,CUUAG(t)=UAG(t)UAG(0)UAG2.
(39)

Here, UXY = VXVY (X, Y = G, D, A) is the energy gap between X and Y PESs. Correspondingly, we can define three different reorganization energies (XY = DA, DG, AG),

ErXY=CUUXY(0)2kBT.
(40)

Following the two-state case, we chose to obtain the frequencies {ωj} of the nuclear modes from CUUDA(t) based on Eq. (24) and assume them to be the same for all three states. Once the frequencies of the nuclear modes are determined, we can turn to assigning values to the equilibrium geometry displacements, {Rjeq} and {Sj}. The equilibrium geometry displacements between the donor and acceptor states, {Rjeq}, can be unambiguously determined from ErDA=Er as in Eq. (28),

Rjeq=RjDA2ErDAN1ωj(j=1,,N).
(41)

Similarly, we also define the equilibrium shifts between DG and between AG,

RjDG2ErDGN1ωj(j=1,,N),
(42)
RjAG2ErAGN1ωj(j=1,,N).
(43)

One can also show that

ErDA=CUUDA(0)2kBT=j12ωj2(RjDA)2,
(44)
ErDG=CUUDG(0)2kBT=j12ωj2(RjDG)2,
(45)
ErAG=CUUAG(0)2kBT=j12ωj2(RjAG)2.
(46)

Now, we need to determine the initial shifts between the donor and ground states, {Sj}. From the nuclear Hamiltonians in Eq. (27), one would expect

ErDA=j12ωj2(Rjeq)2,
(47)
ErDG=j12ωj2Sj2,
(48)
ErAG=j12ωj2(Rjeq+Sj)2.
(49)

Equations (47) and (48) can be satisfied by using Eq. (41) and by choosing

|Sj|=RjDG(j=1,,N).
(50)

It should be noted that only the absolute values of {Sj} are determined by Eq. (48). If we choose the nonequilibrium shifts to be positive and Sj=RjDG (for all j = 1, …, N), the effective AG reorganization energy will be larger than the actual value ErAG,eff=j12ωj2(Rjeq+RjDG)2>ErAG=CUUAG(0)/(2kBT), which means that the displacements between the ground and acceptor minima are too large. It is important to note that the three reorganization energies are not entirely independent (e.g., UDG = UDA + UAG), so forcing the nonequilibrium shifts to be colinear with UDA, for example, will overestimate the nonequilibrium initial shifts. Therefore, the relative signs of {Rjeq,Sj} need to be chosen in a way that satisfies Eq. (49).

Figure 2 illustrates the effect of flipping the sign of Sj for a given Rjeq. The relative signs of Sj and Rjeq can clearly lead to either a positive or a negative contribution of the jth mode to ErAG by the amount of ΔErAG(j)=±2ωj2RjeqRjDG (see Table S3 of the supplementary material). Since ErAG=j12ωj2(Rjeq+Sj)2 can be satisfied for different selections of relative signs of Sj and Rjeq, satisfying Eq. (49) does not dictate a unique choice of the relative signs of {Rjeq,Sj}.

FIG. 2.

Flipping the sign of the minimum position of the jth mode on the ground state (G) with a positive Sj (black) to a negative Sj (gray) would decrease the effective reorganization energy between the acceptor state and the ground state ErAG, whereas the minimum positions of the donor state (D, blue) and the acceptor state (A, red) do not change. Here, ⟨U⟩, Er, ΔE, and εG are the energetic parameters of the jth mode contribution (model 1).

FIG. 2.

Flipping the sign of the minimum position of the jth mode on the ground state (G) with a positive Sj (black) to a negative Sj (gray) would decrease the effective reorganization energy between the acceptor state and the ground state ErAG, whereas the minimum positions of the donor state (D, blue) and the acceptor state (A, red) do not change. Here, ⟨U⟩, Er, ΔE, and εG are the energetic parameters of the jth mode contribution (model 1).

Close modal

In practice, we choose Rjeq=RjDA0(j=1,,N) and randomly flip the signs of {Sj} such that the deviations between ErAG,eff=j12ωj2(Rjeq+Sj)2 and ErAG=CUUAG(0)/(2kBT) are minimal. Alternatively, one can flip the signs of {Sj} relative to {Rjeq} evenly, namely, every few modes, until ErAG,effErAG [the two implementations turn out to give similar results for the system under consideration in this paper (see the supplementary material)]. In what follows, we will refer to this approach as model 1.

An alternative approach for determining {Sj} was inspired by the work of Cho and Silbey.62 The relations between the donor, acceptor, and ground PESs can be depicted in terms of the two-dimensional illustration in Fig. 3. Since the shifts between the equilibrium geometry of different pairs of PESs are proportional to the square root of the corresponding reorganization energy [see Eqs. (41)(43)], i.e., RjXYErXY, we expect the three reorganization energies to form a triangle with the three edges proportional to ErDA,ErDG,ErAG, and the angle θ between DA and DG edges is given by

cosθ=ErDA+ErDGErAG2ErDAErDG.
(51)
FIG. 3.

A schematic representation of the three-state model in two-dimensional space where qCT is the CT axis and qor is the axis orthogonal to the CT axis. The donor (D), acceptor (A), and ground (G) potential energy surfaces (PESs) fall onto the different locations forming a triangle, where θ is the angle between DA and DG PES minima. The effective nonequilibrium shifts {SjCT} are the projections of the donor–ground minimum displacements {RjDG} onto the donor–acceptor CT axis shown as the horizontal dashed line (models 2 and 3).

FIG. 3.

A schematic representation of the three-state model in two-dimensional space where qCT is the CT axis and qor is the axis orthogonal to the CT axis. The donor (D), acceptor (A), and ground (G) potential energy surfaces (PESs) fall onto the different locations forming a triangle, where θ is the angle between DA and DG PES minima. The effective nonequilibrium shifts {SjCT} are the projections of the donor–ground minimum displacements {RjDG} onto the donor–acceptor CT axis shown as the horizontal dashed line (models 2 and 3).

Close modal

As shown in Fig. 3, the DG edge or {RjDG} can be decomposed into a component along the CT axis, qCT (the horizontal dashed line along the donor–acceptor displacement direction), and an orthogonal component qor, which is perpendicular to it,

SjCT=RjDGcosθ(j=1,,N),
(52)
Sjor=RjDGsinθ(j=1,,N).
(53)

In this triangular case, the effective AG reorganization energy is given by

ErAG,eff=j=1N12ωj2Rjeq+SjCT2+Sjor2=j=1N12ωj2Rjeq2+RjDG22RjeqRjDGcosθ=ErDA+ErDG2cosθj=1N12ωj2RjeqRjDG=ErDA+ErDG2cosθj=1N12ωj22ErDAN1ωj2ErDGN1ωj=ErDA+ErDG2cosθErDAErDG=ErAG,
(54)

which agrees with ErAG=CUUAG(0)/(2kBT) from MD simulation. Since the initial shifts along the orthogonal axis are irrelevant to the CT process, it is reasonable to assume that the effective nonequilibrium shifts along the CT axis {Sj} in Eq. (27) are {RjDG} projected onto the CT axis SjCT, i.e.,

Sj=SjCT=RjDGcosθ(j=1,,N).
(55)

In what follows, we will refer to this approach as model 2.

Alternatively, we can extend the number of harmonic modes in model 2 from N to 2N and dividing them into two subgroups of N modes each based on the displacements in their equilibrium geometries [see Eq. (56)]. To this end, the nuclear Hamiltonians for the donor, acceptor, and ground electronic states are redefined as follows:

ĤD=j=12NP^j22+j=1N12ωj2R^j2+R^N+j2+ωDA,ĤA=j=12NP^j22+j=1N12ωj2R^jRjeq2+R^N+j2,ĤG=j=12NP^j22+j=1N12ωj2R^j+SjCT2+R^N+j+Sjor2+εG,

where {SjCT,Sjor} are given by Eqs. (52) and (53), respectively. In what follows, we will refer to this approach as model 3.

Finally, we note that the time-dependent IMT parameters Ut¯ and σt [see Eqs. (36) and (37)], which are necessary for calculating the IMT rate [see Eq. (16)], can be calculated directly from the spectral densities defined as follows:

JDA(ω)=βω40dtCUUDA(t)cos(ωt)=π8j=1Nωj3Rjeq2δ(ωωj),JDG(ω)=βω40dtCUUDG(t)cos(ωt)=π8j=1Nωj3Sj2δ(ωωj),JAG(ω)=βω40dtCUUAG(t)cos(ωt)=π8j=1Nωj3Rjeq+Sj2δ(ωωj).
(56)

We thereby bypass the need to assign values to the equilibrium geometry displacements, {Rjeq,Sj}. Using the spectral densities, we express the TCF in Eq. (38) as follows:

C1L(t)=j=1N1βωj2RjeqSjcos(ωjt)=8πβ0dωωπ8j=1Nω3SjRjeqδ(ωωj)cos(ωt)=4πβ0dωωJDA(ω)+JDG(ω)JAG(ω)cos(ωt).
(57)

The IMT parameters for Eq. (16) are thus given by

Ut¯=4π0dωωJDA(ω)+JDG(ω)JAG(ω)cos(ωt)+U,
(58)
σt2=σ02=8πβ0dωωJDA(ω).
(59)

In this section, we describe the all-atom simulations that were used to provide classical MD inputs for constructing the three-state harmonic models. The all-atom model of the CPC60 triad dissolved in explicit THF was adopted from Ref. 81. For the bent and linear triad conformations, four electronic-state-dependent force fields were constructed for the ground, ππ*, CT1, and CT2 states. The overall PESs of the four electronic states of the triad differ with respect to the excitation energies and the triad’s atomic charges. The non-electrostatic interactions are assumed to be the same for different electronic states. The excitation energies, atomic charges, and electronic coupling coefficients for this system were computed using time-dependent density functional theory (TDDFT) with the range-separated hybrid (RSH) Baer–Neuhauser–Livshits (BNL) functional86–88 using the Q-Chem 4 program package.89 The electronic coupling coefficients between the electronically excited states were calculated via the fragment charge difference (FCD) method.90 The overall PES in the αth excited state is given by

Vα(R)=VαMD(R)+Eα(rTriad)Vα,T(rTriad)=VαMD(R)+Wα.
(60)

Here, Eα(rTriad) represents the αth excited state energy with respect to the ground state for the gas-phase triad, which were obtained from electronic structure calculations in the characteristic bent and linear triad geometries, i.e., rTriad. It is crucial to include the excitation energy in the total potential energy because classical force fields can only describe the interactions or forces between all atoms but not the relative energy between states (e.g., from ground → ππ*, the atomic charges do not change too much, but still there is a sizable energy difference due to the electronic excitation). However, simply adding up the MD potential energy of the entire system VαMD(R) and the excitation energy Eα(rTriad) of the triad would double count the energy of the triad intramolecular interactions, so we need to subtract this double-counted contribution Vα,T(rTriad) evaluated with the force fields where only the triad is present. Finally, we arrive at the energy correction Wα = Eα(rTriad) − Vα,T(rTriad) to the MD potential energy of the αth excited state [Eq. (60)]. Table S1 shows the values of Eα, Vα,T, and Wα for the ππ*, CT1, and CT2 states in the bent and linear conformations. For details of force field parameters and its agreement with experiment on CT rate constants, we refer the reader to Ref. 81.

MD simulations were performed for a system containing one triad molecule and 6741 THF molecules in a 100 Å ×100 Å ×100 Å periodic cubic box, performed using the AMBER 18 package.91 The van der Waals cutoff radius was set to be 9 Å. The SHAKE algorithm92 was used to constrain all covalent bonds involving hydrogen atoms. Particle mesh Ewald summation was used to calculate the electrostatic interactions.93 To maintain the linear triad conformation, the end-to-end distance was constrained at 49.6 Å using a steep harmonic potential with a force constant of 100 kcalmol1Å−2. No constraint was used in the simulation of a flexible bent triad in flexible THF. The MD time step was chosen as δt = 1.0 fs. The system was equilibrated on the ππ* state at 300 K using a Langevin thermostat with a collision frequency of 1.0 ps−1.

Two sets of all-atom equilibrium MD simulations comprising of 200 independent trajectories of length 500 ps each were generated under the NVE ensemble after equilibration under the NVT ensemble, one set propagated on the ground state PES and another on the donor state PES. The initial positions and velocities for the NVE production runs were sampled every 5 ps from an NVT ensemble at a temperature of 300 K equilibrated on either the ground state PES or the excited ππ* state PES, followed by a 50 ps re-equilibration under constant NVE conditions. Configurations were sampled every 5 fs in these 200 NVE trajectories, which means a total sum of 2 × 107 configurations were sampled from equilibrium NVE trajectories amounting to a total propagation duration of 100 ns in each case of different triad conformations. The potential energies on every electronic-state PES were recalculated from these trajectories.

In the all-atom NEMD simulations employed as benchmark, a total of 20 000 nonequilibrium NVE trajectories were initially sampled from an equilibrated ground state in the NVT ensemble at 300 K and propagated for 4 ps under constant NVE conditions with configuration recorded every time step, i.e., 1.0 fs. Thus, a total sum of 8 × 107 nonequilibrium configurations were averaged over to obtain the IMT rate coefficient. The error bars for k(t) of the all-atom NEMD simulations were obtained by averaging five sets of 4000 NVE trajectories.66 

In the construction of model 1, we selected the number of flipped modes, Nf, such that the deviation of the effective reorganization energy ErAG,eff from the accurate value ErAG is minimized,

argminNfErAG,eff(Nf)ErAG.
(61)

In this section, we present the NE-FGR CT rate coefficients obtained using the three three-state harmonic models (models 1–3) as well as IMT via J(ω). The results obtained based on the three-state harmonic models are compared to the IMT rate coefficients calculated via the LSC method from all-atom NEMD simulations. We construct the three-state models for two separate photoinduced CT pathways in the CPC60 triad dissolved in THF. More specifically, we start with a vertical photoexcitation from the equilibrium ground state to the donor state, which is then followed by an electronic transition to one of the two acceptor states [CT1 or CT2, see Eqs. (1) and (2)]. In what follows, we refer to those pathways as ππ* → CT1 and ππ* → CT2. Different sets of results are shown for the two characteristic triad conformations, i.e., bent and linear.

We start out with Table I that summarizes the reorganization energies ErDA, ErDG, and ErAG for the ππ* → CT1 and ππ* → CT2 transitions in the bent and linear triad conformations, as obtained from the all-atom MD simulations. Importantly, any reasonable three-state harmonic model needs to be consistent with those three reorganization energies, ErDA, ErDG, and ErAG. The inputs from MD simulations for constructing the three-state harmonic models are the time autocorrelation functions of the energy gaps between different PESs in the bent and linear triad conformations, which are shown in Fig. 4. It is seen that most of the normalized CUUXY(t) have similar exponential decay profiles with lifetime τ1/e ∼ 1.0 ps except for CUUG,ππ*(t) in the linear triad case that exhibits a fast and non-monotonic decay but with a very small magnitude. The initial values CUUXY(0) are proportional to the corresponding reorganization energies, i.e., CUUXY(0)=2kBTErXY [see Eqs. (47)(49) and reorganization energies in Table I]. The fact that the charge distribution on the ground state is similar to that on the ππ* state results in the smaller amplitude of CUUG,ππ*(t) (black lines). In contrast, since the difference in charge distributions between the ground/ππ* state and the CT2 state is the largest; this leads to the largest amplitude of CUUG/ππ*,CT2(t) (orange and blue lines). As expected, the amplitude of CUUG/ππ*,CT1(t) is intermediate (cyan and red lines). The noticeable difference between the bent and linear triad conformations reflects the fact that the solute–solvent interactions depend on the triad conformation.

TABLE I.

Reorganization energies between different pairs of potential energy surfaces (ErXY in eV) for ππ* → CT1 and ππ* → CT2 transitions in the bent and linear triad conformations.

TransitionConf.ErDAErDGErAG
ππ* → CT1 Bent 0.533 0.0914 0.924 
 Linear 0.687 0.0109 0.666 
ππ* → CT2 Bent 1.46 0.0914 1.75 
 Linear 1.71 0.0109 1.73 
TransitionConf.ErDAErDGErAG
ππ* → CT1 Bent 0.533 0.0914 0.924 
 Linear 0.687 0.0109 0.666 
ππ* → CT2 Bent 1.46 0.0914 1.75 
 Linear 1.71 0.0109 1.73 
FIG. 4.

Un-normalized (top) and normalized (bottom) time correlation functions of the energy gap between different potential energy surfaces, CUUXY(t) for ππ* − G (black), CT1 − G (cyan), CT2 − G (orange), ππ* − CT1 (red), and ππ* − CT2 (blue) in the bent (solid line) and linear (dashed line) triad conformations. The inset of the top panel is CUUXY(t) for the ππ* − G energy gap in the linear triad conformation, and the inset of the bottom panel is for the ππ* − G energy gap in the bent and linear triad conformation.

FIG. 4.

Un-normalized (top) and normalized (bottom) time correlation functions of the energy gap between different potential energy surfaces, CUUXY(t) for ππ* − G (black), CT1 − G (cyan), CT2 − G (orange), ππ* − CT1 (red), and ππ* − CT2 (blue) in the bent (solid line) and linear (dashed line) triad conformations. The inset of the top panel is CUUXY(t) for the ππ* − G energy gap in the linear triad conformation, and the inset of the bottom panel is for the ππ* − G energy gap in the bent and linear triad conformation.

Close modal

Next, Figs. 5 and 6 show the three spectral densities JDA(ω), JAG(ω), and JDG(ω), as obtained via Eq. (56) using energy gap TCFs CUUDA(t), CUUAG(t), and CUUDG(t), in the bent and linear triad conformations (see Fig. 4). Due to the fact that the decay profiles of the normalized CUUXY(t)/CUUXY(0) in the bent triad are relatively insensitive to the PESs, the spectral densities that correspond to different pairs of PESs have a similar shape and only differ by a scaling factor that is proportional to CUUXY(0)=σeq2=2kBTErXY. This justifies using the frequencies obtained from JDA(ω) to construct the three-state models (shown as cyan vertical lines in Fig. 5). However, in the linear triad conformation, the spectral densities JDA(ω) and JAG(ω) have the same profile for the same transition except for a scaling factor, while the spectral density JDG(ω) shows different low frequency distributions, but it is two orders of magnitude smaller than the other spectral densities (bottom of Fig. 6). Likewise, we choose to use the frequencies from discretizing JDA(ω) to construct the three-state models of the linear triad (shown as cyan vertical lines in Fig. 6). From our previous all-atom simulations, the nonequilibrium effect in the linear triad is negligible, which is reflected by the very small magnitude in JDG(ω) that would in turn give rise to rather small equilibrium geometry shifts between the donor and the ground PESs in the three-state models.

FIG. 5.

Three spectral densities JDA(ω), JAG(ω), and JDG(ω), calculated via Eq. (56) using energy gap time correlation functions CUUDA(t), CUUAG(t), and CUUDG(t), respectively. The spectral densities are for the bent triad conformation undergoing donor-to-acceptor transition ππ* → CT1 (left) and ππ* → CT2 (right). The cyan vertical lines indicate the N = 500 mode frequencies extracted from JDA(ω).

FIG. 5.

Three spectral densities JDA(ω), JAG(ω), and JDG(ω), calculated via Eq. (56) using energy gap time correlation functions CUUDA(t), CUUAG(t), and CUUDG(t), respectively. The spectral densities are for the bent triad conformation undergoing donor-to-acceptor transition ππ* → CT1 (left) and ππ* → CT2 (right). The cyan vertical lines indicate the N = 500 mode frequencies extracted from JDA(ω).

Close modal
FIG. 6.

Three spectral densities JDA(ω), JAG(ω), and JDG(ω), calculated via Eq. (56) using energy gap time correlation functions CUUDA(t), CUUAG(t), and CUUDG(t), respectively. The spectral densities are for the linear triad conformation undergoing donor-to-acceptor transition ππ* → CT1 (left) and ππ* → CT2 (right). The cyan vertical lines indicate the N = 500 mode frequencies extracted from JDA(ω).

FIG. 6.

Three spectral densities JDA(ω), JAG(ω), and JDG(ω), calculated via Eq. (56) using energy gap time correlation functions CUUDA(t), CUUAG(t), and CUUDG(t), respectively. The spectral densities are for the linear triad conformation undergoing donor-to-acceptor transition ππ* → CT1 (left) and ππ* → CT2 (right). The cyan vertical lines indicate the N = 500 mode frequencies extracted from JDA(ω).

Close modal

We present our main results in Fig. 7 for the bent triad and Fig. 8 for the linear triad, which show the time-dependent CT rate coefficients for both ππ* → CT1 and ππ* → CT2 transitions computed with NE-FGR and the hierarchy of semiclassical approximations including the IMT using models 1–3 and the IMT via spectral densities [Eqs. (58) and (59)], compared with the IMT rate coefficient obtained with all-atom NEMD simulations. The general observation is that all the three models and IMT via J(ω) can generate accurate time-dependent CT rate coefficients for all the transitions and triad conformations under consideration. We also note that the nuclear quantum effects are small in that the fully quantum-mechanical NE-FGR rate (equivalent to W-AV) is only slightly faster than the IMT values (within a factor of 2). The initial CT rate coefficients (t = 0) agree with the corresponding E-FGR CT rate constants obtained with the equilibrium ground state parameters, and the long-time CT rate coefficients (t = 4 ps) approach the corresponding E-FGR CT rate constants obtained with the equilibrium donor (ππ*) state simulations.

FIG. 7.

Time-dependent charge transfer rate coefficients k(t) for ππ* → CT1 (top) and ππ* → CT2 (bottom) transitions, in the bent triad conformation, calculated with NE-FGR and different levels of approximation using models 1–3 (N = 500) as well as IMT via spectral densities [Eqs. (58) and (59)], compared with the IMT result obtained via NEMD (black lines). In model 1, the total number of flips, Nf, are 80 and 180 for ππ* → CT1 and ππ* → CT2 transitions, respectively. The corresponding E-FGR charge transfer rate constants obtained with different levels of approximation using equilibrium ground state parameters (t = 0) are shown on the left of each panel, and those obtained using equilibrium ππ* state parameters (t = ) are shown on the right of each panel. The insets show the comparison of the long-time CT rate coefficients and the E-FGR CT rate constants.

FIG. 7.

Time-dependent charge transfer rate coefficients k(t) for ππ* → CT1 (top) and ππ* → CT2 (bottom) transitions, in the bent triad conformation, calculated with NE-FGR and different levels of approximation using models 1–3 (N = 500) as well as IMT via spectral densities [Eqs. (58) and (59)], compared with the IMT result obtained via NEMD (black lines). In model 1, the total number of flips, Nf, are 80 and 180 for ππ* → CT1 and ππ* → CT2 transitions, respectively. The corresponding E-FGR charge transfer rate constants obtained with different levels of approximation using equilibrium ground state parameters (t = 0) are shown on the left of each panel, and those obtained using equilibrium ππ* state parameters (t = ) are shown on the right of each panel. The insets show the comparison of the long-time CT rate coefficients and the E-FGR CT rate constants.

Close modal
FIG. 8.

Time-dependent charge transfer rate coefficients k(t) for ππ* → CT1 (top) and ππ* → CT2 (bottom) transitions in the linear triad conformation, calculated with NE-FGR and different levels of approximation using models 1–3 (N = 500) as well as IMT via spectral densities [Eqs. (58) and (59)], compared with the IMT result obtained via NEMD (black lines). In model 1, the total number of flips, Nf, is 295 and 245 for ππ* → CT1 and ππ* → CT2 transitions, respectively. Other settings of E-FGR rate constants and the insets are the same as Fig. 7 except for using the parameters of the linear triad.

FIG. 8.

Time-dependent charge transfer rate coefficients k(t) for ππ* → CT1 (top) and ππ* → CT2 (bottom) transitions in the linear triad conformation, calculated with NE-FGR and different levels of approximation using models 1–3 (N = 500) as well as IMT via spectral densities [Eqs. (58) and (59)], compared with the IMT result obtained via NEMD (black lines). In model 1, the total number of flips, Nf, is 295 and 245 for ππ* → CT1 and ππ* → CT2 transitions, respectively. Other settings of E-FGR rate constants and the insets are the same as Fig. 7 except for using the parameters of the linear triad.

Close modal

The temporal profiles of k(t) show a significant nonequilibrium effect, caused by the initial nuclear sampling, in the case of the bent triad (Fig. 7). All three three-state harmonic models capture the nonequilibrium effect remarkably well and are seen to accurately reproduce the results from all-atom NEMD simulations. The ππ* → CT1 process exhibits a significant enhancement in the transient CT rate coefficient, which can be traced back to the fact that the vertical photoexcitation brings the triad to a region that is close to the crossing point between the donor and acceptor PESs (see Fig. 1). In contrast, the ππ* → CT2 process exhibits an initial slowdown in the CT rate coefficient, which is due to the fact that the vertical photoexcitation puts the triad further away from the crossing point, compared to the equilibrium state.66 

For the bent triad, model 1 shows larger fluctuations than models 2 and 3, which can be traced back to the fact that model 1 is based on flipping the sign of the displacements along random modes, whereas models 2 and 3 are based on scaling all the displacements uniformly. However, the fact that model 1 works in the first place suggests that the NE-FGR nonequilibrium CT rate coefficient is largely shaped by the global parameters, such as ErDA,ErDG,ErAG, rather than the displacements along specific modes.

The fact that models 2 and 3 produce the same k(t) profiles confirms our hypothesis that only the nuclear motion along the CT axis determines the photoinduced CT dynamics, while the relaxation along the orthogonal axis (included in model 3) does not influence the CT rate coefficient.

In the insets of models 1–3, the orders of the NE-FGR long-time rate coefficients and E-FGR rate constants calculated at different levels of approximation are similar with the W-AV result corresponding to the fastest CT rate and the Marcus-level result corresponding to the slowest CT rate. This suggests that nuclear quantum delocalization can enhance the CT rate in this system somewhat. The IMT expression based on the three continuous spectral densities JDA(ω), JAG(ω), and JDG(ω) [Eqs. (58) and (59)] is also seen to reproduce the all-atom IMT CT rate coefficient rather well.

In contrast to the bent conformation, the linear triad does not show significant nonequilibrium effects in both ππ* → CT1 and ππ* → CT2 transitions (see Fig. 8). This result is consistent with the all-atom IMT result. All the three-state harmonic models and the IMT via J(ω) predict that the photoinduced CT dynamics follows Marcusian kinetics except for the very small initial relaxation in the CT rate coefficient for the ππ* → CT1 process since the nonvanishing ErDG gives rise to slight DG shifts in our three-state harmonic models, which were obscured by the fluctuations in the all-atom NEMD simulations. In the ππ* → CT1 process, the IMT via spectral densities is seen to overestimate the initial CT rate coefficient although it approaches the equilibrium Marcus rate constant in 10 ps, while models 1–3 show better agreement with the NEMD benchmark result. In the ππ* → CT2 transition, model 1 shows larger fluctuations compared with models 2 and 3, and the fluctuation is comparable to the error bar size in the NEMD result. We conclude that all four approaches generate accurate photoinduced CT dynamics in the linear triad.

The donor state population can be computed from the time-dependent CT rate coefficient using Eq. (5), and we show the result of model 2 for all the CT processes in Fig. 9. The deviation between the NE-FGR predictions and that of the Marcus CT rate constants are significant for the ππ* → CT1 in the bent triad conformation, which displays the largest nonequilibrium effects in the photoinduced CT process, whereas for the rest of the cases, the differences between NE-FGR population dynamics and Marcus theory are small. Figures S1–S3 of the supplementary material show similar behaviors for the donor state population calculated with model 1, model 3, and the IMT via the spectral density approach, respectively.

FIG. 9.

Donor state population decay for ππ* → CT1 (left) and ππ* → CT2 (right) transitions in the bent (top) and linear (bottom) triad conformations, calculated with NE-FGR and different levels of approximation using model 2 (N = 500), compared with the IMT results obtained via NEMD (black lines) and the Marcus rate constant predictions (blue dashed lines).

FIG. 9.

Donor state population decay for ππ* → CT1 (left) and ππ* → CT2 (right) transitions in the bent (top) and linear (bottom) triad conformations, calculated with NE-FGR and different levels of approximation using model 2 (N = 500), compared with the IMT results obtained via NEMD (black lines) and the Marcus rate constant predictions (blue dashed lines).

Close modal

Figure 10 shows the IMT parameters including the time-dependent donor–acceptor energy gap Ut¯ and its variance σt2 as well as the IMT rate coefficients for the bent triad calculated using models 1–3, the IMT expression via spectral densities [Eqs. (58) and (59)], and the benchmark IMT result from NEMD. Model 1 exhibits more fluctuations in the energy gap and therefore in k(t). Note that here the variance of the energy gap is time-independent and exactly the same for those obtained using models 1–3, and it differs from the one obtained using IMT expression via spectral densities by only 1%. The IMT parameters for the linear triad also give rise to a reasonably accurate prediction of the time-dependent CT rate coefficient and is shown in Fig. S4. All the Marcus parameters (e.g., ΓDA, |ΔE|, and Er) as well as the Marcus CT rate constants kM for all cases are provided in Table S2.

FIG. 10.

Comparison of Ut¯ (top), σt2 (middle), and IMT rate coefficient k(t) (bottom) for ππ* → CT1 (left) and ππ* → CT2 (right) transitions in the bent triad conformation calculated using models 1–3 (N = 500), IMT via spectral densities (orange lines), and IMT results obtained via NEMD (black lines).

FIG. 10.

Comparison of Ut¯ (top), σt2 (middle), and IMT rate coefficient k(t) (bottom) for ππ* → CT1 (left) and ππ* → CT2 (right) transitions in the bent triad conformation calculated using models 1–3 (N = 500), IMT via spectral densities (orange lines), and IMT results obtained via NEMD (black lines).

Close modal

The generally excellent agreement of the IMT result with the NE-FGR rate coefficient in all models suggests that the nuclear quantum effects on the CT rate coefficient are small in this triad system, which was also observed in our previous study, where mapping to a spin-boson model was used for calculating E-FGR CT rate constants.72 Moreover, when dealing with large complex systems, whenever IMT is valid, using the IMT level of approximation would be remarkably advantageous since the computational cost is many orders of magnitude lower compared with the other LSC-based approximations of NE-FGR.66 

Finally, we discuss the dependence of the results on the number of nuclear modes in the harmonic model, N. Figure 11 depicts the NE-FGR time-dependent CT rate coefficients for the ππ* → CT1 transition in the bent triad conformation obtained using model 1 (top panels, N = 100, 500, and 1000 with random flipped modes, and N = 500 with evenly flipped modes) and model 2 (bottom panels, N = 100, 200, 500, and 1000 with scaled modes). In general, incorporating more nuclear modes in both models 1 and 2 damps the fluctuations in k(t) and yields a smoother relaxation profile. However, the general trend of the predicted nonequilibrium rate coefficient virtually remains unchanged using different N, and all models with different N seem to be able to capture both the transient behavior and the long-time limit in the nonequilibrium CT rate coefficients remarkably well.

FIG. 11.

The effect of nuclear degrees of freedom N on the time-dependent charge transfer rate coefficients k(t) for ππ* → CT1 transition in the bent triad conformation calculated with NE-FGR and different levels of approximation using model 1 (top panels, N = 100, 500, 1000 with random flips Nf = 16, 80, 160, respectively, as well as N = 500 with evenly flips every six modes, Nf = 80) and model 2 (bottom panels, N = 100, 200, 500, 1000). Other settings of E-FGR rate constants and the insets are the same as Fig. 7.

FIG. 11.

The effect of nuclear degrees of freedom N on the time-dependent charge transfer rate coefficients k(t) for ππ* → CT1 transition in the bent triad conformation calculated with NE-FGR and different levels of approximation using model 1 (top panels, N = 100, 500, 1000 with random flips Nf = 16, 80, 160, respectively, as well as N = 500 with evenly flips every six modes, Nf = 80) and model 2 (bottom panels, N = 100, 200, 500, 1000). Other settings of E-FGR rate constants and the insets are the same as Fig. 7.

Close modal

As discussed above, the fact that model 1 works suggests a collective solvent effect, where the contributions from individual modes are relatively small. For example, comparing the N = 500 results with random flipped modes and the N = 500 results with evenly flipped modes (see the top right panels of Fig. 11), we barely observe any significant difference in the k(t) evaluated on all levels of approximation. In addition, Fig. S5 displays the statistical averages of k(t) as obtained from ten independent parameter sets of model 1 constructed using different sets of {Sj} whose signs were flipped randomly, where the error bars are notably small. It is noted that the number of flipped modes Nf is the same for both stochastic and deterministic approaches since the AG reorganization energy change in flipping any mode is the same (see Table S3), as a result of our spectral density discretization procedure that equally partitions the reorganization energy. The fact that model 1 with different random sets of flipped modes produces very similar k(t) indicates that its stochastic uncertainties are generally negligible, and thus, model 1 is robust in the simulation of nonequilibrium CT dynamics. Similar figures showing the influence of changing nuclear degrees of freedom for other combinations of transitions and conformations are given in Figs. S6–S8 for model 1 and Figs. S9–S12 for model 2.

In this paper, we proposed and demonstrated three protocols for constructing effective three-state harmonic model Hamiltonians and an alternative strategy to compute the IMT rate coefficient via the spectral densities in the context of simulating photoinduced CT rates, where the nonequilibrium effects may be significant due to the initial preparation of the nuclear DOF. Our proposed three-state models (models 1–3) are based on mapping of the all-atom anharmonic Hamiltonian onto a three-state harmonic Hamiltonian. The accuracy and reliability of the three models and the IMT via spectral density approach were assessed by using this strategy to calculate the quantum-mechanical NE-FGR time-dependent rate coefficient and the hierarchy of semiclassical approximations leading to the classical IMT in the CPC60 molecular triad dissolved in explicit THF, which were compared with the benchmark time-dependent CT rate coefficient obtained from all-atom NEMD simulations.

The model construction process starts out with taking inputs from MD simulations in the form of time correlation functions of the energy gaps between different PESs, which are then mapped onto three spectral densities. JDA(ω) is further discretized to N nuclear modes constituting the CT axis parameters shared by models 1–3. For model 1, we employ a one-dimensional CT axis and randomly flip the shifts Sj between the donor and ground state PESs along this CT axis so as to obtain consistently the value of ErAG obtained from MD simulations. Models 2 and 3 are based on a two-dimensional PES scaling approach, where in model 2 we project the nonequilibrium shift on the CT axis and in model 3 we extend the nuclear DOF of model 2 from N to 2N by including the projection of the nonequilibrium shift along the orthogonal axis.

We also showed that IMT can be formulated using only the spectral densities, which can be obtained directly from the MD simulations. Our results show that all three models and the IMT via spectral densities produce remarkably accurate NE-FGR rate coefficients for all the transitions and triad conformations investigated here. The results obtained at different levels of LSC-based approximations are all in agreement with each other, which shows that nuclear quantum effects are minimal for the system under consideration. For the IMT level of approximation, the IMT parameters and the time-dependent rate coefficients obtained via different strategies are seen to coincide with each other. To provide in depth coverage on the characteristic behavior of our models, we also provided an assessment concerning the influence of the number of modes on the model performance and accuracy.

In summary, we have demonstrated the feasibility of three strategies to construct three-state harmonic models that could be used to investigate the significance of the nonequilibrium effects due to initial nuclear preparation and the nuclear quantum dynamical effects in the photoinduced charge transfer dynamics in the condensed phase. As a natural extension of this work, it would be desirable to have a comprehensive study for applying the models proposed here to more complicated multi-level molecular systems and combining with other condensed-phase quantum dynamical methods such as tensor-train split-operator Fourier transform,94 mixed quantum-classical Liouville,95 quasiclassical mapping Hamiltonian,96,97 and generalized quantum master equation,34,35 as well as hierarchical equations of motion39–41 and stochastic equations of motion.98,99 It is also desirable to extend the current treatment of NE-FGR to more than three states, such as incorporating CT1 → CT2 transition. Work toward those goals is currently under way and will be reported in future publications.

See the supplementary material for the derivation of the IMT expression for the three-state harmonic model, the tabulated data of energy corrections for the MD potential energy of different electronic states, the Marcus charge transfer rate constants evaluated using both ground and ππ* electronic-state MD simulations, the plateau values of time-dependent NE-FGR rate coefficients calculated using different models for all transitions and conformations, and other auxiliary information for determining the number of flips for model 1. We also provide the figures for donor state population of all cases evaluated using models 1 and 3, IMT parameters and rate coefficients for linear triad conformations obtained with different models, the stability test of the stochastic effect and the effect of changing nuclear DOF for model 1 in all cases, and the effect of changing nuclear DOF for model 2 in all cases.

X.S. acknowledges support from NYU Shanghai, the National Natural Science Foundation of China (Grant No. 21903054), the Hefei National Laboratory for Physical Sciences at the Microscale (Grant No. KF2020008), and the Program for Eastern Young Scholar at Shanghai Institutions of Higher Learning. E.G. and B.D.D. acknowledge support from the Department of Energy (DOE), Basic Energy Sciences through the Chemical Sciences, Geosciences and Biosciences Division (Grant No. DE-SC0016501). Computing resources were provided by NYU Shanghai HPC.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
J.-L.
Brédas
,
J. E.
Norton
,
J.
Cornil
, and
V.
Coropceanu
, “
Molecular understanding of organic solar cells: The challenges
,”
Acc. Chem. Res.
42
,
1691
1699
(
2009
).
2.
S. M.
Falke
,
C. A.
Rozzi
,
D.
Brida
,
M.
Maiuri
,
M.
Amato
,
E.
Sommer
,
A.
De Sio
,
A.
Rubio
,
G.
Cerullo
,
E.
Molinari
, and
C.
Lienau
, “
Coherent ultrafast charge transfer in an organic photovoltaic blend
,”
Science
344
,
1001
1005
(
2014
).
3.
G.
Bottari
,
G.
de la Torre
,
D. M.
Guldi
, and
T.
Torres
, “
Covalent and noncovalent phthalocyanine–carbon nanostructure systems: Synthesis, photoinduced electron transfer, and application to molecular photovoltaics
,”
Chem. Rev.
110
,
6768
6816
(
2010
).
4.
A. C.
Rizzi
,
M.
van Gastel
,
P. A.
Liddell
,
R. E.
Palacios
,
G. F.
Moore
,
G.
Kodis
,
A. L.
Moore
,
T. A.
Moore
,
D.
Gust
, and
S. E.
Braslavsky
, “
Entropic changes control the charge separation process in triads mimicking photosynthetic charge separation
,”
J. Phys. Chem. A
112
,
4215
4223
(
2008
).
5.
H.
Tian
,
Z.
Yu
,
A.
Hagfeldt
,
L.
Kloo
, and
L.
Sun
, “
Organic redox couples and organic counter electrode for efficient organic dye-sensitized solar cells
,”
J. Am. Chem. Soc.
133
,
9413
9422
(
2011
).
6.
J.
Kurpiers
,
T.
Ferron
,
S.
Roland
,
M.
Jakoby
,
T.
Thiede
,
F.
Jaiser
,
S.
Albrecht
,
S.
Janietz
,
B. A.
Collins
,
I. A.
Howard
, and
D.
Neher
, “
Probing the pathways of free charge generation in organic bulk heterojunction solar cells
,”
Nat. Commun.
9
,
2038
(
2018
).
7.
A.
Mishra
,
M. K. R.
Fischer
, and
P.
Bäuerle
, “
Metal-free organic dyes for dye-sensitized solar cells: From structure: Property relationships to design rules
,”
Angew. Chem., Int. Ed.
48
,
2474
2499
(
2009
).
8.
S. M.
Feldt
,
E. A.
Gibson
,
E.
Gabrielsson
,
L.
Sun
,
G.
Boschloo
, and
A.
Hagfeldt
, “
Design of organic dyes and cobalt polypyridine redox mediators for high-efficiency dye-sensitized solar cells
,”
J. Am. Chem. Soc.
132
,
16714
16724
(
2010
).
9.
Y.
Zhao
and
W.
Liang
, “
Charge transfer in organic molecules for solar cells: Theoretical perspective
,”
Chem. Soc. Rev.
41
,
1075
1087
(
2012
).
10.
P. A.
Liddell
,
G.
Kodis
,
A. L.
Moore
,
T. A.
Moore
, and
D.
Gust
, “
Photonic switching of photoinduced electron transfer in a dithienylethene–porphyrin–fullerene triad molecule
,”
J. Am. Chem. Soc.
124
,
7668
7669
(
2002
).
11.
P. A.
Liddell
,
D.
Kuciauskas
,
J. P.
Sumida
,
B.
Nash
,
D.
Nguyen
,
A. L.
Moore
,
T. A.
Moore
, and
D.
Gust
, “
Photoinduced charge separation and charge recombination to a triplet state in a carotene–porphyrin–fullerene triad
,”
J. Am. Chem. Soc.
119
,
1400
1405
(
1997
).
12.
X.
Sun
,
H.
Wang
, and
W. H.
Miller
, “
Semiclassical theory of electronically nonadiabatic dynamics: Results of a linearized approximation to the initial value representation
,”
J. Chem. Phys.
109
,
7064
(
1998
).
13.
W. H.
Miller
, “
The semiclassical initial value representation: A potentially practical way for adding quantum effects to classical molecular dynamics simulations
,”
J. Phys. Chem. A
105
,
2942
2955
(
2001
).
14.
G.
Tao
and
W. H.
Miller
, “
Time-dependent importance sampling in semiclassical initial value representation calculations for time correlation functions
,”
J. Chem. Phys.
135
,
024104
(
2011
).
15.
X.
Andrade
,
A.
Castro
,
D.
Zueco
,
J. L.
Alonso
,
P.
Echenique
,
F.
Falceto
, and
Á.
Rubio
, “
Modified Ehrenfest formalism for efficient large-scale ab initio molecular dynamics
,”
J. Chem. Theory Comput.
5
,
728
742
(
2009
).
16.
C. P.
van der Vegte
,
A. G.
Dijkstra
,
J.
Knoester
, and
T. l. C.
Jansen
, “
Calculating two-dimensional spectra with the mixed quantum-classical Ehrenfest method
,”
J. Phys. Chem. A
117
,
5970
5980
(
2013
).
17.
A. V.
Akimov
,
R.
Long
, and
O. V.
Prezhdo
, “
Coherence penalty functional: A simple method for adding decoherence in Ehrenfest dynamics
,”
J. Chem. Phys.
140
,
194107
(
2014
).
18.
S.
Fernandez-Alberti
,
D. V.
Makhov
,
S.
Tretiak
, and
D. V.
Shalashilin
, “
Non-adiabatic excited state molecular dynamics of phenylene ethynylene dendrimer using a multiconfigurational Ehrenfest approach
,”
Phys. Chem. Chem. Phys.
18
,
10028
10040
(
2016
).
19.
J. C.
Tully
, “
Molecular dynamics with electronic transitions
,”
J. Chem. Phys.
93
,
1061
(
1990
).
20.
J. E.
Subotnik
,
W.
Ouyang
, and
B. R.
Landry
, “
Can we derive Tully’s surface-hopping algorithm from the semiclassical quantum Liouville equation? Almost, but only with decoherence
,”
J. Chem. Phys.
139
,
214107
(
2013
).
21.
J. E.
Subotnik
,
A.
Jain
,
B.
Landry
,
A.
Petit
,
W.
Ouyang
, and
N.
Bellonzi
, “
Understanding the surface hopping view of electronic transitions and decoherence
,”
Annu. Rev. Phys. Chem.
67
,
387
417
(
2016
).
22.
C. C.
Martens
and
J.-Y.
Fang
, “
Semiclassical-limit molecular dynamics on multiple electronic surfaces
,”
J. Chem. Phys.
106
,
4918
(
1997
).
23.
R.
Kapral
and
G.
Ciccotti
, “
Mixed quantum-classical dynamics
,”
J. Chem. Phys.
110
,
8919
(
1999
).
24.
D.
Mac Kernan
,
G.
Ciccotti
, and
R.
Kapral
, “
Trotter-based simulation of quantum-classical dynamics
,”
J. Phys. Chem. B
112
,
424
432
(
2008
).
25.
A. A.
Kananenka
,
X.
Sun
,
A.
Schubert
,
B. D.
Dunietz
, and
E.
Geva
, “
A comparative study of different methods for calculating electronic transition rates
,”
J. Chem. Phys.
148
,
102304
(
2018
).
26.
S.
Nakajima
, “
On quantum theory of transport phenomena: Steady diffusion
,”
Prog. Theor. Phys.
20
,
948
(
1958
).
27.
R.
Zwanzig
, “
Ensemble method in the theory of irreversibility
,”
J. Chem. Phys.
33
,
1338
1341
(
1960
).
28.
Q.
Shi
and
E.
Geva
, “
A new approach to calculating the memory kernel of the generalized quantum master equation for an arbitrary system–bath coupling
,”
J. Chem. Phys.
119
,
12063
(
2003
).
29.
Q.
Shi
and
E.
Geva
, “
A semiclassical generalized quantum master equation for an arbitrary system-bath coupling
,”
J. Chem. Phys.
120
,
10647
(
2004
).
30.
M.-L.
Zhang
,
B. J.
Ka
, and
E.
Geva
, “
Nonequilibrium quantum dynamics in the condensed phase via the generalized quantum master equation
,”
J. Chem. Phys.
125
,
044106
(
2006
).
31.
S.
Jang
, “
Nonadiabatic quantum Liouville and master equations in the adiabatic basis
,”
J. Chem. Phys.
137
,
22A536
(
2012
).
32.
L.
Kidon
,
E. Y.
Wilner
, and
E.
Rabani
, “
Exact calculation of the time convolutionless master equation generator: Application to the nonequilibrium resonant level model
,”
J. Chem. Phys.
143
,
234110
(
2015
).
33.
A.
Montoya-Castillo
and
D. R.
Reichman
, “
Approximate but accurate quantum dynamics from the mori formalism: I. Nonequilibrium dynamics
,”
J. Chem. Phys.
144
,
184104
(
2016
).
34.
E.
Mulvihill
,
A.
Schubert
,
X.
Sun
,
B. D.
Dunietz
, and
E.
Geva
, “
A modified approach for simulating electronically nonadiabatic dynamics via the generalized quantum master equation
,”
J. Chem. Phys.
150
,
034101
(
2019
).
35.
E.
Mulvihill
,
X.
Gao
,
Y.
Liu
,
A.
Schubert
,
B. D.
Dunietz
, and
E.
Geva
, “
Combining the mapping Hamiltonian linearized semiclassical approach with the generalized quantum master equation to simulate electronically nonadiabatic molecular dynamics
,”
J. Chem. Phys.
151
,
074103
(
2019
).
36.
Y.
Tanimura
, “
Stochastic Liouville, Langevin, Fokker-Planck, and master equation Approaches to quantum dissipative systems
,”
J. Phys. Soc. Jpn.
75
,
082001
(
2006
).
37.
R.-X.
Xu
and
Y.
Yan
, “
Dynamics of quantum dissipation systems interacting with Bosonic canonical bath: Hierarchical equations of motion approach
,”
Phys. Rev. E
75
,
031107
(
2007
).
38.
J.
Jin
,
X.
Zheng
, and
Y.
Yan
, “
Exact dynamics of dissipative electronic systems and quantum transport: Hierarchical equations of motion approach
,”
J. Chem. Phys.
128
,
234703
(
2008
).
39.
K.
Song
and
Q.
Shi
, “
Theoretical study of photoinduced proton coupled electron transfer reaction using the non-perturbative hierarchical equations of motion method
,”
J. Chem. Phys.
146
,
184108
(
2017
).
40.
Y.
Yan
,
T.
Xing
, and
Q.
Shi
, “
A new method to improve the numerical stability of the hierarchical equations of motion for discrete harmonic oscillator modes
,”
J. Chem. Phys.
153
,
204109
(
2020
).
41.
Y.
Tanimura
, “
Numerically ‘exact’ approach to open quantum dynamics: The hierarchical equations of motion (HEOM)
,”
J. Chem. Phys.
153
,
020901
(
2020
).
42.
X.
Sun
and
E.
Geva
, “
Equilibrium Fermi’s golden rule charge transfer rate constants in the condensed phase: The linearized semiclassical method vs classical Marcus theory
,”
J. Phys. Chem. A
120
,
2976
2990
(
2016
).
43.
X.
Sun
and
E.
Geva
, “
Non-Condon equilibrium Fermi’s golden rule electronic transition rate constants via the linearized semiclassical method
,”
J. Chem. Phys.
144
,
244105
(
2016
).
44.
A. J.
Leggett
,
S.
Chakravarty
,
A. T.
Dorsey
,
M. P. A.
Fisher
,
A.
Garg
, and
W.
Zwerger
, “
Dynamics of the dissipative two-state system
,”
Rev. Mod. Phys.
59
,
1
85
(
1987
).
45.
N.
Makri
, “
The linear response approximation and its lowest order corrections: An influence functional approach
,”
J. Phys. Chem. B
103
,
2823
2829
(
1999
).
46.
J.
Cao
and
G. A.
Voth
, “
Modeling physical systems by effective harmonic oscillators: The optimized quadratic approximation
,”
J. Chem. Phys.
102
,
3337
3348
(
1995
).
47.
U.
Weiss
,
Quantum Dissipative Systems
(
World Scientific
,
London
,
2012
).
48.
M. H.
Lee
,
B. D.
Dunietz
, and
E.
Geva
, “
Calculation from first principles of intramolecular golden-rule rate constants for photo-induced electron transfer in molecular donor–acceptor systems
,”
J. Phys. Chem. C
117
,
23391
23401
(
2013
).
49.
M. H.
Lee
,
E.
Geva
, and
B. D.
Dunietz
, “
Calculation from first-principles of golden rule rate constants for photoinduced subphthalocyanine/fullerene interfacial charge transfer and recombination in organic photovoltaic cells
,”
J. Phys. Chem. C
118
,
9780
9789
(
2014
).
50.
R. A.
Marcus
, “
On the theory of oxidation–reduction reactions involving electron transfer. I
,”
J. Chem. Phys.
24
,
966
978
(
1956
).
51.
R. A.
Marcus
, “
Electrostatic free energy and other properties of states having nonequilibrium polarization. I
,”
J. Chem. Phys.
24
,
979
989
(
1956
).
52.
R. A.
Marcus
, “
Electron transfer reactions in chemistry. Theory and experiment
,”
Rev. Mod. Phys.
65
,
599
(
1993
).
53.
P. F.
Barbara
,
T. J.
Meyer
, and
M. A.
Ratner
, “
Contemporary issues in electron transfer research
,”
J. Phys. Chem.
100
,
13148
13168
(
1996
).
54.
Y.
Georgievskii
,
C.-P.
Hsu
, and
R. A.
Marcus
, “
Linear response in theory of electron transfer reactions as an alternative to the molecular harmonic oscillator model
,”
J. Chem.Phys.
110
,
5307
5317
(
1999
).
55.
N. R.
Kestner
,
J.
Logan
, and
J.
Jortner
, “
Thermal electron transfer reactions in polar solvents
,”
J. Phys. Chem.
78
,
2148
2166
(
1974
).
56.
J.
Jortner
and
M.
Bixon
, “
Intramolecular vibrational excitations accompanying solvent-controlled electron transfer reactions
,”
J. Chem. Phys.
88
,
167
170
(
1988
).
57.
D.
Kuciauskas
,
P. A.
Liddell
,
S.
Lin
,
S. G.
Stone
,
A. L.
Moore
,
T. A.
Moore
, and
D.
Gust
, “
Photoinduced electron transfer in carotenoporphyrin–fullerene triads: Temperature and solvent effects
,”
J. Phys. Chem. B
104
,
4307
4321
(
2000
).
58.
T.
Nelson
,
S.
Fernandez-Alberti
,
V.
Chernyak
,
A. E.
Roitberg
, and
S.
Tretiak
, “
Nonadiabatic excited-state molecular dynamics modeling of photoinduced dynamics in conjugated molecules
,”
J. Phys. Chem. B
115
,
5402
5414
(
2011
).
59.
H.
Zhu
,
Y.
Yang
,
K.
Wu
, and
T.
Lian
, “
Charge transfer dynamics from photoexcited semiconductor quantum dots
,”
Annu. Rev. Phys. Chem.
67
,
259
281
(
2016
).
60.
Y.
Song
,
A.
Schubert
,
X.
Liu
,
S.
Bhandari
,
S. R.
Forrest
,
B. D.
Dunietz
,
E.
Geva
, and
J. P.
Ogilvie
, “
Efficient charge generation via hole transfer in dilute organic donor-fullerene blends
,”
J. Phys. Chem. Lett.
11
,
2203
2210
(
2020
).
61.
R. D.
Coalson
,
D. G.
Evans
, and
A.
Nitzan
, “
A nonequilibrium golden rule formula for electronic state populations in nonadiabatically coupled systems
,”
J. Chem. Phys.
101
,
436
(
1994
).
62.
M.
Cho
and
R. J.
Silbey
, “
Nonequilibrium photoinduced electron transfer
,”
J. Chem. Phys.
103
,
595
606
(
1995
).
63.
A. F.
Izmaylov
,
D.
Mendive–Tapia
,
M. J.
Bearpark
,
M. A.
Robb
,
J. C.
Tully
, and
M. J.
Frisch
, “
Nonequilibrium Fermi golden rule for electronic transitions through conical intersections
,”
J. Chem. Phys.
135
,
234106
(
2011
).
64.
X.
Sun
and
E.
Geva
, “
Nonequilibrium Fermi’s golden rule charge transfer rates via the linearized semiclassical method
,”
J. Chem. Theory Comput.
12
,
2926
2941
(
2016
).
65.
X.
Sun
and
E.
Geva
, “
Non-Condon nonequilibrium Fermi’s golden rule rates from the linearized semiclassical method
,”
J. Chem. Phys.
145
,
064109
(
2016
).
66.
Z.
Hu
,
Z.
Tong
,
M. S.
Cheung
,
B. D.
Dunietz
,
E.
Geva
, and
X.
Sun
, “
Photoinduced charge transfer dynamics in the carotenoid-porphyrin-C60 triad via the linearized semiclassical nonequilibrium Fermi’s golden rule
,”
J. Phys. Chem. B
124
,
9579
9591
(
2020
).
67.
F.
Gottwald
,
S. D.
Ivanov
, and
O.
Kühn
, “
Applicability of the Caldeira-Leggett model to vibrational spectroscopy in solution
,”
J. Phys. Chem. Lett.
6
,
2722
2727
(
2015
).
68.
M. K.
Lee
,
P.
Huo
, and
D. F.
Coker
, “
Semiclassical path integral dynamics: Photosynthetic energy transfer with realistic environment interactions
,”
Annu. Rev. Phys. Chem.
67
,
639
668
(
2016
).
69.
M. K.
Lee
and
D. F.
Coker
, “
Modeling electronic-nuclear interactions for excitation energy transfer processes in light-harvesting complexes
,”
J. Phys. Chem. Lett.
7
,
3171
3178
(
2016
).
70.
X.
Sun
and
E.
Geva
, “
Exact vs Asymptotic spectral densities in the Garg-Onuchic-Ambegaokar charge transfer model and its effect on Fermi’s golden rule rate constants
,”
J. Chem. Phys.
144
,
044106
(
2016
).
71.
P. L.
Walters
,
T. C.
Allen
, and
N.
Makri
, “
Direct determination of discrete harmonic bath parameters from molecular dynamics simulations
,”
J. Comput. Chem.
38
,
110
115
(
2017
).
72.
Z.
Tong
,
X.
Gao
,
M. S.
Cheung
,
B. D.
Dunietz
,
E.
Geva
, and
X.
Sun
, “
Charge transfer rate constants for the carotenoid-porphyrin-C60 molecular triad dissolved in tetrahydrofuran: The spin-boson model vs the linearized semiclassical approximation
,”
J. Chem. Phys.
153
,
044105
(
2020
).
73.
D.
Carbonera
,
M.
Di Valentin
,
C.
Corvaja
,
G.
Agostini
,
G.
Giacometti
,
P. A.
Liddell
,
D.
Kuciauskas
,
A. L.
Moore
,
T. A.
Moore
, and
D.
Gust
, “
EPR investigation of photoinduced radical pair formation and decay to a triple state in a carotene–porphyrin–fullerene triad
,”
J. Am. Chem. Soc.
120
,
4398
4405
(
1998
).
74.
D.
Gust
,
T. A.
Moore
, and
A. L.
Moore
, “
Mimicking photosynthetic solar energy transduction
,”
Acc. Chem. Res.
34
,
40
48
(
2001
).
75.
N.
Spallanzani
,
C. A.
Rozzi
,
D.
Varsano
,
T.
Baruah
,
M. R.
Pederson
,
F.
Manghi
, and
A.
Rubio
, “
Photoexcitation of a light-harvesting supramolecular triad: A time-dependent DFT study
,”
J. Phys. Chem. B
113
,
5345
5349
(
2009
).
76.
A. C.
Rozzi
,
S.
Maria Falke
,
N.
Spallanzani
,
A.
Rubio
,
E.
Molinari
,
D.
Brida
,
M.
Maiuri
,
G.
Cerullo
,
H.
Schramm
,
J.
Christoffers
, and
C.
Lienau
, “
Quantum coherence controls the charge separation in a prototypical artificial light–harvesting system
,”
Nat. Commun.
4
,
1602
1607
(
2013
).
77.
G.
Su
,
A.
Czader
,
D.
Homouz
,
G.
Bernardes
,
S.
Mateen
, and
M. S.
Cheung
, “
Multiscale simulation on a light-harvesting molecular triad
,”
J. Phys. Chem. B
116
,
8460
8473
(
2012
).
78.
D.
Balamurugan
,
A. J. A.
Aquino
,
F.
de Dios
,
L.
Flores
, Jr.
,
H.
Lischka
, and
M. S.
Cheung
, “
Multiscale simulation of the ground and photo-induced charge-separated states of a molecular triad in polar organic solvent: Exploring the conformations, fluctuations, and free energy landscapes
,”
J. Phys. Chem. B
117
,
12065
12075
(
2013
).
79.
A. K.
Manna
,
D.
Balamurugan
,
M. S.
Cheung
, and
B. D.
Dunietz
, “
Unraveling the mechanism of photoinduced charge transfer in carotenoid–porphyrin–C60 molecular triad
,”
J. Phys. Chem. Lett.
6
,
1231
1237
(
2015
).
80.
O. N.
Starovoytov
,
P.
Zhang
,
P.
Cieplak
, and
M. S.
Cheung
, “
Induced polarization restricts the conformational distribution of a light–harvesting molecular triad in the ground state
,”
Phys. Chem. Chem. Phys.
19
,
22969
22980
(
2017
).
81.
X.
Sun
,
P.
Zhang
,
Y.
Lai
,
K. L.
Williams
,
M. S.
Cheung
,
B. D.
Dunietz
, and
E.
Geva
, “
Computational study of charge-transfer dynamics in the carotenoid–porphyrin–C60 molecular triad solvated in explicit tetrahydrofuran and its spectroscopic signature
,”
J. Phys. Chem. C
122
,
11288
11299
(
2018
).
82.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
, “
VMD: Visual molecular dynamics
,”
J. Mol. Graphics
14
,
33
38
(
1996
).
83.
S. S.
Khohlova
,
V. A.
Mikhailova
, and
A. I.
Ivanov
, “
Three-centered model of ultrafast photoinduced charge transfer: Continuum dielectric approach
,”
J. Chem. Phys.
124
,
114507
(
2006
).
84.
A. I.
Ivanov
and
V. G.
Tkachev
, “
Exact solution of three-level model of excited state electron transfer symmetry breaking in quadrupolar molecules
,”
J. Chem. Phys.
151
,
124309
(
2019
).
85.
D.
Brian
and
X.
Sun
, “
Linear-response and nonlinear-response formulations of the instantaneous Marcus theory for nonequilibrium photoinduced charge transfer
,”
J. Chem. Theory Comput.
17
,
2065
2079
(
2021
).
86.
R.
Baer
and
D.
Neuhauser
, “
Density functional theory with correct long-range asymptotic behavior
,”
Phys. Rev. Lett.
94
,
043002
(
2005
).
87.
E.
Livshits
and
R.
Baer
, “
A well-tempered density functional theory of electrons in molecules
,”
Phys. Chem. Chem. Phys.
9
,
2932
2941
(
2007
).
88.
T.
Stein
,
H.
Eisenberg
,
L.
Kronik
, and
R.
Baer
, “
Fundamental gaps in finite systems from eigenvalues of a generalized Kohn-Sham method
,”
Phys. Rev. Lett.
105
,
266802
(
2010
).
89.
Y.
Shao
,
Z.
Gan
,
E.
Epifanovsky
,
A. T.
Gilbert
,
M.
Wormit
,
J.
Kussmann
,
A. W.
Lange
,
A.
Behn
,
J.
Deng
,
X.
Feng
 et al, “
Advances in molecular quantum chemistry contained in the Q-chem 4 program package
,”
Mol. Phys.
113
,
184
215
(
2015
).
90.
A. A.
Voityuk
and
N.
Rösch
, “
Fragment charge difference method for estimating donor–acceptor electronic coupling: Application to DNA π-stacks
,”
J. Chem. Phys.
117
,
5607
5616
(
2002
).
91.
D. A.
Case
,
I. Y.
Ben-Shalom
,
S. R.
Brozell
,
D. S.
Cerutti
,
T. E.
Cheatham
 III
,
V. W. D.
Cruzeiro
,
T. A.
Darden
,
R. E.
Duke
,
D.
Ghoreishi
,
M. K.
Gilson
,
H.
Gohlke
,
A. W.
Goetz
 et al,
AMBER 2018
(
University of California
,
San Francisco
,
2018
).
92.
G.
Ciccotti
and
J. P.
Ryckaert
, “
Molecular dynamics simulation of rigid molecules
,”
Comput. Phys. Rep.
4
,
346
392
(
1986
).
93.
T.
Darden
,
D.
York
, and
L.
Pedersen
, “
Particle mesh Ewald: An N-log(N) method for Ewald sums in large systems
,”
J. Chem. Phys.
98
,
10089
10092
(
1993
).
94.
S. M.
Greene
and
V. S.
Batista
, “
Tensor-train split-operator fourier transform (TT-SOFT) method: Multidimensional nonadiabatic quantum dynamics
,”
J. Chem. Theory Comput.
13
,
4034
4042
(
2017
).
95.
A. A.
Kananenka
,
C.-Y.
Hsieh
,
J.
Cao
, and
E.
Geva
, “
Accurate long-time mixed quantum-classical Liouville dynamics via the transfer tensor method
,”
J. Phys. Chem. Lett.
7
,
4809
4814
(
2016
).
96.
X.
Gao
,
M. A. C.
Saller
,
Y.
Liu
,
A.
Kelly
,
J. O.
Richardson
, and
E.
Geva
, “
Benchmarking quasiclassical mapping Hamiltonian methods for simulating electronically nonadiabatic molecular dynamics
,”
J. Chem. Theory Comput.
16
,
2883
2895
(
2020
).
97.
Y.
Liu
,
X.
Gao
,
Y.
Lai
,
E.
Mulvihill
, and
E.
Geva
, “
Electronic dynamics through conical intersections via quasiclassical mapping Hamiltonian methods
,”
J. Chem. Theory Comput.
16
,
4479
4488
(
2020
).
98.
Y.-A.
Yan
and
J.
Shao
, “
Stochastic description of quantum Brownian dynamics
,”
Front. Phys.
11
,
110309
(
2016
).
99.
Y.
Ke
and
Y.
Zhao
, “
An extension of stochastic hierarchy equations of motion for the equilibrium correlation functions
,”
J. Chem. Phys.
146
,
214105
(
2017
).

Supplementary Material