Theoretical estimation of the rate of photoinduced charge transfer reactions in triphenylamine C60 donor–acceptor conjugate

Fullerene‐based molecular heterojunctions such as the [6,6]‐pyrrolidine‐C60 donor–acceptor conjugate containing triphenylamine (TPA) are potential materials for high‐efficient dye‐sensitized solar cells. In this work, we estimate the rate constants for the photoinduced charge separation and charge recombination processes in TPA‐C60 using the unrestricted and time‐dependent DFT methods. Different schemes are applied to evaluate excited state properties and electron transfer parameters (reorganization energies, electronic couplings, and Gibbs energies). The use of open‐shell singlet or triplet states, several density functionals, and continuum solvation models is discussed. Strengths and limitations of the computational approaches are highlighted. The present benchmark study provides an overview of the expected performance of DFT‐based methodologies in the description of photoinduced charge transfer reactions in fullerene heterojunctions. © 2016 Wiley Periodicals, Inc.


Introduction
Fullerenes are probably one of the most promising materials for the construction of organic dyesensitized solar cells (DSSCs). [1][2][3][4][5][6][7][8] The interest in fullerene DSSCs is mainly due to the fact that they are a low-cost replacement for silicon-based solar cells and they are associated to profitable manufacturing and negligible toxicity. [9][10][11][12][13][14] The reference system for fullerene DSSCs is the blend of poly(3-hexylthiophene) (P3HT) and [6,6]phenyl-C61-butyric acid methyl ester (PCBM). It has been extensively studied experimentally and theoretically. [12,[15][16][17][18][19][20][21] For instance, rate constants (kct, subscript ct stands for charge transfer; it can be turned into subscripts cs and cr standing for charge separation and charge recombination, respectively) for photoinduced charge transfer (CT) reactions occurring in P3HT-PCBM reveal the potential of such materials in the construction of efficient DSSCs. Based on quantum chemistry methods, these kct have been demonstrated [22] to be tunable by a few parameters so that the power conversion efficiency may be increased; [23][24][25] which highlights, in fact, the importance of theoretical developments for the understanding and design of DSSCs with improved efficiency. [26] At the molecular level, blends such as P3HT-PCBM are classified into the so called bulk heterojunctions [27] (BHJs), where the donor and acceptor moieties are not covalently connected, and molecular heterojunctions [28] (mHJs) with covalently linked donor and acceptor. The use of mHJs allows a better control of their structure and the charge mobility. Also, mHJs represent a remarkable advantage for computational modeling of underlying CT processes. [28] In particular, the theoretical computation of the different parameters determining kct in fullerene DSSCs can be conveniently performed for these systems.
In this work, we study the mHJ of triphenylamine (TPA) (electron donor) and C60 (electron acceptor) (see Fig. 1) in order to evaluate the performance of density functional theory (DFT) quantum chemical methodology in predicting CT parameters and kct in the donor-acceptor dyads. TPA and its derivatives have been used as photoactive molecules during several decades and their applicability in the construction of DSSCs is also of technological interest. [29][30][31][32][33][34] . The molecules were considered in several experimental and (TD)DFT computational studies related to the mHJ approach. TPA is able to confine cationic charge and hamper aggregation between molecules, which induces self-quenching and reduces the electron injection efficiency. [31,33] On the other hand, C60 is an excellent electron donor [35][36][37][38] that can react with a variety of chemical agents. [39][40][41][42][43][44] Recently, based on electrochemical and photophysical studies, Echegoyen et al. [45] determined kct for the TPA-C60 and TPA-Sc3N@Ih-C80 donor-acceptor conjugates and showed that (i) ultrafast charge separation (CS) reactions occur at the molecular interface, thus giving evidence of CT activity; (ii) TPA-Sc3N@Ih-C80 generates longer-lived photoinduced charge transfer states (CTSs) than does TPA-C60; (iii) increasing the distance of the donor-acceptor interface produces longer-lived CTSs; (iv) nonpolar solvents bring about lack of CT activity; on the contrary, more polar solvents can slow down the charge recombination (CR) reaction; (v) a smaller value of kct for the CR reaction with Sc3N@Ih-C80 in comparison with C60 was found, which was attributed to a smaller value of the driving force.
In this work we aim to assess the performance of DFT-based methodology for the description of CT reactions occurring in the TPA-C60 donor-acceptor complex. In view of that, different schemes to estimate the reorganization energy and driving forces are analyzed. We are aware of the challenges concerning the understanding of CT processes at the molecular level. [46] In particular, electronic couplings and driving forces can be significantly affected by structural fluctuations on the donoracceptor interface and the accurate description of localization/delocalization dynamics of a photogenerated electron-hole pair is quite difficult. Additionally, some studies have underlined a possible role of two-electron excitations in excited states of cyanine-like structures such as TPA, which may represent a limitation for DFT. [47,48] Nevertheless, our results may be a useful reference for future investigations of CT phenomena in (endohedral) fullerene DSSCs.

General
All of the quantum chemical calculations were completed with the program Gaussian09. [49] Unless otherwise specified, the (U)CAM-B3LYP/6-31G* and TD-CAM-B3LYP/6-31G* level of theories [50] coupled to the conductor-like screening model [51] with benzonitrile as the solvent (COSMO:bzn) in equilibrium solvation were used for the calculations of single-point energies, geometry optimizations, and molecular properties of both solvated ground and excited states. Henceforth, let us refer these methods by using the abbreviations U-DFT and TD-DFT, respectively. Despite the good performance of hybrid functionals in the calculation of intramolecular electronic excitations, [52,53] the range separated functional CAM-B3LYP is expected to provide a more reliable description of CTSs. In particular, CAM-B3LYP yields reasonable oscillator strength values and relatively accurate excitation energies for Rydberg and intermolecular CT transitions (within an error of 0.5 eV). [54,55] Although this work is mainly focused on U-DFT and TD-DFT via CAM-B3LYP coupled to COSMO in equilibrium solvation, through the manuscript we discuss the use of other theoretical frameworks, density functionals, and continuum solvation models, as well as the effect of different solvents on excitation energies; all of them are explicitly specified in due course for reasons of clarity.

Structure of the interface TPA-C60
Only one optimized geometry of TPA-C60 was taken into consideration. It corresponds to the groundstate (GS) equilibrium geometry optimized in the implicit presence of benzonitrile. In a previous work, [56] we demonstrated that the structural variability of the TPA-C60 interface modeled by molecular dynamics simulations does not have a significant impact on the intrinsic nature of excited states. Accordingly, all our estimations of excited states are always done at the geometry of GS schematized in Fig. 1, unless otherwise specified.

Analysis of excited states
For the analysis of excited states, [57,58] we make use of the extent of excitation localization in the excited state ψi which can be measured as: where the sums run over all atomic orbitals centered in fragment F, S is the atomic orbital overlap matrix and P 0i is the symmetrized one-electron transition density matrix between the GS ψ0 and the excited state ψi in the atomic orbital representation. The overall charge separation ∆q between fragments F1 and F2 is estimated as: These quantities Xii(F) (being F either C60 or TPA) and ∆q are computed by means of linearresponse [59] TD-CAM-B3LYP/6-31G* in the gas phase.

Rate Constant Expression for Charge Transfer Reactions, kct
Assuming a weak electronic coupling between the donor and acceptor moieties, the theory of nonadiabatic electron transfer is taken into account for the description of the CT processes. In view of that, we use the standard formulation developed by Marcus: [60,61]   where λ is the total reorganization energy, ΔG 0 is the Gibbs energy difference between the initial and final states, Vij is the electronic coupling between the initial and final states, T is the temperature, kB is the Boltzmann constant, and ħ is the reduced Planck constant.

Total Reorganization Energy, λ
The total reorganization energy λ is usually decomposed into two terms, the internal reorganization energy λint and the external reorganization energy λext, so that λ = λint + λext. [60] Internal reorganization energy, λint The energy required to rearrange all the nuclei of the system due to a CT reaction leading from a neutral state to a CTS is computed as follows: where En(CTS) is the energy of the neutral (CTS) state computed in the equilibrium geometry of the neutral (CTS) state and E'n(CTs) is the energy of the neutral (CTS) state computed in the equilibrium geometry of the CTS (neutral state). Nevertheless, λint can be also calculated by considering isolated donor and acceptor fragments as a separated contribution to the internal organization of each fragment: where λD(A) is the reorganization energy of either the donor or the acceptor. En(ion) is the energy of the neutral (ionic) state computed in the equilibrium geometry of the neutral (ionic) state; and E'n(ion) is the energy of the neutral (ionic) state computed in the equilibrium geometry of the ionic (neutral) state. The ionic states are +1ecation and -1eanion states for the donor and acceptor fragments, respectively.

External reorganization energy, λext
Changes in solvent polarization, λext, is one of the parameters that introduces more uncertainty to the estimation of kct. In this work, we present three different approaches to calculate λext as follows: (i) Two-sphere Marcus model, λext (M) [62][63][64] where εo and Δe are respectively the vacuum permittivity and the amount of transferred charge; a1, a2, and R are respectively the donor and acceptor sphere radii and the distance between the sphere centers. We use the contact distance, defined as R = a1 + a2; wherein the sphere radii were calculated by computing the molecular volume of each fragment, defined as the volume inside a contour of 0.001 electrons/Bohr 3 density which was determined from a gas-phase CAM-B3LYP/6-31G* calculation. εOP and εs are the optical and zero-frequency dielectric constants of the surrounding media obtained from the literature. [49] (ii) Non-equilibrium versus equilibrium solvation, λext(Enoneq-Eeq) The difference between non-equilibrium and equilibrium solvation naturally arises from excited state calculations in solution. That is to say, a very fast process occurs when the solvent polarizes its electron distribution (electronic polarization); but the solvent may be reoriented (dipole polarization; e.g. rotation), which is a much slower process. If the solvent does not have enough time to fully respond, like in a vertical electronic excitation, then the situation is described by non-equilibrium solvation (Enoneq). In contrast, if the solvent does have enough time to fully respond to the solute in both ways (electronic and dipole polarizations) then the situation is described by equilibrium solvation (Eeq). Accordingly, λext is simply the energy difference between Enoneq and Eeq. [49] (iii) Dynamic polarization response, λext(Enoneq(ω)-Eeq) In this approach, Eeq is exactly the same as explained in (ii); that is, the time or frequency (ω) to reach Eeq is such that the electronic and dipole polarizations can be completed. In this regard, molecular movements depend upon ω of the change of direction of the electric field; in fact, as ω increases the dipole polarization effect tends to vanish and εs shall be only dependent on the electronic polarization, as a result the non-equilibrium effects shall appear subsequently in the different contributions to the polarization (i.e. orientational, atomic, and electronic polarization). The vanishing dipole polarization effect can be conveniently simulated by setting a calculation in which εs = εOP; thus obtaining Enoneq(ω). Then λext is simply the energy difference between Enoneq(ω) and Eeq. [49,65] Electronic coupling There are a number of studies focused on the development of theoretical frameworks for the evaluation of the electronic coupling. [66][67][68] Previous studies discussed the implications of using a particular level of theory, [69,70] basis set, [71] and methods. [72] We applied the Fragment Charge Difference method (FCD) to derive the coupling values from excited states of the system calculated with TD DFT. [73] Changes in Gibbs energy, ΔGct The charge transfer reactions under study are given as follows: The initial formation of a molecular exciton state is represented by Eq. (8), the CS and CR reactions are respectively Eq. (9) and Eq. (10). Observe that Eq. (8) is not the common approach wherein the photon absorption mainly occurs in the donor fragment. [22,74] The reason is due to the many lowest excited states that are mostly localized in C60. Assuming a negligible entropic component, changes in electronic energy ΔEct for the CT reactions nearly correspond to ΔGct (ΔEct ≈ ΔGct). Franck-Condon excitations and excited-state properties of every excited state represented in Eq. (9) and Eq. (10) were computed via the TD-DFT framework. On the other hand, we also use the U-DFT formalism to obtain Franck-Condon and relaxed CTSs by interchanging one alpha HOMO localized in TPA with one alpha LUMO located in C60 such that a self-consistent field solution is searched; [49] in other words, constraining the self-consistent field process to a solution of the type TPA + -C60 -. To this end, we first perform a RCAM-B3LYP/6-31G* calculation to obtain the equilibrium geometry and relaxed wavefunction of GS; then, we use the Gaussian09 standard keyword Guess=(Read,Alter) in order to obtain either a Franck-Condon or relaxed CTS by following the next process: (i) Guess=Read recovers the GS wavefunction, (ii) Guess=Alter specifies orbitals to be alternated from the GS wavefunction (from Figure 2, it should be understood that interchanging alpha HOMO by alpha LUMO must lead to a CT situation, that is TPA + -C60 -), and (iii) we perform a UCAM-B3LYP/6-31G* calculation keeping the modified wavefunction so as to reach a self-consistent field solution for both cases: the geometry-relaxed or the Franck-Condon CTS. As mentioned before, we make use of the abbreviation U-DFT for every calculation performed in the unrestricted formalism for either CTSs or some other exciton state. That is, by following the same procedure, neutral excited states can be also modeled by U-DFT by interchanging alpha HOMO-1 and alpha LUMO since this alteration does not lead to a CT situation in view of both orbitals are fully localized in C60 (see TPA-C60* in Eq. 10).
Observe that the restricted formalism is only used for GS.

Results and Discussion
This section initiates with the introduction of the electronic transitions under study with their respective oscillator strength value and orbital contributions, along with the extent of exciton delocalization Xii(F) and amount of charge separation Δq, so as to provide a comprehensive description of the excited states of interest. Then, a detailed evaluation of reorganization energies and electronic couplings is reported; wherein the external reorganization energy is the main subject to discuss in view of the difficulties associated to the calculation of such a parameter. To account for Gibbs energies, we firstly aim to describe the initial formation of a molecular exciton from which charge separation can occur; after that Gibbs energies of charge separation and recombination reactions are given. The variation of the dielectric constant as a function of the excitation energy is also analyzed so that the lack of charge transfer activity in nonpolar solvents is explained. In line with this subject, we briefly discuss alternative approaches for the description of CT activity, such as the use of other density functionals, other continuum solvation models, triplet states, or even other theoretical approximations. Finally, the most relevant results are collected to calculate the rate constants for the photoinduced charge transfer reactions under consideration; then the results are compared to the experimental evidence. At the end, a brief justification of the lack of nuclear degrees of freedom in the current study is explicated.

Extent of charge separation and exciton delocalization for electronically excited states
Our discussion begins with the analysis of the most relevant electronic excited states. In Fig. 2 some frontier molecular orbitals (MOs) in TPA-C60 are depicted. As a matter of fact, most frontier MOs are entirely localized at the fullerene cage; therefore, the lowest-energy electronic transitions are expected to be mainly localized at the fullerene moiety. An electronic transition exclusively occurring in the cage leads to the formation of a localized excited state (LES). On the other hand, since HOMO is fully localized at the TPA fragment, an electronic transition in which the main contribution is a HOMO to LUMO+X (in C60) excitation leads to the formation of a CTS. In Table 1 the two lowest in energy LESs and CTSs are reported as calculated via gas-phase TD-CAM-B3LYP/6-31G*. The HOMO-1-to-LUMO transition leads to the excited state lowest in energy (LES1) that lies at 2.46 eV with respect to the gas-phase GS. It is an orbitally forbidden transition as indicated by the vanishing value of the oscillator strength. Furthermore, we can observe that TPA has no contribution in the formation of LES1 because this molecular exciton state is not extended through this fragment, Xii(TPA) = 0.0. However, this state is clearly originated from a transition fully localized at the C60 moiety as shown by the Xii(C60) = 1.0 value and the lack of charge separation (Δq = 0.0). Exactly the same conclusions are attained for LES2; which are in good agreement with the experimental spectrum of C60, wherein the 2.0 to 3.0 eV energy range corresponds only to forbidden transitions. [75,76] In view of that, other LESs are expected to be inaccessible by photon absorption within the mentioned energy range. The excitation corresponding to a HOMO-to-LUMO transition lies at 3.48 eV and leads to the CTS lowest in energy. Nevertheless, CTSs can be only accessed upon dissociation of a Frenkel exciton; from which an electron may be migrated to a phase boundary through spatial translation of such an exciton. As a result, electronic transitions leading to a CTS by direct photon absorption are either completely forbidden or only weakly allowed. Observe that this is indeed the case for CTS1 and CTS2 in view of their vanishing oscillator strengths. These exciton states with strong charge transfer character are mainly characterized by the amount of donor-to-acceptor charge separation, Δq = 1.0, and their notable exciton delocalization through the whole complex; that is to say, CTS1 and CTS2 are neither localized at TPA nor C60 since Xii(TPA) and Xii(C60) are both null.
Further discussion is focused in the electronic excited states reported in Table 1. For our purposes, CT reactions can be consistently described by those states; although the reader should keep in mind that other states of different nature might contribute to the CS process; like those exciton states with partial charge separation where some amount of charge (Δq ~ 0.5) is transferred from the donor to the acceptor. Such partial CTSs were analyzed in detail by us for the TPA-C60 dyad, [56] although other works describing such partial CTSs in fullerene DSSCs can be found in the literature. [77][78][79] In addition, conformational changes in TPA-C60 were demonstrated to have no significant impact in the intrinsic nature of excited states, therefore the analysis done in this subsection in a single (the most stable) conformation for the GS remains valid. [56]

Internal and external reorganization energies
The appropriate selection of CTSs is crucial for the evaluation of the CT parameters affecting the rate constant. Nuclei are rearranged to compensate the structural stress caused by the excess of charge generated through the entire system. However, geometry relaxation of CTSs is a challenging task for large-sized systems; as a matter of fact, in TPA-C60, excited states are scarcely separated by ca. 0.1 eV (see complete details in Table S1), which hinders the geometry optimization of a CTS since the procedure rather leads to a neutral excited state because LESs are more abundant. In view of that, we make use of a low-cost methodology to perform a geometry relaxation of CTSs with no geometrical constraints; i.e. open-shell singlet CTSs are relaxed by constraining an excess of charge in the donor and acceptor fragments. This strategy can be accomplished by U-DFT as explained in the Methodology section; where we have reached a geometry-relaxed CTS by interchanging alpha HOMO in TPA by alpha LUMO in C60 (see Fig. 2); [49] then the calculation of λint through Eq. (4) can be conveniently done. In Table 2 we compare this result with the ones obtained through Eq. (5). In this latter case, we optimized the isolated fragments C60, C60 -, TPA, and TPA + in the gas phase and in solution. Observe that there is no significant difference between all λint results; therefore we use the average value for λint for the estimation of kct in view of the small standard deviation.  Table 2 lists also λext values. Unlike λint, the uncertainty of the estimations of λext is large since our results diverge within a 0.4 eV energy range. As mentioned before, the reliable evaluation of this parameter is difficult to achieve. Some authors set the value of λext as an adjustable parameter; [80,81] on the other hand, the direct evaluation of λext can be consistently performed via QM/MM [82] or a polarizable force field. [83] It is expected that QM/MM would provide a more accurate value for λext, yet it is required that the polarizability of the molecular fragments can be validated through experiments; a condition that is not fulfilled for the system under study. Additionally, the lack of detailed information on the interface and surroundings makes us consider the methods reported in Table 2; indeed, we keep in mind that, in some cases, neither implicit nor explicit solvent models are sufficient enough to reproduce experimental observations as shown for paracyclophane chromophores and, therefore, one might expect different result from different modeling of the solvent. [84] However, it is worth to mention that the different approaches that we use to calculate λext are completely different (two-sphere against continuum solvation models), but the missing reproducibility must be solved; this issue shall be analyzed at the end of the manuscript.

Electronic couplings
The calculated values of electronic couplings are listed in Table 3. Because the values are sufficiently small, the non-adiabatic approach, Eq. (3), may be applied to estimate the corresponding kct. The coupling is rather sensitive to localization of excited state density. In particular, because of the coupling difference, the LES1 → CTS1 charge separation is expected to be faster than LES1 → CTS2 by a factor of ~ 25.

Initial formation of a molecular exciton state upon light absorption
The initial formation of LES1 represented in Eq. (8) by direct light absorption corresponds to a dipoleforbidden transition as discussed in section 3.1. Dipole-allowed LESs of higher-energy populated by light absorbance are expected to relax very fast to LES1. In order to prove that charge separation can indeed occur from LES1 quantum dynamics simulations are required. [59] This type of calculation is out of the scope of the present work. There are some experimental and theoretical works wherein real-time dynamics of CS states in fullerene DSSCs were described. [85,86] For our purposes we highlight the study carried out by Li et al. [87,88] where simulations based on real-time dynamics were performed for a C60 fragment (the acceptor) covalently linked to an aniline derivative (the donor). Regardless the oscillator strength, they show the population of LUMO, depopulation of HOMO and population/depopulation of HOMO-1 (coincidently, those MOs resemble the ones in Fig. 2) in a femtosecond (fs) scale after photoinduction. The interesting point of their study is that the population of LUMO (in C60) by one electron can be actually originated from a HOMO-1→LUMO transition; however, while the time evolves HOMO-1 (the photogenerated hole in C60) is populated again at the time HOMO (in the donor) is depopulated by one electron bringing about charge separation, a process that occurs coherently as the hole in C60 is generated/dissociated in an oscillating way (see Fig. 3). Moreover, they observed that the CR reaction takes place after ca. 80 fs as a result of the declination of the LUMO population. In higher-energy and short-time scales they also analyzed the initial molecular excitation formed in the donor fragment, where the charge separation barely lasted 5 fs. In both processes the electron-hole recombination was demonstrated to be dominant. Based on the similitudes in the molecular and electronic structure of the system studied by Li et al. and ours (the donors are cyanine-like structures linked to a pyrrolidine-like ring as the bridge, and the acceptor is C60), we extrapolate their results to our system and infer that charge separation can occur from the initial formation of LES1.

Charge transfer reactions
The relative energy of LES1 lies at 2.68 eV as calculated through gas-phase linear-response TD-CAM-B3LYP/6-31G* and 2.42 eV as computed via state-specific (COSMO:bzn)TD-CAM-B3LYP/6-31G* in equilibrium solvation (recall that these and subsequent Franck-Condon excitation energies are relative to the energy of GS; unless otherwise specified). Continuum solvation models implemented under the state-specific framework are claimed to be more accurate than the linearresponse approach; we refer the reader to the seminal papers dealing with this subject. [89][90][91][92][93][94] Unfortunately, the improvements of time-dependent state-specific calculations could not be observed in the calculation of CTSs since we had many difficulties trying to obtain a physically reasonable excitation energy for a solvated CTS in TPA-C60, with no success (CTSs become more stable than the GS). In view of that, we had to resort to the U-DFT approach in order to obtain the Gibbs energy of solvation ΔGsol of CTSs.
The U-DFT Franck-Condon excited states reported in Table 4 were modeled by moving one alpha electron of the HOMO-1 (in C60) or HOMO (in TPA) to the LUMO+X (in C60). We verified that dipole moments (μ), summations of all the atomic Mulliken charges on a fragment (QF, where the values +1, 0 or -1 erespectively means that the complete fragment F is in a cation, neutral or anion state), and the expected values of <S 2 > correspond to the excited state under consideration (see Table  S2). That is to say, in neutral states, μ for GS are 4.1 and 4.8 debye in the gas phase and in solution, respectively; and analogously 5.1 and 6.7 debye for LES1. Besides, we distinguish three fragments: the donor TPA, the bridge, and the acceptor C60. For a neutral state, QF are +0.1, 0.0 and -0.1 efor TPA, the bridge, and C60, respectively. On the other hand, μ for CTSs ranges from 47.7 to 49.4 and from 56.1 to 59.5 debye in the gas phase and in solution, respectively. Moreover, QF are +1.0, 0.0 and -1.0 efor TPA, the bridge, and C60; which confirms the correct simulation of the CTSs since TPA appears as a cation and C60 as an anion. In all the cases, the open-shell singlet character was also verified since <S 2 > is ca. 1.05 (see also the frontier MOs depicted in Table S3 for every excited state under consideration). U-DFT therefore predicts the average value of -1.58 eV for ΔGsol of CTSs with a standard deviation of 0.03 eV. We use this value for the evaluation of the solvated TD-DFT CTS1. The energy of CTS1 lies at 3.69 eV as calculated through gas-phase linear-response TD-CAM-B3LYP/6-31G*; as a result the energy of CTS1 in solution lies at 2.11 eV (compare this value with the one of U-DFT CTS1 in Table 4). In Figure 3 we provide a schematization of our calculations and CT processes. Our results lead to a ΔGcs of -0.31 eV, which is exactly the experimentally reported value. [45] This value is obtained from the 2.42 eV corresponding to LES1 obtained with state-specific (COSMO:bzn)TD-CAM-B3LYP/6-31G* and the 2.11 eV of the CTS1 calculated from the gas-phase TD-CAM-B3LYP/6-31G* together with the U-DFT solvation energy of the singlet open-shell CTS. However, the neutral singlet excited state energy is reported to be 1.76 eV, which results in an overestimation of 0.66 eV when compared to the energy calculated for LES1. As a matter of fact, ΔGcr = 2.11 eV is also overestimated as compared to the experimental value by exactly the same 0.66 eV. [45] This overestimation constitutes another source of error for the theoretical prediction of kcr. This issue shall be analyzed in detail at the end of the manuscript together with λext.

Quenched charge transfer reactions
The solvent plays an important role in the CT activity. So far, all our calculations are performed by taking into account benzonitrile (εs = 25.6). Nonetheless, in Fig. 4 the variation of the excitation energy in solution of LES1 and CTS1 as a function of εs is examined via U-DFT. As εs decrease, CTS1 becomes less stable than LES1 and the CS reaction does not occur because it corresponds to an endergonic process. These results concur with the experimental evidence in which steady-state fluorescence measurements confirmed the absent of CT activity in toluene (εs = 2.4) and carbon disulfide (εs = 2.6). Moreover, CT activity can be observed even in tetrahydrofuran (THF, εs = 7.4) despite of the slightly favorable driving force, which is calculated to be only 0.02 eV (0.09 eV experimental); unlike less polar solvents. [45]  The use of hybrid functionals for gas-phase TD-DFT calculations should be avoided since they predict CTS1 to be the lowest-energy excited state, which is in contradiction with Fig. 4 and experimental results. Indeed, the similar results obtained through linear-response TD-B3LYP/6-31G* coupled to the PCM, [65] COSMO, and SMD [95] solvation models suggest that, when compared to the gas-phase results, failures in the determination of CTS energies are not originated from the continuum solvation model, but from the inherent problems of the theory (see Table S4). The use of range separated functionals are supposed to be a better alternative in the determination of CTSs. [96] Accordingly, our results at the TD-ωB97XD/6-31G* level of theory [97] were compared to the analogous TD-CAM-B3LYP ones. Unlike hybrid functionals, range separated functionals do reproduce the correct order of excited state energies since CTSs are calculated to be lying ~1.5 eV above the lowest localized excited state in the gas phase; in agreement with the U-DFT ΔGsol estimated for CTSs in benzonitrile (see Table S5). This trend was also observed by using the Tamm-Dancoff approximation [59] (see Table  S1). Unfortunately, solvent effects could not be observed by means of linear-response TD-DFT coupled to continuum solvation models since energies in solution result in exactly the same quantity as in the gas phase; and in the case of state-specific TD-DFT, CTS energies were highly underestimated so that negative excitation energies were calculated [49] (i.e. CTSs were determined to be lying below the GS). Both approaches therefore could not reproduce the trends observed in Fig. 4 despite the use of range separated density functionals. We assume that a possible reason of the failure of linear-response TD-DFT to calculate solvated CTSs is due to the vanishing predicted transition dipole moments.
As a final remark, it has been purposed the alternative use of U-DFT triplet excited states instead of TD-DFT in the determination of excitation energies and properties of CTSs in fullerene DSSCs. [22,79] However, the formation of low-energy triplet excited states in donor-acceptor pairs have been proved to be a loss mechanism; [98,99] which may be attributed to the dipole-forbidden transitions. The justification of the use of U-DFT triplet states lies in the fact that it is a low-cost alternative to reach reliable energies of CTSs in fullerene DSSCs; that is, the unpaired electrons are well separated such that the exchange interaction is negligible, therefore singlet and triplet CTSs are expected to be nearly structurally and energetically equivalent. In fact, the relative energy with respect to GS of a geometryrelaxed triplet CTS is 1.88 eV; only 0.15 eV lower in energy than the relaxed singlet CTS1. Moreover, the triplet CTS has λint = 0.16 eV in good agreement with our previous results for singlet CTSs (see Table S6 for more details).

Rate constants of photoinduced charge transfer reactions
In previous sections we discussed the evaluation of λint, λext, Vij, and ΔGct concerning different drawbacks related to the intrinsic problems of the theoretical methodologies under consideration. We were able to anticipate the sources of error with a bigger impact in the estimation of kct; which are λext and ΔGcr (i.e. the energy of CTS1); therefore, let us conduct this section around these two parameters.
In Table 5, kct is estimated as a function of λext. In the case of the CS reaction, the value of ΔGcs is in excellent agreement with the experimental reference; as a result, it seems that a two-sphere model (i.e. λext(M)) might be the best alternative in the estimation of kct since our theoretical values resemble the experimental value of kcs. However, to avoid misleading observations, we consider that consistent conclusions related to the appropriate evaluation of λext can be only attained by simultaneously taking into account both CT reactions. In this regard, by using a two-sphere model, kcr drastically vanishes; as a result this approach is not suitable for the calculation of kct in the system under consideration. ΔGcr also introduces uncertainty. In Table 5, kcr is reported as a function of ΔGcr. We note that our estimation of ΔGcr can be used together with λext(Enoneq(ω)-Eeq) so as to better reproduce the experimental trends of kcr within deviations of less than two order of magnitudes. Similarly, the experimental value of ΔGcr can be used together with λext(Enoneq-Eeq). Those observations suggest that solvent reorganization can be better simulated [49] by continuum solvation models than by the twosphere model. The limitation of the latter may be attributed to the fact that the condition R >> a1+a2 is not fulfilled for our system.
As a final remark, it might be proposed that the lack of structural variability could alter our estimations; and also the use of the Marcus-Levich-Jortner equation might be a better approach to include the effect of quantum vibrational modes, wherein the nuclear relaxation caused by CT is expected to involve intramolecular (high-frequency) modes. [100,101] Nevertheless, we infer that the role of nuclear degrees of freedom in the system under study can be ignored and still the physics of the CT process can be described satisfactorily. Our deduction is supported by the fact that we have previously described the intrinsic nature of excited states for a number of thermally accessible structures of the TPA-C60 interface as simulated by molecular dynamics; from which we concluded that geometrical variations do not significantly change the intrinsic nature and distribution of excited states. [56] As a matter of fact, through the manuscript we show the similitudes in structure and energetics between Franck-Condon and geometry-relaxed CTSs. Moreover, in an interesting discussion by Troisi, [102] it is suggested that CS reactions in organic photovoltaic interfaces are mostly purely electronic processes which is in line with our conclusions.

Conclusions
We have assessed the performance of U-CAM-B3LYP and TD-CAM-B3LYP with the 6-31G* basis set for estimating the rate constant of photoinduced charge transfer reactions using the molecular heterojunction triphenylamine C60 donor-acceptor conjugate as the reference system with the following conclusions: (i) Electronic couplings have been computed using multi-configurational excited states as generated by TD-DFT. The applied scheme does not use orbital approximation to represent locally excited and CT states and thus provides more accurate values.
(ii) Franck-Condon and relaxed charge transfer states can be conveniently modeled by a U-DFT approach coupled to a continuum solvation model to explain experimental data or predict trends in CT properties of related systems. Due to the large separation of the unpaired electrons, either singlet or triplet states with charge transfer character can be used to estimate the internal and external reorganization energies as well as the energetics of charge transfer reactions. On the other hand, both linear-response and state-specific solvation models were unable to provide a reasonable estimate of the energy of solvated charge transfer states.
(iii) The evaluation of external reorganization energies can be properly performed by continuum solvation models. The use of the two-sphere model may be limited for systems with relatively small donor-acceptor distance.
(iv) The lack of charge transfer activity experimentally observed in nonpolar solvents is attributed to destabilization of charge transfer states relative to the lowest excited state localized on the dye (TPA).
This result can be only reproduced by time-dependent range separated functionals (CAM-B3LYP).
In addition, a good performance of U-DFT has been found: the obtained variation of excitation energies as a function of the dielectric constant reproduces well the experimental data.
(v) We have provided a benchmark for the strengths and weaknesses of computational schemes used for modelling photoinduced charge separation and charge recombination processes in fullerene based molecular heterojunctions.