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–C_{60} (CPC_{60}) 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 CPC_{60}/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 CPC_{60}/THF can be described accurately by the effective harmonic three-state models and that nuclear quantum effects are small in this system.

## I. INTRODUCTION

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 theory^{50–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 Hamiltonian^{72} has been recently found to lead to remarkably accurate equilibrium FGR (E-FGR) CT rate constants in the case of the carotenoid–porphyrin–C_{60} (CPC_{60}) molecular triad^{11,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 CPC_{60} triad conformations in explicit THF (see Fig. 1). The triad is initially equilibrated on the ground (G) state, CPC_{60}, before it is photoexcited impulsively to the P-localized *ππ*^{*} state, CP^{*}C_{60}. Following the impulsive photoexcitation at time 0, there could be a nonradiative transition to the excited P-to-C_{60} CT state, $CP+C60\u2212$, which is denoted as CT1,

or to the excited C-to-C_{60} charge separated state, $C+PC60\u2212$, which is denoted as CT2,

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-C_{60} [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 CPC_{60} 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.

## II. THEORY

### A. Nonequilibrium Fermi’s golden rule (NE-FGR)

Consider a charge transfer system with the Hamiltonian

Here, |*D*⟩ and |*A*⟩ represent the diabatic donor and acceptor electronic states, respectively, and $\u0124D/A$ are the corresponding nuclear Hamiltonians, $\u0124D/A=P^2/2+VD/A(R^)$, where $R^=R^1,\u2026,R^N$ and $P^=P^1,\u2026,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 $\rho ^G=e\u2212\beta \u0124G/Trn[e\u2212\beta \u0124G]$. Here, $\u0124G=P^2/2+VG(R^)$ is the ground state nuclear Hamiltonian, $VG(R^)$ is the ground state PES, *β* = 1/*k*_{B}*T* is the inverse temperature, Tr_{n}(·) denotes the trace over the nuclear Hilbert space, and Tr_{n} Tr_{e} denotes the trace over both the nuclear and electronic Hilbert spaces. The donor state population, *P*_{D}(*t*), at a later time *t* is given by

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

where the *time-dependent rate coefficient* is defined as

Here,

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, $\rho ^D=e\u2212\beta \u0124D/Trn[e\u2212\beta \u0124D]$, *C*(*t*, *τ*) reduces to

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

Under those conditions, the donor state population could be estimated by the exponential decay *P*_{D}(*t*) = exp(−*k*_{D→A}*t*). It should be noted that the rate constant *k*_{D→A} 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 *k*_{D→A} are significantly different and the timescale for reaching thermal equilibrium on the donor PES is comparable to or longer than the timescale of CT, $\u223ckD\u2192A\u22121$.

### B. The hierarchy of LSC-based approximations

The linearized semiclassical (LSC) approximation provides a computationally feasible route for calculating NE-FGR rates. The main idea behind LSC is to express $Ct,\tau $ 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}

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

and “C” denotes the classical nuclear sampling using the corresponding classical distribution *ρ*_{G,Cl}(**R**, **P**). Following the initial sampling of (**R**_{0}, **P**_{0}) 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 **R**_{t}, and then backward in time from time *t* to *t* − *τ* (*τ* = 0 → *t*), integrating the donor–acceptor energy gap *U*(**R**) = *U*_{DA}(**R**) = *V*_{D}(**R**) − *V*_{A}(**R**) in the phase factor for each (*t*, *τ*). “AV” denotes that the dynamics over *τ* is propagated on the average PES, *V*_{av}(**R**) = (*V*_{D}(**R**) + *V*_{A}(**R**))/2, leading to trajectories of $R\tau av$, “*D*” denotes that the *τ*-dynamics is on the donor surface, leading to trajectories of $R\tau D$, and “0” denotes no *τ*-dynamics beyond the *t*-relaxation on the donor PES until **R**_{t}. 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}

where we assume that the instantaneous distribution of the energy gap *U*_{t} is Gaussian with a time-dependent mean $Ut\xaf$ and the corresponding variance $\sigma t2=Ut2\xaf\u2212Ut\xaf2$ at time *t* after the photoexcitation so that $Prob(Ut)=1/2\pi \sigma t2exp\u2212(Ut\u2212Ut\xaf)2/(2\sigma t2).$ Here, $Ut\xaf$ and $\sigma t2$ can be calculated using energies obtained from all-atom molecular dynamics simulations.^{66} The IMT CT rate coefficient *k*^{IMT}(*t*) can also be expressed equivalently in terms of time-dependent reaction free energy Δ*E*(*t*) and reorganization energy *E*_{r}(*t*),

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 by^{42}

The corresponding reorganization energy and reaction free energy satisfy the relation $Er=\sigma eq2/(2kBT)=\u2212\Delta E\u2212\u27e8U\u27e9$ and the activation energy $Ea=kBT\u27e8U\u27e92/(2\sigma eq2)$, where ⟨·⟩ is the equilibrium average on the donor state PES and the equilibrium variance of the energy gap $\sigma eq2=\u27e8U2\u27e9\u2212\u27e8U\u27e92$.

## III. CONSTRUCTING THREE-STATE HARMONIC MODEL HAMILTONIANS

### A. Mapping an anharmonic all-atom Hamiltonian onto the spin-boson model

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 by^{47}

Here, $\sigma ^x=|D\u2009\u3009\u3008A|+|A\u2009\u3009\u3008D|$ and $\sigma ^z=|D\u2009\u3009\u3008D|\u2212|A\u2009\u3009\u3008A|$ are the Pauli operators; Γ_{DA} is the electronic coupling coefficient; Δ*E* = −*ℏω*_{DA} is the donor-to-acceptor reaction free energy; ${R^j,P^j,\omega j|j=1,\u2026,N}$ are the mass-weighted coordinates, momenta, and frequencies associated with the nuclear normal modes, respectively; and {*c*_{j}} are the coupling coefficients between the electronic DOF and nuclear normal modes. The frequencies and the coupling coefficients, {*ω*_{j}, *c*_{j}}, are often given in terms of spectral density, which is defined by

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, *C*_{UU}(*t*), from an equilibrium MD simulation on the donor PES,^{72}

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):

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)=\sigma eq2$,

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 *j*th mode (*j* = 1, 2, …, *N*):^{71}

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 *j*th mode is determined, the corresponding coupling coefficients, *c*_{j}, are obtained via the relation

### B. Three-state harmonic models for photoinduced charge transfer

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:

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

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 *j*th mode (*j* = 1, 2, …, *N*) is given by [see Eqs. (25) and (19) for the definition of *c*_{j}]

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

Expressed in terms of ${Rjeq}$ instead of {*c*_{j}}, the spectral density can be written in the following form:

Finally, {*S*_{j}} 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}

The time-dependent IMT parameters $Ut\xaf$ and $\sigma 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)

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

It should be noted that $\sigma t2=\sigma 02$ for the harmonic model, which implies that any significant time-dependence of $\sigma 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).

### C. Mapping an anharmonic all-atom Hamiltonian onto a three-state harmonic model

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,

Here, *U*_{XY} = *V*_{X} − *V*_{Y} (*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*),

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 {*S*_{j}}. The equilibrium geometry displacements between the donor and acceptor states, ${Rjeq}$, can be unambiguously determined from $ErDA=Er$ as in Eq. (28),

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

One can also show that

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

It should be noted that only the absolute values of {*S*_{j}} 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=\u2211j12\omega 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., *U*_{DG} = *U*_{DA} + *U*_{AG}), so forcing the nonequilibrium shifts to be colinear with *U*_{DA}, 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 *S*_{j} for a given $Rjeq$. The relative signs of *S*_{j} and $Rjeq$ can clearly lead to either a positive or a negative contribution of the *j*th mode to $ErAG$ by the amount of $\Delta ErAG(j)=\xb12\omega j2RjeqRjDG$ (see Table S3 of the supplementary material). Since $ErAG=\u2211j12\omega j2(Rjeq+Sj)2$ can be satisfied for different selections of relative signs of *S*_{j} and $Rjeq$, satisfying Eq. (49) does not dictate a unique choice of the relative signs of ${Rjeq,Sj}$.

In practice, we choose $Rjeq=RjDA\u22650\u2009(j=1,\u2026,N)$ and randomly flip the signs of {*S*_{j}} such that the deviations between $ErAG,eff=\u2211j12\omega j2(Rjeq+Sj)2$ and $ErAG=CUUAG(0)/(2kBT)$ are minimal. Alternatively, one can flip the signs of {*S*_{j}} relative to ${Rjeq}$ evenly, namely, every few modes, until $ErAG,eff\u2248ErAG$ [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 {*S*_{j}} 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., $RjXY\u221dErXY$, 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

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

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

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 {*S*_{j}} in Eq. (27) are ${RjDG}$ projected onto the CT axis $SjCT$, i.e.,

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 2*N* 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:

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\xaf$ 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:

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:

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

## IV. SIMULATION TECHNIQUES

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 CPC_{60} 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) functional^{86–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

Here, *E*_{α}(**r**^{Triad}) 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., **r**^{Triad}. 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\alpha MD(R)$ and the excitation energy *E*_{α}(**r**^{Triad}) of the triad would double count the energy of the triad intramolecular interactions, so we need to subtract this double-counted contribution *V*_{α,T}(**r**^{Triad}) evaluated with the force fields where only the triad is present. Finally, we arrive at the energy correction *W*_{α} = *E*_{α}(**r**^{Triad}) − *V*_{α,T}(**r**^{Triad}) 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$ Å $\xd7\u2009100$ Å $\xd7\u2009100$ Å periodic cubic box, performed using the AMBER 18 package.^{91} The van der Waals cutoff radius was set to be 9 Å. The SHAKE algorithm^{92} 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 $kcal\u2009mol\u22121$Å^{−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 × 10^{7} 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 × 10^{7} 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, *N*_{f}, such that the deviation of the effective reorganization energy $ErAG,eff$ from the accurate value $ErAG$ is minimized,

## V. RESULTS AND DISCUSSION

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 CPC_{60} 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.

### A. All-atom simulation of CPC_{60} in THF and its harmonic models

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,\pi \pi *(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,\pi \pi *(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/\pi \pi *,CT2(t)$ (orange and blue lines). As expected, the amplitude of $CUUG/\pi \pi *,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.

Transition . | Conf. . | $ErDA$ . | $ErDG$ . | $ErAG$ . |
---|---|---|---|---|

ππ^{*} → 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 |

Transition . | Conf. . | $ErDA$ . | $ErDG$ . | $ErAG$ . |
---|---|---|---|---|

ππ^{*} → 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 |

Next, Figs. 5 and 6 show the three spectral densities *J*_{DA}(*ω*), *J*_{AG}(*ω*), and *J*_{DG}(*ω*), 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)=\sigma eq2=2kBTErXY$. This justifies using the frequencies obtained from *J*_{DA}(*ω*) to construct the three-state models (shown as cyan vertical lines in Fig. 5). However, in the linear triad conformation, the spectral densities *J*_{DA}(*ω*) and *J*_{AG}(*ω*) have the same profile for the same transition except for a scaling factor, while the spectral density *J*_{DG}(*ω*) 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 *J*_{DA}(*ω*) 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 *J*_{DG}(*ω*) that would in turn give rise to rather small equilibrium geometry shifts between the donor and the ground PESs in the three-state models.

### B. Photoinduced CT dynamics of CPC_{60} in THF via three-state harmonic models

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.

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 *J*_{DA}(*ω*), *J*_{AG}(*ω*), and *J*_{DG}(*ω*) [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.

Figure 10 shows the IMT parameters including the time-dependent donor–acceptor energy gap $Ut\xaf$ and its variance $\sigma 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 *E*_{r}) as well as the Marcus CT rate constants *k*^{M} for all cases are provided in Table S2.

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.

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 {*S*_{j}} whose signs were flipped randomly, where the error bars are notably small. It is noted that the number of flipped modes *N*_{f} 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.

## VI. CONCLUDING REMARKS

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 CPC_{60} 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. *J*_{DA}(*ω*) 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 *S*_{j} 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 2*N* 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 motion^{39–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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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

## REFERENCES

_{60}triad via the linearized semiclassical nonequilibrium Fermi’s golden rule

_{60}molecular triad dissolved in tetrahydrofuran: The spin-boson model vs the linearized semiclassical approximation

_{60}molecular triad

_{60}molecular triad solvated in explicit tetrahydrofuran and its spectroscopic signature