Charge transfer rate constants were calculated for the carotenoid-porphyrin-C_{60} (CPC_{60}) molecular triad dissolved in explicit tetrahydrofuran. The calculation was based on mapping the all-atom anharmonic Hamiltonian of this system onto the spin-boson Hamiltonian. The mapping was based on discretizing the spectral density from the time correlation function of the donor–acceptor potential energy gap, as obtained from all-atom molecular dynamics simulations. Different spin-boson Hamiltonians were constructed for each of the possible transitions between the three excited electronic states in two different triad conformations. The rate constants of three possible transitions were calculated via the quantum-mechanically exact Fermi’s golden rule (FGR), as well as a progression of more approximate expressions that lead to the classical Marcus expression. The advantage of the spin-boson approach is that once the mapping is established, the quantum-mechanically exact FGR and the hierarchy of approximations are known in closed form. The classical Marcus charge transfer rate constants obtained with the spin-boson Hamiltonians were found to reproduce those obtained from all-atom simulations with the linearized semiclassical approximation, thereby confirming the equivalence of the two approaches for this system. Within the spin-boson Hamiltonian, we also found that the quantum-mechanically exact FGR rate constants were significantly enhanced compared to the classical Marcus theory rate constants for two out of three transitions in one of the two conformations under consideration. The results confirm that mapping to the spin-boson model can yield accurate predictions for charge transfer rate constants in a system as complex as CPC_{60} dissolved in tetrahydrofuran.

## I. INTRODUCTION

Photoinduced charge transfer (CT) lies at the heart of processes fundamental to natural and synthetic systems that convert light energy into other forms of energy (electrical, chemical, mechanical, etc.).^{1–12} Quantitative modeling of the underlying CT dynamics based on explicitly molecular models is necessary for understanding the operational principles of such systems, as well as developing rational design principles toward enhancing their performance.^{13–20}

Marcus theory has served as the workhorse for estimating CT rate constants in a wide range of condensed-phase systems.^{21–24} Within this approach, CT rate constants can be expressed in terms of only three parameters: the donor–acceptor electronic coupling coefficient, the reorganization energy, and the reaction free energy. As such, Marcus theory offers an intuitive approach toward understanding CT rates, as well as fitting and interpreting experimentally measured CT rate constants. However, the validation of the assumptions underlying Marcus theory requires being able to compare the Marcus theory rate constants to rate constants obtained via a less restrictive approach.^{14,17,19,25–28}

The linearized semiclassical (LSC) approximation for the equilibrium Fermi’s golden rule (E-FGR) CT rate constant^{29–33} represents such a practical and less restrictive approach for calculating CT rate constants.^{14,17} Importantly, the LSC expression for the E-FGR CT rate constant can be applied to complex molecular systems described by general anharmonic force fields.^{19} The LSC E-FGR expression for the CT rate constant can also be related to the classical Marcus expression for the rate constant.^{14} More specifically, Marcus theory is seen to be based on three major approximations: (1) treating the nuclear degrees of freedom (DOF) as classical, (2) assuming that the electronic dephasing time [as measured by the lifetime of *C*(*t*) in Eq. (4)] is shorter than the time scale of nuclear motion, and (3) assuming that the equilibrium distribution of the donor–acceptor potential energy gap is Gaussian.^{14} However, those assumptions can break down in the deep inverted and normal regimes and/or at sufficiently low temperatures.^{13,14,17,24,34–36}

An alternative and commonly used approach for going beyond Marcus theory is based on mapping the all-atom anharmonic Hamiltonian onto an effective spin-boson Hamiltonian.^{25–28,37–42} The mapping onto the spin-boson Hamiltonian is often based on the equilibrium donor–acceptor potential energy gap time correlation function as obtained from classical molecular dynamics (MD) simulations performed on an all-atom anharmonic molecular model when the system is in the donor state.^{37}

The main advantage of this effective spin-boson Hamiltonian approach is that once the mapping is established, the quantum mechanically exact E-FGR rate constant, as well as the hierarchy of approximations that lead from it to the corresponding classical Marcus limit, are known in closed form.^{14,16} Interestingly, it has been shown that the LSC-based E-FGR expression coincides with the quantum-mechanically exact E-FGR rate constant when the system can be described by the spin-boson Hamiltonian.^{14} However, the validity of mapping a highly anharmonic all-atom Hamiltonian onto a spin-boson Hamiltonian is rarely tested.

In this paper, we pursue a practical strategy for testing the self-consistency of mapping onto the spin-boson Hamiltonian by comparison to the rate constants obtained based on the all-atom anharmonic Hamiltonian via the LSC approximation. We demonstrate this approach in the context of the carotenoid-porphyrin-C_{60} (CPC_{60}) molecular triad^{11,43–47} dissolved in explicit tetrahydrofuran (THF) solvent^{19,48–51} (see Fig. 1), for which CT rate constants were recently calculated based on an all-atom anharmonic Hamiltonian via the LSC approach.^{19}

The remainder of this paper is organized as follows: The theory underlying the E-FGR CT rate constant and its LSC approximation is outlined in Sec. II. The mapping onto the spin-boson Hamiltonian is outlined in Sec. III. The all-atom anharmonic model for CPC_{60} dissolved in liquid THF and MD simulation techniques are described in Sec. IV. The results are reported in Sec. V. Conclusions and outlook are provided in Sec. VI.

## II. THE EQUILIBRIUM FERMI’S GOLDEN RULE (E-FGR) CT RATE CONSTANT AND ITS LSC APPROXIMATION

### A. The E-FGR CT rate constant

Consider a two-state system with the Hamiltonian *Ĥ* = *Ĥ*_{0} + *Ĥ*_{1}, where the zeroth-order Hamiltonian, *Ĥ*_{0}, and the perturbation, *Ĥ*_{I}, are given by

Here, |*D*⟩ and |*A*⟩ are the diabatic donor and acceptor states, respectively, $\u0124D/A=P^2/2+VD/A(R^)$ are the corresponding nuclear Hamiltonians, $R^=(R^1,R^2,\u2026,R^N)$ and $P^=(P^1,P^2,\u2026,P^N)$ are the mass-weighted nuclear coordinates and momenta, $VD/A(R^)$ are the donor/acceptor potential energy surfaces (PESs), and Γ is the donor–acceptor diabatic electronic coupling coefficient within the Condon approximation (in what follows, boldfaced quantities correspond to vectors and quantities capped with $\u02c6$ correspond to quantum-mechanical operators, while the same quantities without the $\u02c6$ correspond to their classical counterparts).

It should be noted that the model system may have more than two states. In such a case, the donor and acceptor states correspond to the two electronic states associated with a given transition. Thus, different choices of donor and acceptor states will give rise to different transition rate constants (see below).

In the case of the E-FGR, we assume that the system starts out at the donor electronic state with the nuclear DOF at thermal equilibrium,

Here, Tr_{N}[·] is the trace over the nuclear DOF, *k*_{B} is the Boltzmann constant, and *T* is the absolute temperature.

The E-FGR donor-to-acceptor CT rate constant is given by^{52}

where

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

The E-FGR expression in Eq. (4) is fully quantum-mechanical. The LSC approximation for the correlation function *C*(*t*) is given by *C*^{W-AV}(*t*) in Eq. (5). A hierarchy of increasingly more approximate and easier to calculate expressions that eventually lead to the classical Marcus expression, *C*^{M}(*t*), are given by Eqs. (6)–(10),^{14}

Here, *U*(**R**) = *V*_{D}(**R**) − *V*_{A}(**R**) is the donor–acceptor potential energy gap; “W” denotes the initial sampling of the nuclear DOF based on the corresponding Wigner distribution, $\rho D,Weq(R0,P0)$; “C” denotes the initial sampling of the nuclear DOF based on the corresponding classical Boltzmann distribution, $\rho D,Cleq(R0,P0)$; $R\tau av$ is obtained via classical dynamics of the nuclear DOF on the average PES, *V*_{av} = (*V*_{D}(**R**) + *V*_{A}(**R**))/2; $R\tau D$ is obtained via classical dynamics of the nuclear DOF on the donor PES, *V*_{D}(**R**); and “0” denotes the approximation that nuclear motion can be neglected over the lifetime of *C*(*t*) so that **R**_{τ} = **R**_{0}.

The Marcus limit, *C*^{M}(*t*), corresponds to a second-order cumulant expansion of *C*^{C-0}(*t*) and depends on the average donor–acceptor potential energy gap, ⟨*U*⟩, and the corresponding standard deviation *σ*_{U}, defined by

The Gaussian form of *C*^{M}(*t*) also makes it possible to perform the integral in Eq. (3) explicitly, to yield the following expression for the Marcus CT rate constant in terms of Γ^{2}, ⟨*U*⟩, and *σ*_{U}:

Comparing Eq. (13) with the more familiar expression for the classical Marcus CT rate constant when given in terms of the reorganization energy, *E*_{r}, and reaction free energy, Δ*E*,

then yields the following relations between {*E*_{r}, Δ*E*}, the corresponding activation energy, $EA=\Delta E+Er2/4Er$, and {⟨*U*⟩, *σ*_{U}}:^{14}

It should be noted that ⟨*U*⟩ < 0 corresponds to the normal Marcus regime (|Δ*E*| < *E*_{r}), while ⟨*U*⟩ > 0 corresponds to the inverted Marcus regime (|Δ*E*| > *E*_{r}).

Unlike *E*_{r} and Δ*E*, obtaining ⟨*U*⟩ and *σ*_{U} from equilibrium MD simulations based on an all-atom anharmonic Hamiltonian is relatively straightforward. It should be noted that although truncating the cumulant expansion at second order would be exact when the donor and acceptor PESs are harmonic and identical, except for a shift in equilibrium energy and geometry, no explicit mapping to a harmonic Hamiltonian is necessary for calculating the Marcus rate constant in Eq. (13) [although the relations in Eq. (15) can be shown to be self-consistent in such a case].

## III. MAPPING ONTO A SPIN-BOSON HAMILTONIAN

One of the most popular approaches for simulating quantum dynamics in the condensed phase is based on mapping onto the spin-boson model.^{28} The spin-boson Hamiltonian is given by

Here, $\sigma ^x=|D\u2009\u2009A|+|A\u2009\u2009D|$ and $\sigma ^z=|D\u2009\u2009D|\u2212|A\u2009\u2009A|$ are the Pauli operators; Γ is the electronic coupling coefficient; *ℏω*_{DA} = −Δ*E*, where Δ*E* is the donor-to-acceptor reaction free energy; ${x^j,p^j,\omega j|j=1,\u2026,N}$ are the mass-weighted nuclear coordinates, momenta, and frequencies, respectively; and {*c*_{j}} are the coupling coefficients between the electronic DOF and nuclear DOF.

The frequencies and the coupling coefficients {*ω*_{j}, *c*_{j}} for a spin-boson model are often expressed in terms of a spectral density function, which is defined by

A popular procedure for mapping an all-atom anharmonic Hamiltonian like that of the molecular triad in liquid THF solution described in Sec. IV, onto the spin-boson Hamiltonian, is based on the average donor–acceptor potential energy gap, ⟨*U*⟩, and time correlation function, *C*_{UU}(*t*),

Here, *U*(*t*) is the instantaneous donor–acceptor potential energy gap at time *t* and the average ⟨⋯⟩ is based on classical thermal equilibrium when the system is in the donor state. Note that *C*_{UU}(*t*) is different from *C*(*t*). The relationship between the spectral density, *J*(*ω*), and *C*_{UU}(*t*) is given by^{37}

To obtain a finite set of nuclear modes from the continuous spectral density function in Eq. (20), we follow the procedure in Ref. 37. More specifically, given the time correlation function *C*_{UU}(*t*), as obtained from classical MD simulations performed on the all-atom anharmonic system, and assuming that the overall number of nuclear modes is *N* (*N* ≫ 1), the frequency of the *j*th nuclear mode in the spin-boson Hamiltonian is dictated by the following equation:^{37}

The coupling coefficients, {*c*_{j}}, are obtained via the relation

where the reorganization energy is *E*_{r} = *C*_{UU}(0)/2*k*_{B}*T* and the reaction free energy is thus Δ*E* = −*ℏω*_{DA} = −⟨*U*⟩ − *C*_{UU}(0)/2*k*_{B}*T* [see Eq. (15) and note that $CUU(0)\u2261\sigma U2$].

Here, $Rjeq=2cj/\omega j2$ is the displacement between the donor and acceptor PESs along the *j*th mode coordinate. Importantly, for the spin-boson model, the approximations in Eqs. (5)–(10) can be obtained in closed form,^{14}

## IV. AN ALL-ATOM ANHARMONIC HAMILTONIAN FOR CPC_{60} DISSOLVED IN THF

The all-atom anharmonic Hamiltonian of CPC_{60} dissolved in THF was adopted from Refs. 19 and 50. The model takes into account four electronic states of the triad: (1) the ground state, CPC_{60}; (2) the P-localized excitonic *ππ*^{*} state, $CP*C60$; (3) the excited P-to-C_{60} CT state, $CP+C60\u2212$, which is referred to as CT1; and (4) the excited C-to-C_{60} charge-separated state, $C+PC60\u2212$, which is referred to as CT2.

Following photoexcitation from the ground state to the *ππ*^{*} state, the triad is believed to undergo a nonradiative transition from *ππ*^{*} to CT1, followed by another nonradiative transition from CT1 to CT2,^{50}

The potential energy surfaces (PESs) for the entire system of the triad and solvent differ from one electronic state to another by the excitation energies of the triad and the partial charges of the triad atoms. Thus, the atom–atom Coulomb pair interaction terms in the force fields differ from one electronic state to another. We also consider two characteristic triad conformations: *bent* and *linear* (see Fig. 1). Each conformation has its own set of force field parameters for the four electronic states under consideration.

The excitation energies, partial charges, and electronic coupling coefficients for this system were obtained from time-dependent density functional theory (TDDFT) using the range-separated hybrid BNL functional.^{53–55} The electronic coupling coefficients between the electronically excited states were assumed to be constant within the Condon approximation and calculated via the fragment charge difference (FCD) method.^{56} For details of force field parameters, we refer the reader to Refs. 19, 49, and 51.

MD simulations for a system containing one triad molecule and 6741 THF molecules in a 100 Å × 100 Å × 100 Å periodic cubic box at 300 K were performed within the AMBER 16 package.^{57} Particle mesh Ewald summation was used to calculate the electrostatic interactions.^{58} The van der Waals cutoff radius was set to 12 Å. The SHAKE algorithm^{59} was used to constrain the covalent bonds involving hydrogen atoms. The time step was set to 1.0 fs. The system was equilibrated at temperature 300 K using a Langevin thermostat with 1.0 ps^{−1} collision frequency. 200 trajectories of 100 ps length each were generated on the *ππ*^{*} and CT1 PESs. ⟨*U*⟩, *σ*_{U}, and *C*_{UU}(*t*) were calculated by averaging over those trajectories. In the case of the linear triad, a harmonic potential of 100 kcal mol^{−1} Å^{−2} force constant was applied to constrain the end-to-end distance at 54.4 Å. No such constraint was used in the case of the bent triad.

It should be noted that based on our previous study, CT in this system is induced by the solute–solvent interactions rather than by the triad’s intramolecular vibrations.^{19} In addition, the CT1 → CT2 process was previously found to be significantly more efficient in the linear conformation.^{19,50}

It should also be noted that different donor-to-acceptor transitions will give rise to different spin-boson Hamiltonians (see Fig. 2). Thus, the triad in liquid solution model under consideration in this paper needs to be mapped onto six different spin-boson Hamiltonians—one for each of the three transitions (*ππ*^{*} → CT1, *ππ*^{*} → CT2, and CT1 → CT2) and the two triad conformations (bent and linear).

## V. RESULTS

The characteristic behavior of the donor–acceptor potential energy gap as a function of time is shown in Fig. 3(a) for different transitions in the bent and linear triad conformations. The average, ⟨*U*⟩, and range, as measured by *σ*_{U}, are seen to be strongly dependent on the electronic state and conformation. Figure 3(b) shows that the energy gap distribution is Gaussian.

The values of ⟨*U*⟩ for the different transitions in the bent conformation are significantly different but similar in the linear conformation. This behavior is consistent with Fig. 2, which shows that the *ππ*^{*} → CT1, *ππ*^{*} → CT2, and CT1 → CT2 transitions correspond to a positive, a slightly negative, and a strongly negative energy gap when the system is at the equilibrium of the parent state in the bent conformation [remember that *U*(**R**) = *V*_{D}(**R**) − *V*_{A}(**R**)]. At the same time, all three transitions correspond to slightly negative energy gaps when the system is at the equilibrium of the parent state in the linear conformation.

Since $Er=\sigma U2/2kBT$, the values of energy gap fluctuation *σ*_{U} obtained from MD simulation (see Fig. 3) are used to calculate *E*_{r}, whose trends in different CT transitions are shown in Fig. 2. More specifically, *E*_{r}, and thereby *σ*_{U}, for the *ππ*^{*} → CT1 transition is seen to be smaller than for the *ππ*^{*} → CT2 and CT1 → CT2 transitions. This can be traced back to the fact that the CT2 state involves the carotenoid as the donor. More specifically, the fact that carotenoid is spatially more extended than both *C*_{60} and porphyrin, and thereby exposed to a larger number of solvent molecules, implies that transitions to the CT2 state involve a larger reorganization energy.

The donor–acceptor potential energy gap correlation functions, *C*_{UU}(*t*), for the three transitions and two conformations are shown in Fig. 4. The initial values, $CUU(0)=\sigma U2=2kBTEr$, are seen to be consistent with the abovementioned trends in *σ*_{U} and *E*_{r} (see the top panel in Fig. 4). The lifetime of *C*_{UU}(*t*) is ∼1.0 ps (much longer than the decay time of the E-FGR electronic dephasing correlation function *C*(*t*) that is ∼20 fs) and seen to be relatively insensitive to the transition and conformation (see the bottom panel in Fig. 4). This time scale can be traced back to orientational motion of the THF molecules that dominate the reorganization in this system and is expected to be insensitive to the transition and triad conformation.

The spectral densities for the different transitions and conformations as obtained from the correlation functions shown in Fig. 4 are shown in Fig. 5. Since the decay profile of the normalized *C*_{UU}(*t*)/*C*_{UU}(0) is relatively insensitive to the transition and conformation, the spectral densities for different transitions and conformations hence have a similar shape and only differ by a scaling factor determined by the relative values of the corresponding $CUU(0)=\sigma U2=2kBTEr$. The major difference between the different spin-boson models therefore lies in the reorganization energy *E*_{r} and the average energy gap ⟨*U*⟩. Consequently, the coupling coefficients {*c*_{j}} differ by a scaling factor [see Eq. (22)], and the reaction free energy, Δ*E*, is determined by Δ*E* = −⟨*U*⟩ − *E*_{r}.

Also shown in Fig. 4 are the frequencies of the 200 nuclear modes as obtained based on Eq. (21). As expected, the vast majority of the modes have frequencies below ∼300 cm^{−1} and can be traced back to the orientational motion of the solvent molecules.

The hierarchy of CT rate constants [see Eqs. (5)–(10)] for the three transitions and two conformations, as obtained based on the corresponding spin-boson models, are shown in Fig. 6. The CT rate constants, *k*, are shown as a function of |Δ*E*| on a semilog plot, with the arrows indicating the values of |Δ*E*| that correspond to the actual transition at each conformation. The corresponding rate constants are shown in Table I.

Transition . | Conf. . | k^{W-AV/exact}
. | k^{W-0}
. | k^{C-AV}
. | k^{C-D}
. | k^{C-0}
. | k^{M}
. |
---|---|---|---|---|---|---|---|

ππ^{*} → CT1 | Bent | 1.8 ± 0.1 × 10^{12} | 2.0 ± 0.2 × 10^{12} | 9.7 ± 0.8 × 10^{11} | 3.1 ± 0.5 × 10^{11} | 1.6 ± 0.5 × 10^{11} | 1.2 ± 0.2 × 10^{11} |

Linear | 9.1 ± 0.1 × 10^{11} | 7.9 ± 0.1 × 10^{11} | 1.4 ± 0.1 × 10^{12} | 9.0 ± 0.2 × 10^{11} | 1.1 ± 0.1 × 10^{11} | 1.1 ± 0.1 × 10^{12} | |

ππ^{*} → CT2 | Bent | 1.2 ± 0.1 × 10^{7} | 1.0 ± 0.1 × 10^{7} | 2.1 ± 0.1 × 10^{7} | 6.3 ± 0.9 × 10^{6} | 7 ± 1 × 10^{6} | 7 ± 1 × 10^{6} |

Linear | 2.3 ± 0.1 × 10^{8} | 2.0 ± 0.1 × 10^{8} | 4.1 ± 0.1 × 10^{8} | 1.6 ± 0.5 × 10^{8} | 1.9 ± 0.1 × 10^{8} | 1.9 ± 0.1 × 10^{8} | |

CT1 → CT2 | Bent | 1.4 ± 0.3 × 10^{5} | 1.4 ± 0.2 × 10^{6} | 1.2 ± 0.4 × 10^{7} | 1.2 ± 0.2 × 10^{4} | 5.2 ± 0.8 × 10^{2} | 6 ± 1 × 10^{2} |

Linear | 3.4 ± 0.5 × 10^{9} | 3.0 ± 0.3 × 10^{9} | 1 ± 2 × 10^{9} | 6 ± 2 × 10^{8} | 6 ± 2 × 10^{8} | 6 ± 2 × 10^{8} |

Transition . | Conf. . | k^{W-AV/exact}
. | k^{W-0}
. | k^{C-AV}
. | k^{C-D}
. | k^{C-0}
. | k^{M}
. |
---|---|---|---|---|---|---|---|

ππ^{*} → CT1 | Bent | 1.8 ± 0.1 × 10^{12} | 2.0 ± 0.2 × 10^{12} | 9.7 ± 0.8 × 10^{11} | 3.1 ± 0.5 × 10^{11} | 1.6 ± 0.5 × 10^{11} | 1.2 ± 0.2 × 10^{11} |

Linear | 9.1 ± 0.1 × 10^{11} | 7.9 ± 0.1 × 10^{11} | 1.4 ± 0.1 × 10^{12} | 9.0 ± 0.2 × 10^{11} | 1.1 ± 0.1 × 10^{11} | 1.1 ± 0.1 × 10^{12} | |

ππ^{*} → CT2 | Bent | 1.2 ± 0.1 × 10^{7} | 1.0 ± 0.1 × 10^{7} | 2.1 ± 0.1 × 10^{7} | 6.3 ± 0.9 × 10^{6} | 7 ± 1 × 10^{6} | 7 ± 1 × 10^{6} |

Linear | 2.3 ± 0.1 × 10^{8} | 2.0 ± 0.1 × 10^{8} | 4.1 ± 0.1 × 10^{8} | 1.6 ± 0.5 × 10^{8} | 1.9 ± 0.1 × 10^{8} | 1.9 ± 0.1 × 10^{8} | |

CT1 → CT2 | Bent | 1.4 ± 0.3 × 10^{5} | 1.4 ± 0.2 × 10^{6} | 1.2 ± 0.4 × 10^{7} | 1.2 ± 0.2 × 10^{4} | 5.2 ± 0.8 × 10^{2} | 6 ± 1 × 10^{2} |

Linear | 3.4 ± 0.5 × 10^{9} | 3.0 ± 0.3 × 10^{9} | 1 ± 2 × 10^{9} | 6 ± 2 × 10^{8} | 6 ± 2 × 10^{8} | 6 ± 2 × 10^{8} |

One case where the deviation between the quantum mechanically exact FGR CT rate constant and the classical Marcus theory CT rate constant is significant (one order of magnitude) is the *ππ*^{*} → CT1 transition in the bent conformation. This transition is also the only one that falls in the inverted regime (|Δ*E*| > *E*_{r}), where such deviations are expected and are attributed to tunneling.

Another case where the deviation between the quantum mechanically exact FGR CT rate constant and the classical Marcus theory CT rate constant is significant (three orders of magnitude) is the CT1 → CT2 transition in the bent conformation. Interestingly, this transition falls in the deep normal regime (|Δ*E*| ≪ *E*_{r}). Hence, the enhancement is likely due to the contribution of tunneling becoming increasingly more important relative to that of the competing classical activated process with increasing activation energy.

Far more modest enhancement of the E-FGR CT rate constant relative to the classical Marcus theory CT rate constant is observed in the remaining four cases, which all fall in the normal regime and close vicinity of the Marcus turnover (|Δ*E*| ∼ *E*_{r}). These results confirm the validity of classical Marcus theory for describing the CT kinetics in this system when close to the Marcus turnover region.

It should also be noted that the abovementioned significant quantum enhancement of *ππ*^{*} → CT1 and CT1 → CT2 in the bent conformation does not alter the overall CT kinetics as obtained from classical Marcus theory. This is because the rate determining step is CT1 → CT2, which is still much faster in the linear conformation, even after the abovementioned quantum enhancement of the CT1 → CT2 rate in the bent conformation.

Finally, in Table II, we compare the classical Marcus rate constants obtained based on mapping to a spin-boson Hamiltonian via Eq. (30) (*k*^{M}) to those obtained based on the all-atom anharmonic Hamiltonian via Eq. (10) ($kMDM$). The excellent agreement between *k*^{M} and $kMDM$ for all three transitions in both conformations suggests that the two approaches are essentially equivalent for the molecular system under consideration.^{19}

Transition . | Conf. . | Γ . | ⟨U⟩
. | σ_{U}
. | ΔE . | E_{r}
. | E_{A}
. | $kMDM$ . | k^{M}
. |
---|---|---|---|---|---|---|---|---|---|

ππ^{*} → CT1 | Bent | 2.4 × 10^{−2} | 0.500 | 0.163 | −1.014 | 0.513 | 0.122 | 1.6 ± 0.1 × 10^{11} | 1.2 ± 0.2 × 10^{11} |

Linear | 9.0 × 10^{−3} | −0.163 | 0.181 | −0.472 | 0.636 | 0.011 | 1.1 ± 0.1 × 10^{12} | 1.1 ± 0.1 × 10^{12} | |

ππ^{*} → CT2 | Bent | 4.5 × 10^{−5} | −0.445 | 0.273 | −0.981 | 1.426 | 0.035 | 7.7 ± 0.3 × 10^{6} | 7 ± 1 × 10^{6} |

Linear | 2.0 × 10^{−4} | −0.413 | 0.296 | −1.280 | 1.693 | 0.025 | 1.9 ± 0.1 × 10^{8} | 1.9 ± 0.1 × 10^{8} | |

CT1 → CT2 | Bent | 8.6 × 10^{−5} | −1.333 | 0.271 | −0.086 | 1.418 | 0.313 | 6 ± 1 × 10^{2} | 6 ± 1 × 10^{2} |

Linear | 1.0 × 10^{−3} | −0.685 | 0.267 | −0.691 | 1.376 | 0.085 | 6 ± 2 × 10^{8} | 6 ± 2 × 10^{8} |

Transition . | Conf. . | Γ . | ⟨U⟩
. | σ_{U}
. | ΔE . | E_{r}
. | E_{A}
. | $kMDM$ . | k^{M}
. |
---|---|---|---|---|---|---|---|---|---|

ππ^{*} → CT1 | Bent | 2.4 × 10^{−2} | 0.500 | 0.163 | −1.014 | 0.513 | 0.122 | 1.6 ± 0.1 × 10^{11} | 1.2 ± 0.2 × 10^{11} |

Linear | 9.0 × 10^{−3} | −0.163 | 0.181 | −0.472 | 0.636 | 0.011 | 1.1 ± 0.1 × 10^{12} | 1.1 ± 0.1 × 10^{12} | |

ππ^{*} → CT2 | Bent | 4.5 × 10^{−5} | −0.445 | 0.273 | −0.981 | 1.426 | 0.035 | 7.7 ± 0.3 × 10^{6} | 7 ± 1 × 10^{6} |

Linear | 2.0 × 10^{−4} | −0.413 | 0.296 | −1.280 | 1.693 | 0.025 | 1.9 ± 0.1 × 10^{8} | 1.9 ± 0.1 × 10^{8} | |

CT1 → CT2 | Bent | 8.6 × 10^{−5} | −1.333 | 0.271 | −0.086 | 1.418 | 0.313 | 6 ± 1 × 10^{2} | 6 ± 1 × 10^{2} |

Linear | 1.0 × 10^{−3} | −0.685 | 0.267 | −0.691 | 1.376 | 0.085 | 6 ± 2 × 10^{8} | 6 ± 2 × 10^{8} |

## VI. CONCLUDING REMARKS

In this paper, we tested the self-consistency of calculating CT rate constants based on the spin-boson Hamiltonian by comparison to CT rate constants obtained based on the all-atom anharmonic Hamiltonian via the LSC approximation. We did so in the context of a system consisting of a CPC_{60} molecular triad dissolved in liquid THF.

The mapping of the all-atom anharmonic Hamiltonian onto the spin-boson Hamiltonian was based on obtaining the spectral densities from the correlation functions of the donor–acceptor potential energy gap. Six spectral densities were obtained for the molecular system under consideration (one for each choice of transition and conformation). However, it was observed that the six different spectral densities only differ by a scaling factor, which can be related to the corresponding reorganization energy. The advantage of mapping to the spin-boson Hamiltonian is that once the mapping is established, the quantum-mechanically exact FGR and the hierarchy of approximations are known in closed form.

For the solvated triad system under consideration, we found that the classical Marcus CT rate constants obtained via mapping to a spin-boson Hamiltonian reproduce the classical Marcus CT rate constants obtained via the LSC approximation rather accurately. We also found that two out of the six CT rate constants were significantly enhanced by accounting for quantum effects within the E-FGR framework, either in the inverted regime or deep normal regime.

The strategy outlined here is general and can be used in many other systems where CT, and more generally electronic transitions, occurs. The approach can also be extended beyond E-FGR to account for nonequilibrium and strong electronic coupling effects. Work on such extensions is currently underway in our groups and will be reported in future publications.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

## ACKNOWLEDGMENTS

X.S. acknowledges support from NYU Shanghai, the National Natural Science Foundation of China (Grant No. 21903054), the Shanghai Sailing Program (Grant No. 19YF1435600), and the Program for Eastern Young Scholar at Shanghai Institutions of Higher Learning. E.G., B.D.D., and M.S.C. 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 the NYU Shanghai High Performance Computing Center.