Charge transfer rate constants were calculated for the carotenoid-porphyrin-C60 (CPC60) 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 CPC60 dissolved in tetrahydrofuran.

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 constant29–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-C60 (CPC60) molecular triad11,43–47 dissolved in explicit tetrahydrofuran (THF) solvent19,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 

FIG. 1.

The two main conformations of the carotenoid-porphyrin-C60 (CPC60) molecular triad dissolved in explicit tetrahydrofuran (THF): (a) the bent triad conformation and (b) the linear triad conformation.

FIG. 1.

The two main conformations of the carotenoid-porphyrin-C60 (CPC60) molecular triad dissolved in explicit tetrahydrofuran (THF): (a) the bent triad conformation and (b) the linear triad conformation.

Close modal

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 CPC60 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.

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

Ĥ0=ĤD|DD|+ĤA|AA|,ĤI=Γ(|DA|+|AD|).
(1)

Here, |D⟩ and |A⟩ are the diabatic donor and acceptor states, respectively, ĤD/A=P^2/2+VD/A(R^) are the corresponding nuclear Hamiltonians, R^=(R^1,R^2,,R^N) and P^=(P^1,P^2,,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 ˆ correspond to quantum-mechanical operators, while the same quantities without the ˆ 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,

ρ^(0)=ρ^Deq|DD|,ρ^Deq=eĤD/kBT/TrNeĤD/kBT.
(2)

Here, TrN[·] is the trace over the nuclear DOF, kB is the Boltzmann constant, and T is the absolute temperature.

The E-FGR donor-to-acceptor CT rate constant is given by52 

kDA=22Re0dtC(t),
(3)

where

C(t)=Γ2TrNρ^DeqeiĤAt/eiĤDt/.
(4)

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

CW-AV(t)=Γ2dR0dP0ρD,Weq(R0,P0)expi0tdτU(Rτav),
(5)
CW-0(t)=Γ2dR0dP0ρD,Weq(R0,P0)expiU(R0)t,
(6)
CC-AV(t)=Γ2dR0dP0ρD,Cleq(R0,P0)expi0tdτU(Rτav),
(7)
CC-D(t)=Γ2dR0dP0ρD,Cleq(R0,P0)expi0tdτU(RτD),
(8)
CC-0(t)=Γ2dR0dP0ρD,Cleq(R0,P0)expiU(R0)t,
(9)
CM(t)=Γ2expiUt122σU2t2.
(10)

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

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

U=dRdPρD,Cleq(R,P)U(R),
(11)
σU2=U2U2.
(12)

The Gaussian form of CM(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:

kDAM=Γ22πσU2expU22σU2.
(13)

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

kDAM=Γ2πkBTErexpΔE+Er24kBTEr,
(14)

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

Er=σU2/(2kBT)=ΔEU,
(15)
EA=kBT2U2σU2.
(16)

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

Unlike Er 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].

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

Ĥ=Γσ^x+ωDA2σ^z+j=1Np^j22+12ωj2x^j2cjx^jσ^z.
(17)

Here, σ^x=|DA|+|AD| and σ^z=|DD||AA| 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,ωj|j=1,,N} are the mass-weighted nuclear coordinates, momenta, and frequencies, respectively; and {cj} are the coupling coefficients between the electronic DOF and nuclear DOF.

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

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

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, CUU(t),

CUU(t)=U(0)U(t)U2.
(19)

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 CUU(t) is different from C(t). The relationship between the spectral density, J(ω), and CUU(t) is given by37 

J(ω)=βω0dtCUU(t)cos(ωt).
(20)

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 CUU(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 jth nuclear mode in the spin-boson Hamiltonian is dictated by the following equation:37 

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

The coupling coefficients, {cj}, are obtained via the relation

cj=Er/(2N)ωj  (j=1,2,,N),
(22)

where the reorganization energy is Er = CUU(0)/2kBT and the reaction free energy is thus ΔE = −ℏωDA = −⟨U⟩ − CUU(0)/2kBT [see Eq. (15) and note that CUU(0)σU2].

The spin-boson Hamiltonian in Eq. (17) can be written in the form of Eq. (1),

Ĥ=ĤD|DD|+ĤA|AA|+Γ|DA|+|AD|,
(23)
ĤD=j=1Np^j22+12j=1Nωj2x^j2+ωDA,
(24)
ĤA=j=1Np^j22+12j=1Nωj2x^jRjeq2.
(25)

Here, Rjeq=2cj/ωj2 is the displacement between the donor and acceptor PESs along the jth mode coordinate. Importantly, for the spin-boson model, the approximations in Eqs. (5)–(10) can be obtained in closed form,14 

CW-AV(t)=C(t)=Γ2expiωDAtj=1Nωj(Rjeq)22×cothβωj21cos(ωjt)+isin(ωjt),
(26)
CW-0(t)=Γ2expiωDAtj=1Nωj(Rjeq)22×cothβωj2ωj2t22+iωjt,
(27)
CC-AV(t)=Γ2expiωDAtj=1Nωj(Rjeq)22×2βωj1cos(ωjt)+isin(ωjt),
(28)
CC-D(t)=Γ2expiωDAtj=1Nωj(Rjeq)22×2βωj1cos(ωjt)+iωjt,
(29)
CC-0(t)=CM(t)=Γ2expiωDAtj=1Nωj(Rjeq)22ωjt2β+iωjt.
(30)

It should be noted that in the case of the spin-boson model, CW-AV(t) coincides with C(t) {where C(t) is the quantum mechanically exact correlation function [see Eq. (4)]} and CC-0(t) coincides with CM(t).14 

The all-atom anharmonic Hamiltonian of CPC60 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, CPC60; (2) the P-localized excitonic ππ* state, CP*C60; (3) the excited P-to-C60 CT state, CP+C60, which is referred to as CT1; and (4) the excited C-to-C60 charge-separated state, C+PC60, 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

CPC60hνCP*C60(ππ*)CP+C60(CT1)C+PC60(CT2).
(31)

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 algorithm59 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 CUU(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).

FIG. 2.

Potential energy surfaces for different donor-to-acceptor transitions between excited states in the bent (solid) and linear (dashed) triad conformations: (a) ππ* → CT1, (b) ππ* → CT2, and (c) CT1 → CT2. Er is the reorganization energy. The relative positions of the parabolas are quantitatively obtained using the Marcus parameters, Er and ΔE, from MD simulations, and the curvatures and the positions of the donor-state parabolas are chosen to be the same for simplicity.

FIG. 2.

Potential energy surfaces for different donor-to-acceptor transitions between excited states in the bent (solid) and linear (dashed) triad conformations: (a) ππ* → CT1, (b) ππ* → CT2, and (c) CT1 → CT2. Er is the reorganization energy. The relative positions of the parabolas are quantitatively obtained using the Marcus parameters, Er and ΔE, from MD simulations, and the curvatures and the positions of the donor-state parabolas are chosen to be the same for simplicity.

Close modal

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.

FIG. 3.

Donor–acceptor potential energy gap fluctuations (a) and distributions (b) for transitions ππ*CT1 (black), ππ*CT2 (blue), and CT1 → CT2 (red) in the bent (top) and linear (bottom) triad conformations, over a time interval of 100 ps. The average and the standard deviation values for each case are marked on the right of panel (a).

FIG. 3.

Donor–acceptor potential energy gap fluctuations (a) and distributions (b) for transitions ππ*CT1 (black), ππ*CT2 (blue), and CT1 → CT2 (red) in the bent (top) and linear (bottom) triad conformations, over a time interval of 100 ps. The average and the standard deviation values for each case are marked on the right of panel (a).

Close modal

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) = VD(R) − VA(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=σU2/2kBT, the values of energy gap fluctuation σU obtained from MD simulation (see Fig. 3) are used to calculate Er, whose trends in different CT transitions are shown in Fig. 2. More specifically, Er, 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 C60 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, CUU(t), for the three transitions and two conformations are shown in Fig. 4. The initial values, CUU(0)=σU2=2kBTEr, are seen to be consistent with the abovementioned trends in σU and Er (see the top panel in Fig. 4). The lifetime of CUU(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.

FIG. 4.

Un-normalized (top) and normalized (bottom) donor–acceptor potential energy gap correlation functions, CUU(t), for transitions ππ*CT1 (black), ππ*CT2 (blue), and CT1 → CT2 (red) in the bent (solid) and linear (dashed) triad conformations, as obtained from molecular dynamics simulations. Each CUU(t) was calculated by averaging over 4 × 106 configurations.

FIG. 4.

Un-normalized (top) and normalized (bottom) donor–acceptor potential energy gap correlation functions, CUU(t), for transitions ππ*CT1 (black), ππ*CT2 (blue), and CT1 → CT2 (red) in the bent (solid) and linear (dashed) triad conformations, as obtained from molecular dynamics simulations. Each CUU(t) was calculated by averaging over 4 × 106 configurations.

Close modal

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 CUU(t)/CUU(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)=σU2=2kBTEr. The major difference between the different spin-boson models therefore lies in the reorganization energy Er and the average energy gap ⟨U⟩. Consequently, the coupling coefficients {cj} differ by a scaling factor [see Eq. (22)], and the reaction free energy, ΔE, is determined by ΔE = −⟨U⟩ − Er.

FIG. 5.

Spectral density J(ω) calculated from the energy gap correlation function CUU(t) via Eq. (20) and the discrete harmonic mode frequencies (blue lines) for bent (left) and linear (right) triad conformations of transitions ππ*CT1 (top), ππ*CT2 (middle), and CT1 → CT1 (bottom).

FIG. 5.

Spectral density J(ω) calculated from the energy gap correlation function CUU(t) via Eq. (20) and the discrete harmonic mode frequencies (blue lines) for bent (left) and linear (right) triad conformations of transitions ππ*CT1 (top), ππ*CT2 (middle), and CT1 → CT1 (bottom).

Close modal

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.

FIG. 6.

A semilog plot (log base 10) of CT rate constants at different levels of approximation for transitions ππ*CT1 (top), ππ*CT2 (middle), and CT1 → CT1 (bottom) as a function of |ΔE| for the bent (left) and linear (right) triad conformations, as obtained based on the spin-boson Hamiltonian. The arrows indicate the actual value of |ΔE| for the triad in THF.

FIG. 6.

A semilog plot (log base 10) of CT rate constants at different levels of approximation for transitions ππ*CT1 (top), ππ*CT2 (middle), and CT1 → CT1 (bottom) as a function of |ΔE| for the bent (left) and linear (right) triad conformations, as obtained based on the spin-boson Hamiltonian. The arrows indicate the actual value of |ΔE| for the triad in THF.

Close modal
TABLE I.

CT rate constants (in unit of s−1) calculated at different levels of approximation for the spin-boson models that correspond to different donor-to-acceptor transitions in the bent and linear triad conformations.

TransitionConf.kW-AV/exactkW-0kC-AVkC-DkC-0kM
ππ* → CT1 Bent 1.8 ± 0.1 × 1012 2.0 ± 0.2 × 1012 9.7 ± 0.8 × 1011 3.1 ± 0.5 × 1011 1.6 ± 0.5 × 1011 1.2 ± 0.2 × 1011 
 Linear 9.1 ± 0.1 × 1011 7.9 ± 0.1 × 1011 1.4 ± 0.1 × 1012 9.0 ± 0.2 × 1011 1.1 ± 0.1 × 1011 1.1 ± 0.1 × 1012 
ππ* → CT2 Bent 1.2 ± 0.1 × 107 1.0 ± 0.1 × 107 2.1 ± 0.1 × 107 6.3 ± 0.9 × 106 7 ± 1 × 106 7 ± 1 × 106 
 Linear 2.3 ± 0.1 × 108 2.0 ± 0.1 × 108 4.1 ± 0.1 × 108 1.6 ± 0.5 × 108 1.9 ± 0.1 × 108 1.9 ± 0.1 × 108 
CT1 → CT2 Bent 1.4 ± 0.3 × 105 1.4 ± 0.2 × 106 1.2 ± 0.4 × 107 1.2 ± 0.2 × 104 5.2 ± 0.8 × 102 6 ± 1 × 102 
 Linear 3.4 ± 0.5 × 109 3.0 ± 0.3 × 109 1 ± 2 × 109 6 ± 2 × 108 6 ± 2 × 108 6 ± 2 × 108 
TransitionConf.kW-AV/exactkW-0kC-AVkC-DkC-0kM
ππ* → CT1 Bent 1.8 ± 0.1 × 1012 2.0 ± 0.2 × 1012 9.7 ± 0.8 × 1011 3.1 ± 0.5 × 1011 1.6 ± 0.5 × 1011 1.2 ± 0.2 × 1011 
 Linear 9.1 ± 0.1 × 1011 7.9 ± 0.1 × 1011 1.4 ± 0.1 × 1012 9.0 ± 0.2 × 1011 1.1 ± 0.1 × 1011 1.1 ± 0.1 × 1012 
ππ* → CT2 Bent 1.2 ± 0.1 × 107 1.0 ± 0.1 × 107 2.1 ± 0.1 × 107 6.3 ± 0.9 × 106 7 ± 1 × 106 7 ± 1 × 106 
 Linear 2.3 ± 0.1 × 108 2.0 ± 0.1 × 108 4.1 ± 0.1 × 108 1.6 ± 0.5 × 108 1.9 ± 0.1 × 108 1.9 ± 0.1 × 108 
CT1 → CT2 Bent 1.4 ± 0.3 × 105 1.4 ± 0.2 × 106 1.2 ± 0.4 × 107 1.2 ± 0.2 × 104 5.2 ± 0.8 × 102 6 ± 1 × 102 
 Linear 3.4 ± 0.5 × 109 3.0 ± 0.3 × 109 1 ± 2 × 109 6 ± 2 × 108 6 ± 2 × 108 6 ± 2 × 108 

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| > Er), 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| ≪ Er). 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| ∼ Er). 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) (kM) to those obtained based on the all-atom anharmonic Hamiltonian via Eq. (10) (kMDM). The excellent agreement between kM and kMDM for all three transitions in both conformations suggests that the two approaches are essentially equivalent for the molecular system under consideration.19 

TABLE II.

Marcus CT rate constants (in unit of s−1) and parameters for different donor-to-acceptor transitions in the bent and linear triad conformations, as obtained from molecular dynamics simulations via the LSC approximation (kMDM) and spin-boson model (kM). Here, electronic coupling coefficients Γ, average value of the donor–acceptor potential energy gap ⟨U⟩ and its standard deviation σU, reorganization energy Er, reaction free energy ΔE, and activation energy EA are given in unit of eV.

TransitionConf.ΓUσUΔEErEAkMDMkM
ππ* → CT1 Bent 2.4 × 10−2 0.500 0.163 −1.014 0.513 0.122 1.6 ± 0.1 × 1011 1.2 ± 0.2 × 1011 
 Linear 9.0 × 10−3 −0.163 0.181 −0.472 0.636 0.011 1.1 ± 0.1 × 1012 1.1 ± 0.1 × 1012 
ππ* → CT2 Bent 4.5 × 10−5 −0.445 0.273 −0.981 1.426 0.035 7.7 ± 0.3 × 106 7 ± 1 × 106 
 Linear 2.0 × 10−4 −0.413 0.296 −1.280 1.693 0.025 1.9 ± 0.1 × 108 1.9 ± 0.1 × 108 
CT1 → CT2 Bent 8.6 × 10−5 −1.333 0.271 −0.086 1.418 0.313 6 ± 1 × 102 6 ± 1 × 102 
 Linear 1.0 × 10−3 −0.685 0.267 −0.691 1.376 0.085 6 ± 2 × 108 6 ± 2 × 108 
TransitionConf.ΓUσUΔEErEAkMDMkM
ππ* → CT1 Bent 2.4 × 10−2 0.500 0.163 −1.014 0.513 0.122 1.6 ± 0.1 × 1011 1.2 ± 0.2 × 1011 
 Linear 9.0 × 10−3 −0.163 0.181 −0.472 0.636 0.011 1.1 ± 0.1 × 1012 1.1 ± 0.1 × 1012 
ππ* → CT2 Bent 4.5 × 10−5 −0.445 0.273 −0.981 1.426 0.035 7.7 ± 0.3 × 106 7 ± 1 × 106 
 Linear 2.0 × 10−4 −0.413 0.296 −1.280 1.693 0.025 1.9 ± 0.1 × 108 1.9 ± 0.1 × 108 
CT1 → CT2 Bent 8.6 × 10−5 −1.333 0.271 −0.086 1.418 0.313 6 ± 1 × 102 6 ± 1 × 102 
 Linear 1.0 × 10−3 −0.685 0.267 −0.691 1.376 0.085 6 ± 2 × 108 6 ± 2 × 108 

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 CPC60 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.

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

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.

1.
J.-L.
Brédas
,
J. E.
Norton
,
J.
Cornil
, and
V.
Coropceanu
,
Acc. Chem. Res.
42
,
1691
(
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
,
Science
344
,
1001
(
2014
).
3.
G.
Bottari
,
G.
de la Torre
,
D. M.
Guldi
, and
T.
Torres
,
Chem. Rev.
110
,
6768
(
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
,
J. Phys. Chem. A
112
,
4215
(
2008
).
5.
H.
Tian
,
Z.
Yu
,
A.
Hagfeldt
,
L.
Kloo
, and
L.
Sun
,
J. Am. Chem. Soc.
133
,
9413
(
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
,
Nat. Commun.
9
,
2038
(
2018
).
7.
A.
Mishra
,
M. K. R.
Fischer
, and
P.
Bäuerle
,
Angew. Chem., Int. Ed.
48
,
2474
(
2009
).
8.
S. M.
Feldt
,
E. A.
Gibson
,
E.
Gabrielsson
,
L.
Sun
,
G.
Boschloo
, and
A.
Hagfeldt
,
J. Am. Chem. Soc.
132
,
16714
(
2010
).
9.
Y.
Zhao
and
W.
Liang
,
Chem. Soc. Rev.
41
,
1075
(
2012
).
10.
P. A.
Liddell
,
G.
Kodis
,
A. L.
Moore
,
T. A.
Moore
, and
D.
Gust
,
J. Am. Chem. Soc.
124
,
7668
(
2002
).
11.
P. A.
Liddell
,
D.
Kuciauskas
,
J. P.
Sumida
,
B.
Nash
,
D.
Nguyen
,
A. L.
Moore
,
T. A.
Moore
, and
D.
Gust
,
J. Am. Chem. Soc.
119
,
1400
(
1997
).
12.
J.
Tinnin
,
S.
Bhandari
,
P.
Zhang
,
H.
Aksu
,
B.
Maiti
,
E.
Geva
,
B. D.
Dunietz
,
X.
Sun
, and
M. S.
Cheung
,
Phys. Rev. Appl.
13
,
054075
(
2020
).
13.
M. H.
Lee
,
B. D.
Dunietz
, and
E.
Geva
,
J. Phys. Chem. C
117
,
23391
(
2013
).
14.
X.
Sun
and
E.
Geva
,
J. Phys. Chem. A
120
,
2976
(
2016
).
15.
X.
Sun
and
E.
Geva
,
J. Chem. Phys.
144
,
044106
(
2016
).
16.
X.
Sun
and
E.
Geva
,
J. Chem. Theory Comput.
12
,
2926
(
2016
).
17.
X.
Sun
and
E.
Geva
,
J. Chem. Phys.
144
,
244105
(
2016
).
18.
X.
Sun
and
E.
Geva
,
J. Chem. Phys.
145
,
064109
(
2016
).
19.
X.
Sun
,
P.
Zhang
,
Y.
Lai
,
K. L.
Williams
,
M. S.
Cheung
,
B. D.
Dunietz
, and
E.
Geva
,
J. Phys. Chem. C
122
,
11288
(
2018
).
20.
Y.
Song
,
A.
Schubert
,
X.
Liu
,
S.
Bhandari
,
S. R.
Forrest
,
B. D.
Dunietz
,
E.
Geva
, and
J. P.
Ogilvie
,
J. Phys. Chem. Lett.
11
,
2203
(
2020
).
21.
R. A.
Marcus
,
J. Chem. Phys.
24
,
966
(
1956
).
22.
R. A.
Marcus
,
J. Chem. Phys.
24
,
979
(
1956
).
23.
R. A.
Marcus
,
Rev. Mod. Phys.
65
,
599
(
1993
).
24.
P. F.
Barbara
,
T. J.
Meyer
, and
M. A.
Ratner
,
J. Phys. Chem.
100
,
13148
(
1996
).
25.
A. J.
Leggett
,
S.
Chakravarty
,
A. T.
Dorsey
,
M. P. A.
Fisher
,
A.
Garg
, and
W.
Zwerger
,
Rev. Mod. Phys.
59
,
1
(
1987
).
26.
N.
Makri
,
J. Phys. Chem. B
103
,
2823
(
1999
).
27.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
102
,
3337
(
1995
).
28.
U.
Weiss
,
Quantum Dissipative Systems
(
World Scientific
,
London
,
2012
).
29.
R.
Kubo
and
Y.
Toyozawa
,
Prog. Theor. Phys.
13
,
160
(
1955
).
30.
J.
Jortner
,
J. Chem. Phys.
64
,
4860
(
1976
).
31.
K. V.
Mikkelsen
and
M. A.
Ratner
,
Chem. Rev.
87
,
113
(
1987
).
32.
M. D.
Newton
and
N.
Sutin
,
Annu. Rev. Phys. Chem.
35
,
437
(
1984
).
33.
M. D.
Newton
,
Chem. Rev.
91
,
767
(
1991
).
34.
D.
De Vault
and
B.
Chance
,
Biophys. J.
6
,
825
(
1966
).
35.
M. H.
Lee
,
E.
Geva
, and
B. D.
Dunietz
,
J. Phys. Chem. C
118
,
9780
(
2014
).
36.
M. H.
Lee
,
B. D.
Dunietz
, and
E.
Geva
,
J. Phys. Chem. Lett.
5
,
3810
(
2014
).
37.
P. L.
Walters
,
T. C.
Allen
, and
N.
Makri
,
J. Comput. Chem.
38
,
110
(
2017
).
38.
F.
Gottwald
,
S. D.
Ivanov
, and
O.
Kühn
,
J. Phys. Chem. Lett.
6
,
2722
(
2015
).
39.
M. K.
Lee
,
P.
Huo
, and
D. F.
Coker
,
Annu. Rev. Phys. Chem.
67
,
639
(
2016
).
40.
M. K.
Lee
and
D. F.
Coker
,
J. Phys. Chem. Lett.
7
,
3171
(
2016
).
41.
S.
Valleau
,
A.
Eisfeld
, and
A.
Aspuru-Guzik
,
J. Chem. Phys.
137
,
224103
(
2012
).
42.
M.
Aghtar
,
J.
Strümpfer
,
C.
Olbrich
,
K.
Schulten
, and
U.
Kleinekathöfer
,
J. Phys. Chem. Lett.
5
,
3131
(
2014
).
43.
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
,
J. Am. Chem. Soc.
120
,
4398
(
1998
).
44.
D.
Kuciauskas
,
P. A.
Liddell
,
S.
Lin
,
S. G.
Stone
,
A. L.
Moore
,
T. A.
Moore
, and
D.
Gust
,
J. Phys. Chem. B
104
,
4307
(
2000
).
45.
D.
Gust
,
T. A.
Moore
, and
A. L.
Moore
,
Acc. Chem. Res.
34
,
40
(
2001
).
46.
N.
Spallanzani
,
C. A.
Rozzi
,
D.
Varsano
,
T.
Baruah
,
M. R.
Pederson
,
F.
Manghi
, and
A.
Rubio
,
J. Phys. Chem. B
113
,
5345
(
2009
).
47.
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
,
Nat. Commun.
4
,
1602
(
2013
).
48.
G.
Su
,
A.
Czader
,
D.
Homouz
,
G.
Bernardes
,
S.
Mateen
, and
M. S.
Cheung
,
J. Phys. Chem. B
116
,
8460
(
2012
).
49.
D.
Balamurugan
,
A. J. A.
Aquino
,
F.
de Dios
,
L.
Flores
, Jr.
,
H.
Lischka
, and
M. S.
Cheung
,
J. Phys. Chem. B
117
,
12065
(
2013
).
50.
A. K.
Manna
,
D.
Balamurugan
,
M. S.
Cheung
, and
B. D.
Dunietz
,
J. Phys. Chem. Lett.
6
,
1231
(
2015
).
51.
O. N.
Starovoytov
,
P.
Zhang
,
P.
Cieplak
, and
M. S.
Cheung
,
Phys. Chem. Chem. Phys.
19
,
22969
(
2017
).
52.
Q.
Shi
and
E.
Geva
,
J. Phys. Chem. A
108
,
6109
(
2004
).
53.
R.
Baer
and
D.
Neuhauser
,
Phys. Rev. Lett.
94
,
043002
(
2005
).
54.
E.
Livshits
and
R.
Baer
,
Phys. Chem. Chem. Phys.
9
,
2932
(
2007
).
55.
T.
Stein
,
H.
Eisenberg
,
L.
Kronik
, and
R.
Baer
,
Phys. Rev. Lett.
105
,
266802
(
2010
).
56.
A. A.
Voityuk
and
N.
Rösch
,
J. Chem. Phys.
117
,
5607
(
2002
).
57.
D. A.
Case
,
R.
Betz
,
D.
Cerutti
,
T.
Cheatham
 III
,
T.
Darden
,
R.
Duke
,
T.
Giese
,
H.
Gohlke
,
A.
Goetz
,
N.
Homeyer
 et al,
AMBER 16
(
University of California
,
San Francisco
,
2016
).
58.
T.
Darden
,
D.
York
, and
L.
Pedersen
,
J. Chem. Phys.
98
,
10089
(
1993
).
59.
G.
Ciccotti
and
J. P.
Ryckaert
,
Comput. Phys. Rep.
4
,
346
(
1986
).