Excited state electron and hole transfer underpin fundamental steps in processes such as exciton dissociation at photovoltaic heterojunctions, photoinduced charge transfer at electrodes, and electron transfer in photosynthetic reaction centers. Diabatic states corresponding to charge or excitation localized species, such as locally excited and charge transfer states, provide a physically intuitive framework to simulate and understand these processes. However, obtaining accurate diabatic states and their couplings from adiabatic electronic states generally leads to inaccurate results when combined with low-tier electronic structure methods, such as time-dependent density functional theory, and exorbitant computational cost when combined with high-level wavefunction-based methods. Here, we introduce a density functional theory (DFT)-based diabatization scheme that directly constructs the diabatic states using absolutely localized molecular orbitals (ALMOs), which we denote as Δ-ALMO(MSDFT2). We demonstrate that our method, which combines ALMO calculations with the ΔSCF technique to construct electronically excited diabatic states and obtains their couplings with charge-transfer states using our MSDFT2 scheme, gives accurate results for excited state electron and hole transfer in both charged and uncharged systems that underlie DNA repair, charge separation in donor–acceptor dyads, chromophore-to-solvent electron transfer, and singlet fission. This framework for the accurate and efficient construction of excited state diabats and evaluation of their couplings directly from DFT thus offers a route to simulate and elucidate photoinduced electron and hole transfer in large disordered systems, such as those encountered in the condensed phase.

Photoinduced electron transfer (ET) is a fundamental step in important photochemical processes ranging from carrier generation and transport in photovoltaic materials such as organic solar cells and dye–semiconductor heterostructures to photodamage of biomolecules such as DNA and RNA and photo-initiated catalytic processes including reduction of CO2 by photoactive transition metal complexes, water oxidation at photoelectrodes, and hot electron transfer from plasmonic nanoparticles to molecules. Charge-localized diabatic states (diabats) provide a chemically intuitive basis to simulate and elucidate these photoinduced processes since they possess the desirable feature that their chemical character is robust to changes in nuclear configurations. Indeed, diabats form the basis of widely used theories of charge transfer, such as Marcus–Hush theory,1–3 where diabatic state energies and couplings are essential to determine reaction rates. These diabatic states, especially when it is feasible to compute their associated nuclear forces, also permit a fully dynamical description of excited state charge-transfer processes that provides mechanistic details beyond those accessible from rate theories. Accurate and efficient ways to construct such diabats and calculate their couplings therefore open the door to studying the nonadiabatic dynamics of out-of-equilibrium photoinduced processes occurring across multiple time and length scales.

In photoinduced electron transfer processes, the reactant state typically features a photo-excitation on the donor molecule, while the product state features electron–hole separation. The former can be physically represented by a locally excited (LE) diabat and the latter by a charge-transfer (CT) diabat. Density functional theory (DFT)4,5 offers an appealing approach to constructing these diabats due to their computationally efficient treatment of electron correlation. One commonly used approach that can be applied to both wavefunction-based theories and DFT is computing a set of adiabatic states and then transforming them into diabats through a unitary rotation. This adiabatic-to-diabatic (ATD) transformation can be constructed by defining diabats as eigenstates of operators that are closely related to the ET process6–12 or states that are maximally localized on the donor or acceptor moieties13–16 and offer the advantage that they provide a route to construct diabatic states in cases where the chemical process of interest does not involve a physically intuitive charge separation between distinct chemical units. The most commonly employed approach to generate excited states from DFT is through linear-response time-dependent density functional theory (LR-TDDFT).17–19 However, due to self-interaction error, LR-TDDFT underestimates the energy of charger-transfer excited states18,20 that are essential to photoinduced ET processes. Although these issues can somewhat be alleviated by range-separated hybrid (RSH) functionals21–24 and system-specific ω-tuning techniques,25,26 they often require balancing improving the accuracy of CT-type excitation energies with artificial overestimation of LE state energies. Furthermore, while the application of ATD diabatization to ground-state electron or hole transfer (ET/HT) typically only requires computing the ground and the first excited adiabatic states, in photoinduced ET processes many higher-energy excited states of distinct characters can contribute to the relevant diabats. Since ATD-based approaches, such as the generalized Mulliken–Hush (GMH) method,6,7 can usually only be applied to a few adiabatic states, these must be carefully selected through a comprehensive excited-state analysis to achieve reliable results. In addition, diabatic states constructed through an ATD procedure are typically not variationally optimized, and the transformation matrix itself varies with the nuclear positions, which complicates access to forces needed for geometry optimization and nuclear dynamics on diabatic potential energy surfaces.

These drawbacks of ATD-based diabatization schemes have motivated the development of DFT-based approaches that create diabatic states without first computing the excited-state adiabats. However, whereas for ground-state ET and HT there are many DFT-based methods to directly construct diabats and evaluate their couplings,27–37 fewer approaches exist to generate the LE and CT states and their couplings needed to describe photoinduced ET. One such method38 constructs LE states from an LR-TDDFT calculation and the CT states using constrained DFT (CDFT)39,40 and evaluates LE–CT couplings using the CDFT configuration interaction (CDFT-CI) prescription.31,41 However, its applicability is limited by its strict requirement that the global system TDDFT excited state be well-localized on the donor moiety such that it can be directly treated as the LE state. To enforce the locality of the LE diabat, one can perform an LR-TDDFT calculation on the donor fragment in isolation and then construct the global system diabat by forming an antisymmetrized product with the ground state wavefunctions of the other fragments. A recent work42 used this procedure followed by a Thouless transformation38,43 to construct the LE state and then constructs the CT diabat using the block-localized Kohn–Sham (BLKS)35,44 approach. However, the use of two different theories (LR-TDDFT and BLKS) for the construction of the LE and CT states means that these are not constructed on an equal footing and the use of LR-TDDFT for the LE state means that this state is not variationally relaxed.

Here, we present a new approach to construct variationally optimized LE and CT diabats using absolutely localized molecular orbitals (ALMOs) and demonstrate its accuracy in treating photoinduced electron and hole transfer in model complexes inspired by excited-state processes in biochemical and photovoltaic applications. Our approach exploits ΔSCF45,46 to target molecular excited states while maintaining a balanced accuracy for different types of electronic excited states46,47 as well as the ability of ALMOs37,48 to produce diabatic states that are variationally optimized at the supersystem KS-DFT level. We then couple these diabatic states using a generalization of our ALMO(MSDFT2) scheme, which we have previously demonstrated to give errors below 5% over a broad range of ground state electron/hole transfer chemical systems.37 The resulting approach circumvents the trade-off between the accuracy of the LE and CT states in TDDFT, and the use of variationally optimized ALMOs allows one to easily obtain gradients with respect to nuclear configurations that are not easily accessible using ATD-type diabatization methods. We then show how our new method performs in comparison to benchmark values obtained from excited state calculations at the equation-of-motion coupled-cluster singles and doubles (EOM-CCSD)49,50 level of theory followed by GMH diabatization6,7 on a variety of systems relevant to DNA photo-damage/protection, charge separation in photovoltaics, photoinduced electron transfer from the chromophore to solvent, and singlet fission. By doing this, we show that our approach provides an accurate and efficient method for constructing DFT-based diabatic states and evaluating their couplings for photoinduced electron and hole transfer that can easily provide access to nuclear forces for quantum dynamics simulations.

To provide a DFT-based method to create the electronic states that govern photoinduced electron transfer processes, here we show how to construct locally excited (LE) and charge-transfer (CT) diabats and compute their couplings using ALMOs. We have previously shown37 for electron and hole transfer in systems in their ground electronic states that using ALMOs prevents charge delocalization between donor and acceptor moieties under the Mulliken definition of charge population,51 making it suitable for the construction of charge-localized diabats. Here, we consider photoinduced electron transfer in a donor–acceptor complex DA, where the donor (D) has nD electrons and the acceptor (A) has nA electrons. The ground state of this system, |DA⟩, can be represented as

|DA=NdetϕD1,,ϕDi,ϕDnD,ϕA1,,ϕAnA,
(1)

where “det” denotes the Slater determinant and N=1/(nD+nA)! is the normalization constant. Each molecular orbital (MO) in the determinant is assigned to either the donor or the acceptor fragment and is expanded using atomic orbital (AO) basis functions located only on that fragment, yielding orbitals that are “absolutely localized.”

For systems in their ground electronic states that can be readily partitioned into donor and acceptor moieties, one can use a “bottom-up” approach to construct the ALMO diabats.37 This consists of performing SCF calculations for the isolated fragments D and A first and then generating the initial guess to the ALMO state by assembling the resulting fragment orbitals via Eq. (1). This state is also known as the frozen state in ALMO-based energy decomposition analysis.52,53 One can then variationally optimize the ALMOs by minimizing the Kohn–Sham (KS) energy functional of the full system, EKS[P], with respect to the occupied-virtual mixings on each fragment, where P is the one-particle density matrix (1PDM) that can be constructed from the occupied ALMOs (Co) through

P=Co(σoo)1Co,σoo=CoTSCo.
(2)

Here, σoo is the overlap metric of occupied orbitals, which is obtained by transforming the AO overlap matrix (S) into the ALMO basis.

To create diabats with the donor moiety in its electronically excited state, one can utilize the ΔSCF approach45,46 within ALMO-based DFT calculations. The ΔSCF approach has been shown to accurately capture both valence and charge-transfer excited states in cases where the target excited state is dominated by the transition between a specific pair of orbitals.46,47 This procedure is typically applied by taking the electronic ground state orbitals and then promoting electrons from the occupied to the virtual MOs to generate a non-aufbau electronic configuration. The MOs are then relaxed within this non-aufbau electronic configuration, resulting in an SCF solution that corresponds to an electronic excited state. To prevent the variational optimization from collapsing to the ground state, one can employ the maximum overlap method (MOM)54 or its initial MOM (IMOM)47 variant. In these methods, the occupied orbitals at each SCF iteration are selected to be the ones that have the largest overlap with the span of the occupied orbitals in the previous (MOM) or initial (IMOM) non-aufbau electronic configuration. Hence, they are compatible with SCF algorithms performing diagonalization of the Fock matrix, such as the widely used direct inversion of the iterative subspace (DIIS)55 method.

To generate a variationally optimized LE state from the fragment orbitals of the excited-state donor and ground-state acceptor, one brings together these orbitals and then relaxes them via an additional ΔSCF-type procedure, yielding a global system ALMO state that corresponds to the mutually relaxed LE diabat |D*A⟩ (illustrated in the top row of Fig. 1). One begins by calculating the ground state of the acceptor using standard SCF and the excited state of the donor using ΔSCF and then concatenates the fragment orbitals obtained in the form of a single Slater determinant given by

|D*A0=NdetϕD1(0),,ϕDa(0),ϕDnD(0),ϕA1(0),,ϕAnA(0),
(3)

which we denote as |D*A0 and corresponds to the “Initial” unrelaxed state in the top row of Fig. 1. This initial LE state corresponds to an excitation on the donor dominated by the transition from its ith occupied orbital to ath virtual orbital.

FIG. 1.

Procedure of our ALMO-based approach for the calculation of LE–CT diabatic couplings. Top row: construction of the LE diabat |D*A⟩ using the ALMO-ΔSCF method. Bottom row: construction of the CT diabat |D+A⟩ using a ground state ALMO calculation. The solid blocks represent unrelaxed fragment orbitals and the shaded blocks represent orbitals variationally relaxed in the full system.

FIG. 1.

Procedure of our ALMO-based approach for the calculation of LE–CT diabatic couplings. Top row: construction of the LE diabat |D*A⟩ using the ALMO-ΔSCF method. Bottom row: construction of the CT diabat |D+A⟩ using a ground state ALMO calculation. The solid blocks represent unrelaxed fragment orbitals and the shaded blocks represent orbitals variationally relaxed in the full system.

Close modal

To variationally relax the ALMOs in the LE state, we then solve the locally projected SCF (Stoll) equation48,56 for each fragment,

ISP+SPXTFIPS+PXSXXCX=SXXCXϵX.
(4)

Here, X ∈ {D, A} is the fragment index, FF[P] is the AO Fock matrix constructed from the 1PDM (P), PX is the projector for fragment X, and ϵX and CX are the diagonal matrix of eigenvalues and the matrix composed of eigenvectors, which correspond to the orbital energies and AO expansion coefficients of ALMOs on fragment X, respectively. ϵX and CX are obtained by diagonalizing the locally projected Fock matrix, defined as the quantity in the square bracket on the left-hand side of Eq. (4). In Stoll’s equation, the projector for fragment X has the form

PX=Coσoo1XCo,XT.
(5)

Here, the subscript X denotes matrix columns that belong to the partition corresponding to fragment X. Explicitly, when X = D, one takes columns 1 to nD of the resulting matrix, and when X = A, one takes columns nD + 1 to nD + nA. The matrices [Co(σoo1)]X and Co,X are thus the contravariant and covariant versions of the occupied ALMO coefficients belonging to fragment X, respectively.48,57 The outer product of [Co(σoo1)]X and Co,X then gives the Stoll projector to fragment X.48,56 While iterating the Stoll equation [Eq. (4)] to self-consistency with only fragments in their ground states yields the relaxed ALMOs, when any of the fragments is in its excited state, this procedure alone would lead to the collapse of this local excited state to the ground state. Hence, the Stoll iteration must be performed in concert with a procedure to prevent this excited state collapse.

Here, we employ the IMOM procedure47 to prevent the excited state from collapsing to the ground state when solving the Stoll equation for the excited fragment D*. This procedure requires one to calculate the overlap (O) between ALMOs on fragment D* obtained at a given iteration and the initial occupied MOs on the same fragment, denoted as CD*,o(0),

O=CD*,o0TSCD*.
(6)

To select the orbitals to be occupied for the relaxed fragments, one has to choose the orbitals that have the largest overlap with the initial set of occupied orbitals in the non-aufbau electronic configuration. These orbitals are selected as the nD orbitals that possess the largest projection onto the span of the initial set of occupied orbitals, PD*j,

PD*j=i=1nDOD*iD*j21/2,
(7)

where i and j denote the molecular orbitals in the initial and current iterations, respectively. The selected orbitals are then employed to generate the updated 1PDM. Iterating Eq. (4) with IMOM to self-consistency, one obtains the relaxed LE state, denoted as the “Relaxed” state in the top row of Fig. 1,

|D*A=NdetϕD1(LE),,ϕDa(LE),ϕDnD(LE),ϕA1(LE),,ϕAnA(LE),
(8)

where the superscript (LE) indicates that the ALMOs are optimized within the LE state in contrast to the unrelaxed fragment orbitals in |D*A0. We denote the procedure to variationally optimize ALMOs in the excited state introduced in the previous paragraphs as ALMO-ΔSCF in the following discussion.

To prepare the lowest-energy intermolecular CT diabat, we employ the procedure introduced in our previous work.37 This CT diabat usually corresponds to transferring an electron from the highest occupied molecular orbital (HOMO) of the donor to the lowest unoccupied molecular orbital (LUMO) of the acceptor, which we denote as |D+A⟩. One can prepare this diabat by performing an unrestricted ALMO calculation in this charge-separated state.57 Starting from SCF calculations for fragments D+ and A in isolation, we construct the CT diabat and then relax its orbitals using the standard ALMO-based SCF procedure, yielding

|D+A=NdetϕD1(CT),,ϕDnD1(CT),ϕA1(CT),,ϕAnA(CT),ϕAnA+1(CT).
(9)

This procedure is illustrated in the bottom row of Fig. 1. As in Eq. (8), the superscript (CT) indicates that the ALMOs are relaxed within the specific electronic configuration of the CT state. While the diabat given in Eq. (9) corresponds to an aufbau CT configuration, one could generate higher-energy CT states by employing the ALMO-ΔSCF procedure used to prepare the LE states described above.

To evaluate the electronic coupling between LE and CT diabats, D*A|Ĥ|D+A, one can utilize the original multistate DFT (MSDFT) approach35,36 or the MSDFT2 scheme that we proposed in our previous work.37 Defining |ψa⟩ ≡ |D*A⟩ and |ψb⟩ ≡ |D+A⟩, one can construct the diabatic Hamiltonian in the non-orthogonal ALMO basis,

H=HaaHabHbaHbb.
(10)

To extract the electronic coupling, Hab=D*A|Ĥ|D+A, between these two ALMO diabats, one transforms H′ into an orthogonal basis. Using Löwdin’s symmetric orthogonalization scheme,58 one obtains H=S1/2HS1/2, where the interstate overlap matrix is [S]xyψx|ψy and x, y ∈ {a, b}. In the two-state case, this yields

Hab=11Sab2HabHaa+Hbb2Sab,
(11)

where Sab is given by

Sabψa|ψb=detCoaTSCob.
(12)

Here, Co(a) and Co(b) are ALMO coefficient matrices associated with the two diabats and S is the AO overlap matrix. In the following discussion, we use Soo(ab)(Co(a))TSCo(b) to denote the overlap matrix between occupied orbitals from diabats |ψa⟩ and |ψb⟩.

In both MSDFT35,36 and MSDFT237 approaches, the diagonal elements of H′ are the KS energies of each diabat, i.e., Haa=EaKS[P(a)] and Hbb=EbKS[P(b)], where P(a) and P(b) are the 1PDMs associated with |ψa⟩ and |ψb⟩, respectively. These two schemes differ in the way they approximate the off-diagonal element Hab. In the original MSDFT scheme, Hab is given by a KS-DFT correction to the Hartree–Fock (HF) interstate coupling59,60 result,

Hab=SabVnn+Pabh+12PabIIPab+12(ΔEac+ΔEbc),
(13)

where the first three terms on the right-hand side arise from the nuclear repulsion energy (Vnn), one-electron Hamiltonian (h), and electron-repulsion integrals (II) as in the HF theory. Pab is the transition density matrix between diabats |ψa⟩ and |ψb⟩,

Pab=Co(a)Soo(ba)1Co(b)T.
(14)

The last term in Eq. (13) is a correction intended to incorporate the contribution from the exchange-correlation (XC) functional to the interstate coupling, which is given by the average of the difference between the KS and HF energies of each diabat calculated using their individual 1PDMs,

ΔEac=EaKSP(a)EaHFP(a),
(15)
ΔEbc=EbKSP(b)EbHFP(b).
(16)

As shown in our previous work,37 the treatment of interstate coupling in the original MSDFT approach [Eq. (13)] effectively couples the KS determinants using only the HF Hamiltonian once substituted into Eq. (11) and is thus equivalent to how the interstate couplings are evaluated in some recently introduced ab initio exciton models.61–63 Our MSDFT2 approach explicitly incorporates the XC contribution in the interstate coupling by utilizing the KS energy functional of the symmetrized transition density matrix (P̃ab),37 

Hab=SabVnn+Pabh+12PabIIPab+ExcP̃ab,
(17)

where P̃ab=(Pab+Pba)/2. Importantly, the electron-repulsion integrals, II, in Eq. (17) account for Coulomb integrals as well as a fraction of HF exchange when hybrid functionals are used, which differ from those in Eq. (13) where II denotes the full HF Coulomb and exchange integrals.

A notable feature of the MSDFT2 approach is that the off-diagonal element of the diabatic Hamiltonian given by Eq. (17) reduces to the KS energy of a diabat when a = b, ensuring internal consistency of the theory. The same formula [Eq. (17)] was also employed to evaluate the coupling between DFT-based diabats in the frozen density embedding (FDE) method for ground-state ET/HT.33,34,64 We have previously shown37 that for ground-state ET and HT processes, the ALMO(MSDFT2) method can accurately predict diabatic couplings between donor and acceptor states and that its performance is superior to the original MSDFT approach and other popular DFT-based diabatization schemes.27–32 In this work, we assess the ability of both MSDFT2 and MSDFT in coupling LE diabats constructed from ALMO-ΔSCF calculations with ALMO CT states. The overall schemes are denoted as Δ-ALMO(MSDFT2) and Δ-ALMO(MSDFT) in the following discussion.

In the following, we discuss the weak-coupling regime where the diabatic coupling between two diabats approaches zero, which presents a challenge when calculating the diabatic coupling using ALMO(MSDFT2). We then provide a way to determine when one is in this regime and a physically intuitive alternative coupling scheme that yields stable results in this limit. Two situations that lead to weak coupling are as follows: when the overlap between two diabats approaches zero due to the increase in the intermolecular distance and when there is destructive interference between the orbitals involved in the charge transfer process.

The former situation, which leads to near-zero overlap between the orbitals relevant to the charge transfer process, precludes the construction of the transition density matrix between two ALMO diabats. In particular, the matrix inversion of the overlap matrix between the diabats’ occupied orbitals (Soo(ba)) required in Eq. (14) becomes numerically unstable when the matrix becomes near-singular. In the latter situation, the destructive interference between the two diabats results in underestimated off-diagonal elements of the non-orthogonal diabatic Hamiltonian (|Hab|). This is most likely due to the amplified effects of the self-interaction error in this regime since the results become very sensitive to the choice of the XC functional in these cases. When this situation causes a sign flip upon Löwdin orthogonalization of the ALMO diabats using Eq. (11), i.e., |Hab|<|Sab(Ea+Eb)/2|, the coupling value obtained becomes unreliable.

To address this challenge, we exploit an analogy to the generalized Slater–Condon rules to suggest an alternative that is applicable in the weak-coupling limit. As we noted previously, the MSDFT2 scheme can be viewed as a KS-DFT analogue to the generalized Slater–Condon rule58,60 for two non-orthogonal HF states (determinants) that have a non-zero overlap.37 Based on this realization, here, we propose a KS-DFT analogue to the generalized Slater–Condon rule for two non-orthogonal determinants whose interstate occupied orbital overlap matrix possesses one vanishing singular value.

We achieve this by first symmetrically orthogonalizing the occupied orbitals in two ALMO-based diabats, which yield C¯o(a) and C¯o(b). We then perform a singular value decomposition (SVD) on the interstate occupied orbital overlap matrix, S¯oo(ab),

S¯oo(ab)C¯o(a)TSC¯o(b)=UsVT,
(18)

and then use the U and V matrices obtained to transform the orthogonal orbitals C¯o(a) and C¯o(b) to the Löwdin-paired orbitals,58 C̃o(a)=C¯o(a)U and C̃o(b)=C¯o(b)V. We proceed to order the Löwdin-paired orbitals using their singular values such that the first vector in each set of orbitals (C̃o(a) and C̃o(b)) corresponds to the pair of orbitals that overlap the least (associated with the smallest singular value, s1), denoted here as C̃o,1(a) and C̃o,1(b). The remaining N − 1 vectors in each set, C̃o,N1(a) and C̃o,N1(b), correspond to the remaining occupied orbitals (where N denotes the total number of electrons) and sN−1 corresponds to the set of N − 1 singular values. This results in two re-ordered sets of Löwdin-paired orbitals: C̃o(a)=[C̃o,1(a),C̃o,N1(a)] and C̃o(b)=[C̃o,1(b),C̃o,N1(b)].

With this notation in place, we obtain the electronic coupling between diabats |ψa⟩ and |ψb⟩ as

Hab=S̃abP1(ab)FKS(ab)PN1(ab)=S̃abP1(ab)h+IIPN1(ab)+VxcP̃N1(ab),
(19)

where

P1ab=C̃o,1aC̃o,1bT,
(20)
PN1ab=C̃o,N1asN11C̃o,N1bT
(21)

are two transition density matrix-like objects, P̃N1(ab)=(PN1(ab)+PN1(ba))/2 is the symmetrized version of PN1(ab), and S̃ab=i1Nsi is the reduced interstate overlap. The KS Fock matrix (FKS) is constructed from the transition density matrix between N − 1 pairs of orbitals, which consists of contributions from the core-Hamiltonian (h), two-electron integrals (II), and the exchange-correlation potential (Vxc), as shown in the second line of Eq. (19).

Our prescription in Eq. (19) for calculating the diabatic coupling in the weak-coupling regime is a KS-DFT analogue to the generalized Slater–Condon rule for two HF determinants whose interstate occupied orbital overlap, Soo(ba), has one zero singular value,58,60

Hab=S̃abP1(ab)h+IIPN1(ab),
(22)

where II denotes the contribution from the Coulomb interaction and full HF exchange. The approximation central to the application of our KS-DFT analogue to this Slater–Condon rule is that one can employ Eq. (19) when the overlap between orbitals C̃o,1(a) and C̃o,1(b) is small but nonzero, which is the case in the weak-coupling regime.

In the following, we refer to this approach to evaluate the diabatic coupling in the weak-coupling regime given in Eq. (19) as the MSDFT2-wc scheme. We apply MSDFT2-wc in the cases discussed above where either orbital overlap is vanishingly small or destructive interference between the diabats occurs. Our identification of the weak-coupling regime in which Eq. (19) is used to compute the diabatic coupling is based on the magnitude of the smallest singular value in s and a “sign-flipping” criterion and is described in detail in Sec. III.

We extended our previous implementation of ALMO(MSDFT2)37 in Q-Chem65 to construct the diabats involved in photoinduced ET processes using Δ-ALMO(MSDFT2). To construct the LE diabat |D*A⟩, one first obtains the ground-state SCF solutions for both fragments. A non-aufbau electronic configuration of the donor is then prepared by swapping the occupied and virtual orbitals to reflect the nature of the excitation that is being constructed. As shown in Fig. 1, one first computes the ΔSCF solution for the excited donor (D*) and then uses the orbitals on D* and ground-state A to construct the initial guess for the full-system LE state [Eq. (3)], which is followed by a variational optimization via the ALMO-ΔSCF procedure introduced in Sec. II A. To prevent the collapse onto the ground electronic state, the IMOM method47 was employed in combination with the DIIS algorithm55 in both the fragment ΔSCF and the global system ALMO-ΔSCF calculations. The CT states involved in this work correspond to the transition from the donor’s HOMO to the acceptor’s LUMO and therefore were all constructed from ground-state unrestricted ALMO calculations. We note that the LE and CT states in this work are spin-symmetry broken, i.e., the commonly adopted spin-projection scheme45 for singlet states in ΔSCF is not applied. Since spin-projected (LE and CT diabatic) states are no longer variationally optimized and the same broken-symmetry states were used in the previous studies using the MOM or IMOM methods and provided good accuracy,47,54 the choice to work with spin-symmetry broken states provides a convenient way to compute LE-CT couplings, which would require the evaluation of additional matrix elements if spin-projection were employed.

The MSDFT and MSDFT2 coupling schemes, which we have implemented in the released version of the Q-Chem 5.3 package,65 are as introduced in our previous work.37 To evaluate the off-diagonal element of the non-orthogonal diabatic Hamiltonian (Hab), our implementation first orthogonalizes the occupied orbitals constituting two ALMO diabats separately and then transforms these orbitals to be biorthogonal using Löwdin’s orbital pairing scheme,58 allowing one to identify the pair of orbitals having the smallest overlap. Our implementation identifies the weak-coupling regime based on two criteria: when the smallest singular value of the interstate overlap is smaller than 10−4 and/or when “sign-flipping” occurs in Eq. (11), i.e., when |Hab|<|Sab(Ea+Eb)/2|. When either of these conditions occurs, we employ the MSDFT2-wc scheme given in Eq. (19) to evaluate the diabatic coupling. We have employed the standard MSDFT2 scheme that uses the full sets of Löwdin-paired orbitals from both diabats (C̃o(a) and C̃o(b)) to construct the transition density matrix in all other cases. For calculations based on MSDFT, we have employed the original prescription given in Eq. (13) in all scenarios.

To diabatize the adiabatic states obtained from TDDFT calculations, we employed the generalized Mulliken–Hush (GMH)6,7 and fragment charge difference (FCD)8 approaches. These methods first require one to select a set of adiabatic states that are relevant to a given electron transfer process. One then constructs the adiabatic-to-diabatic transformation by defining diabatic states as eigenstates of a given molecular property: dipole moment in GMH and difference in fragment charge populations in FCD.

The GMH method utilizes the dipole operator,6,7 which requires calculation of both the dipole moment of each adiabat and the transition dipole between each pair of them. When only two adiabats (denoted with indices “1” and “2”) are considered, the coupling between two resulting diabats (denoted with indices “a” and “b”) is given by

Hab=|μ^12|(E2E1)μ1μ22+4μ^1221/2,
(23)

where E1 and E2 are the energies of the two adiabats, μ1 and μ2 are their dipole moments, and μ^12 denotes the projection of the transition dipole (μ12) along the charge-transfer direction: e^CT=(μ1μ2)/|μ1μ2|.

The FCD method,8 on the other hand, defines the diabats as eigenstates of the fragment charge difference matrix (ΔQij = Qij(D) − Qij(A), where D denotes the donor and A denotes the acceptor). The calculation of the diagonal and off-diagonal elements of ΔQ requires performing charge population analysis with the electron density of each adiabat and the transition density between each pair of them. In the two-state case, the diabatic coupling is given by

Hab=|ΔQ12|(E2E1)(ΔQ11ΔQ22)2+4ΔQ122.
(24)

Here, we evaluated the differences in fragment charge populations subject to the Mulliken charge population scheme.51 

The reference values for the diabatic couplings were obtained from GMH diabatization of EOM-CCSD states. Unless otherwise specified, we used the two-state version of the GMH and FCD diabatization schemes given by Eqs. (23) and (24). The procedures for the three-state GMH diabatization required for the indole–guanine complex and FCD for a general number of states that were used to diabatize the TDDFT adiabats for the pentacene dimer are described in Secs. S1 and S2 of the supplementary material, respectively.

The calculations in this work were all performed with our locally developed version of the Q-Chem 5.3 package.65 Unless otherwise specified, all the DFT-based calculations for diabatic couplings (based on ALMO or TDDFT) were performed with the ωB97X-D functional66 and the 6-31+G(d) basis set67,68 on a (99, 590) grid, which has 99 radical shells and 590 Lebedev points in each. To test the performance of each method on the choice of the functional, we also generated results for the indole–guanine complex (Sec. IV A) with pure GGA (BLYP69,70 and PBE71), global hybrid (B3LYP72 and PBE073), and two other range-separated hybrid (RSH) functionals (CAM-B3LYP22 and LRC-ωPBEh24). Using the same system, we also investigated the basis set dependence of the diabatic couplings calculated from ALMO- and TDDFT-based methods and revealed that the results are largely insensitive to the choice of basis set (see Fig. S1 of the supplementary material). Both the standard and ALMO-based SCF calculations were converged to a DIIS error below 10−8 a.u. The linear-response TDDFT calculations were performed within the Tamm–Dancoff approximation,74 and only the singlet excited states were considered. The reference values for the diabatic couplings were obtained from diabatization of the adiabatic states obtained from EOM-CCSD50,75 calculations, where the IP (ionization potential) variant was applied to the indole–guanine complex while the EE (electronic excitation) variant was applied to the others. Both EOM-CCSD variants employ a closed-shell reference for the systems that we investigate here, and Hartree–Fock (HF) theory was used as the default method to calculate the reference state. While it is often possible to identify the characters of excited states using the canonical orbitals from HF calculations, allowing one to identify the excited states relevant to a photoinduced ET process, the HF virtual orbitals of the 1-naphthol–CHCl3 complex (Sec. IV C) are highly diffuse, rendering it difficult to identify the relevant states contributing to the ET process. To address this issue, we generated the reference orbitals from a ground-state SCF calculation using the LRC-ωPBEh functional for this system, which resulted in canonical KS orbitals that allow for transparent assignment of relevant excited states.

The geometries of the systems used here were obtained from previous studies: the cationic indole–guanine complex from Ref. 9, naphthalene–tetracyanoethylene complex from Ref. 25, 1-naphthol–CHCl3 complex from Ref. 76, and the pentacene dimer from Ref. 77. The scan of inter-monomer distances and monomer rotation angles for the pentacene dimer was performed with respect to the “reference” dimer structure given in Ref. 77 in which dCC (shown in the inset of Fig. 10) is 5.98 Å. For completeness, we provide these geometries in a ZIP file available in the supplementary material.

Here, we benchmark a set of excited-state electron and hole transfer systems to demonstrate the performance of Δ-ALMO(MSDFT2) in obtaining the diabatic couplings involved in these photoinduced processes by comparing to the reference values obtained from a GMH diabatization of states computed at the EOM-CCSD level. Because TDDFT is commonly employed to calculate adiabatic excited states, we also compare our Δ-ALMO(MSDFT2) results to diabatic couplings obtained from the application of adiabatic-to-diabatic (ATD) transformation schemes (GMH and FCD) to the TDDFT states. The four systems on which we benchmark these methods were chosen to highlight four fundamentally important applications of excited-state electron and hole transfer as well as the specific challenges associated with their diabatization.

In particular, we first consider ground- and excited-state hole transfer between indole and guanine, which is relevant to DNA repair mechanisms and where the approaches based on the ATD transformation require more than two states to obtain accurate results. We then treat photoexcited electron transfer in naphthalene–tetracyanoethylene, which mimics the photoinduced charge separation process within a donor–acceptor dyad in organic electronics and illustrates a case where Δ-ALMO(MSDFT2) can be used when the LE state involves more than one orbital transition. Next, we focus on the photoexcited electron transfer in the 1-naphthol–CHCl3 complex, which serves as a paradigm of chromophore-to-solvent electron transfer and demonstrates the need to use the MSDFT2-wc coupling scheme introduced in Sec. II C in the weak-coupling regime. Finally, we address the diabatic coupling between the LE and CT states in the pentacene dimer, which is relevant for the superexchange mechanism that has been shown to be essential for efficient singlet fission42,77,78 and presents a challenge to the commonly used schemes based on ATD transformation such as GMH due to the involvement of multiple adiabatic excited states and lack of a uniform charge-transfer direction.

The cationic indole–guanine complex, [Ind-G]+, serves as a prototypical example of hole transfer from nucleobases to amino acid residues in DNA–protein complexes, a process that can inhibit the oxidative damage of DNA.9,79 Here, we consider two hole transfer processes from indole, whose HOMO and HOMO−1 orbitals are close in energy, to guanine. These two processes are illustrated in Fig. 2 where for the ground state process, the hole is on the HOMO of indole (Ind+) and for the excited state process, it is on the HOMO−1 (Ind+*). We note that physically, the reverse process, i.e., hole transfer from guanine to indole, is relevant to DNA protection. However, considering the forward process allows us to use the language that aligns with that of Sec. II A, with the Ind+-G and Ind+*-G configurations assigned as the ground (GS) and locally excited (LE) initial states, respectively, and their diabatic couplings to the CT state denoted as HGS-CT and HLE-CT. Since the diabatic states involved in the forward and reverse processes are identical, the diabatic couplings are the same.

FIG. 2.

The ground- and excited-state hole transfer processes in the cationic indole–guanine ([Ind-G]+) complex. The hole is transferred from the HOMO (ground-state HT) or HOMO-1 (excited-state HT) of indole to the HOMO of guanine. These orbitals are obtained from ALMO-based diabats for Ind+-G and Ind+*-G and visualized with an isovalue of 0.02 a.u.

FIG. 2.

The ground- and excited-state hole transfer processes in the cationic indole–guanine ([Ind-G]+) complex. The hole is transferred from the HOMO (ground-state HT) or HOMO-1 (excited-state HT) of indole to the HOMO of guanine. These orbitals are obtained from ALMO-based diabats for Ind+-G and Ind+*-G and visualized with an isovalue of 0.02 a.u.

Close modal

To obtain the reference values for the GS–CT and LE–CT couplings, we first calculated the three lowest-lying adiabatic states for the cationic system using EOM-IP-CCSD and then performed GMH diabatization on them (see Sec. III). The EOM-IP-CCSD calculation starts from the closed-shell charge-neutral reference system. The three lowest-lying resultant adiabatic states are dominated by ionization from the highest three occupied orbitals (HOMO, HOMO−1, and HOMO−2) of the global system. These orbitals, as shown in the middle column of Fig. 3, are approximately linear combinations of orbitals in the diabatic picture (ALMOs), namely, the HOMO and HOMO-1 of indole and the HOMO of guanine. The delocalized nature of the adiabatic orbitals and particularly the entangled character of HOMO−1 of the global system (in which all three diabatic orbitals possess non-negligible weights) make it necessary to invoke three-state GMH diabatization80 rather than using the more commonly used two-state approximation (see Secs. S1 and S3.A in the supplementary material for the detailed procedure). The reference values obtained are HGS-CT = 277 meV and HLE-CT = 60 meV.

FIG. 3.

Depiction of the connection between the adiabatic and diabatic molecular orbitals involved in the ground- and excited-state hole transfer in the cationic [Ind–G]+ complex. The adiabatic orbitals that are plotted (shown in the middle) are obtained from the closed-shell HF reference for the EOM-IP-CCSD calculation, and the diabatic orbitals plotted (shown in the two sides and also in Fig. 2) are obtained from ALMO-based diabats for Ind+–G and Ind+*–G. The orbitals are visualized with an isovalue of 0.02 a.u.

FIG. 3.

Depiction of the connection between the adiabatic and diabatic molecular orbitals involved in the ground- and excited-state hole transfer in the cationic [Ind–G]+ complex. The adiabatic orbitals that are plotted (shown in the middle) are obtained from the closed-shell HF reference for the EOM-IP-CCSD calculation, and the diabatic orbitals plotted (shown in the two sides and also in Fig. 2) are obtained from ALMO-based diabats for Ind+–G and Ind+*–G. The orbitals are visualized with an isovalue of 0.02 a.u.

Close modal

Since TDDFT is commonly used to investigate excited-state processes, we also calculated the GS–CT and LE–CT couplings for the HT in this system using TDDFT followed by three-state GMH and FCD diabatization. Unlike EOM-IP-CCSD, the TDDFT calculations employ the unrestricted KS-DFT description of the cationic complex as the reference. With global hybrid (B3LYP and PBE0) and range-separated hybrid (RSH) functionals (CAM-B3LYP, LRC-ωPBEh, and ωB97X-D), the first TDDFT excited state is dominated by a local excitation involving a transition from indole’s HOMO−1 to HOMO, while the second excited state exhibits marked charge transfer from guanine to indole. We assign the LE and CT characters to these states using the detachment and attachment densities81 associated with them, as shown in Fig. S2 of the supplementary material. When pure GGA functionals (BLYP and PBE) are employed, the state featuring G → Ind CT is shifted to become the fourth excited state.

The resulting GS–CT and LE–CT couplings obtained from applying the adiabatic-to-diabatic transformation schemes, GMH and FCD, to the TDDFT states are shown in Fig. 4. For the GS–CT coupling shown in Fig. 4(a), both TDDFT/GMH and TDDFT/FCD yield similar results that systematically overestimate HGS-CT. The relative errors are over 100% when paired with pure and global hybrid functionals and reduce to 30%–40% when RSH functionals are employed. While Fig. 4(b) demonstrates that GMH and FCD also give similar results for the LE–CT coupling when using RSH and global hybrid functionals with relative errors of 40%–50%, they differ dramatically when paired with the two GGA functionals. In particular, when using BLYP or PBE, TDDFT/GMH gives more accurate LE–CT couplings, likely due to fortuitous error cancellation, whereas TDDFT/FCD significantly overestimates the coupling by ∼150%. Hence, diabatization based on TDDFT excited states is ill-suited for treating hole transfer in this explicitly charged system.

FIG. 4.

Diabatic couplings for the (a) ground and (b) excited state hole transfer in the [Ind–G]+ complex calculated using TDDFT with the 3-state GMH and FCD diabatization schemes as well as with two ALMO-based approaches. The calculations are performed with density functionals ranging from pure GGA to global and range-separated hybrid functionals. The horizontal dashed lines mark the reference values for HGS-CT and HLE-CT, obtained at the EOM-IP-CCSD/6-31+G(d) level, and the vertical dotted lines indicate the grouping of pure GGA, global hybrid, and range-separated hybrid functionals.

FIG. 4.

Diabatic couplings for the (a) ground and (b) excited state hole transfer in the [Ind–G]+ complex calculated using TDDFT with the 3-state GMH and FCD diabatization schemes as well as with two ALMO-based approaches. The calculations are performed with density functionals ranging from pure GGA to global and range-separated hybrid functionals. The horizontal dashed lines mark the reference values for HGS-CT and HLE-CT, obtained at the EOM-IP-CCSD/6-31+G(d) level, and the vertical dotted lines indicate the grouping of pure GGA, global hybrid, and range-separated hybrid functionals.

Close modal

In contrast to the large overestimation of TDDFT combined with GMH and FCD schemes, ALMO(MSDFT2) gives accurate results (with the lowest error being 5% and the largest one being 17%) for the diabatic coupling of the GS–CT process, even when combined with pure GGA functionals. This observation is consistent with our previous work37 where we observed for a range of ground-state ET and HT processes that ALMO(MSDFT2) could give highly accurate results even when combined with lower-tier functionals. Whereas the Δ-ALMO(MSDFT2) scheme that we have introduced here suffers from poor performance for the LE–CT coupling when paired with pure functionals producing an error of 50%–60%, it yields accurate results when combined with global hybrid and RSH functionals, with relative errors of 10%–20% that are significantly lower than those obtained using the other methods, except for TDDFT/GMH when combined with GGA functionals, which has very likely benefited from fortuitous cancellation of error. In contrast to MSDFT2, using the MSDFT coupling scheme leads to systematic overestimation of both the GS–CT and LE–CT diabatic couplings when combined with any tier of functionals, whose results degrade markedly in the excited-state process with errors of 45%–85%. The systematic overestimation of MSDFT couplings is in agreement with the trend revealed by our previous ground-state ET and HT benchmarks.37 Thus, when combined with global hybrid or RSH functionals, ALMO(MSDFT2) provides a good balance of general applicability and accuracy in obtaining the diabatic couplings for the ground- and excited-state hole transfer processes in this system.

The naphthalene–tetracyanoethylene (Np–TCNE) complex is a prototypical model system for a donor–acceptor dyad resembling those in photovoltaic devices. In contrast to the [Ind–G]+ complex investigated above, this complex possesses a charge-neutral closed-shell ground state. The naphthalene molecule has two low-lying singlet excited states, which we denote as LE1 and LE2 (also denoted as 1La and 1Lb states in the literature,82 respectively). A TDDFT calculation for naphthalene using an RSH functional ωB97X-D reveals that the LE1 state is dominated by the HOMO → LUMO transition, while the LE2 state a 1:1 mix of the HOMO−1 → LUMO and HOMO → LUMO+1 transitions, as illustrated in Fig. 5. The calculations reveal that the two states are 4.66 eV and 4.92 eV above the ground state, respectively. The charge transfer excitation from the HOMO of naphthalene to the LUMO of TCNE has a lower energy (2.6 eV)25 and thus can be accessed via photoinduced electron transfer from the LE states on naphthalene.

FIG. 5.

The excited-state electron transfer processes in the naphthalene–tetracyanoethylene (Np-TCNE) complex. The electron can be transferred from Np in the LE1 or LE2 locally excited state to TCNE giving rise to a CT state corresponding to an electron transfer from Np’s HOMO to TCNE’s LUMO. The former process is symmetry-forbidden, while the latter is symmetry-allowed. The LE1 state is dominated by the HOMO → LUMO transition. The LE2 state consists of a 1:1 mix of the HOMO−1 → LUMO and HOMO → LUMO+1 transitions, which has non-zero coupling with the CT state.

FIG. 5.

The excited-state electron transfer processes in the naphthalene–tetracyanoethylene (Np-TCNE) complex. The electron can be transferred from Np in the LE1 or LE2 locally excited state to TCNE giving rise to a CT state corresponding to an electron transfer from Np’s HOMO to TCNE’s LUMO. The former process is symmetry-forbidden, while the latter is symmetry-allowed. The LE1 state is dominated by the HOMO → LUMO transition. The LE2 state consists of a 1:1 mix of the HOMO−1 → LUMO and HOMO → LUMO+1 transitions, which has non-zero coupling with the CT state.

Close modal

The equilibrium geometry of this complex is of C2v symmetry, and the inter-fragment distance (defined as the distance between the midpoints of the central C–C bonds in each monomer) is 3.9 Å. Owing to this symmetry, the LE1 state and the lowest CT state belong to the B2 and B1 irreducible representations, respectively, causing their coupling to vanish. This zero-coupling result is exactly reproduced (to numerical precision) by both the Δ-ALMO(MSDFT) and Δ-ALMO(MSDFT2) approaches. The LE2 state, on the other hand, belongs to the B1 irreducible representation leading to a non-zero coupling with the CT state. Using EOM-EE-CCSD to compute the first two B2 states followed by a two-state GMH diabatization gives a reference value of 128 meV for the diabatic coupling between the LE2 and CT states.

The LE2 state of Np, in principle, cannot be described by standard ΔSCF-based methods since it corresponds to a linear superposition of two equally weighted transitions: HOMO → LUMO+1 and HOMO−1 → LUMO. However, exploiting the equal weights of these two configurations in the LE2 state, one can calculate their respective couplings with the CT state first and then evaluate the coupling between the LE2 and CT states using

HLE-CT=[H(H1L,CT)+H(HL+1,CT)]/2.

The resulting LE–CT couplings obtained by Δ-ALMO(MSDFT) and Δ-ALMO(MSDFT2) are 164 meV and 127 meV, respectively (see Table S3 of the supplementary material for the two intermediate couplings between a single ALMO-ΔSCF state and the CT state), with the latter in excellent agreement with the EOM-EE-CCSD reference (128 meV).

We also obtained the results of TDDFT-based diabatization for the Np–TCNE complex at the equilibrium geometry (see Table S4 of the supplementary material for details about the obtained TDDFT states). Using the lowest two B1 excited states as the adiabatic basis, the result of FCD diabatization (126 meV) also agrees very well with the EOM-EE-CCSD reference although TDDFT underestimates the energy of the lowest B1 (CT) state by ∼1 eV (see Table S4), while GMH underestimates the coupling slightly (114 meV). The improved accuracy of approaches based on the ATD transformation of TDDFT states for this closed-shell system is in stark contrast with the more inaccurate results obtained for the open-shell [Ind–G]+ complex investigated in Sec. IV A.

We now investigate the ability of TDDFT- and ALMO-based approaches to capture the distance dependence of the LE–CT coupling over inter-fragment separations from 3.5 Å to 5.0 Å with fixed monomer geometries. The results in Fig. 6 show that all methods examined are able to qualitatively capture the exponential decay of the LE–CT coupling, with Δ-ALMO(MSDFT2) and TDDFT/FCD capturing the benchmark result within graphical accuracy. In contrast, Δ-ALMO(MSDFT) systematically overestimates the coupling with the relative error growing from 24% to 48% as the inter-fragment distance increases, whereas TDDFT/GMH underestimates the coupling over the entire range with an (unsigned) relative error increasing from 9% to 18% with the increase in distance. The increase in the size of the relative error with the increase in donor–acceptor separation implies that, in addition to inaccurately capturing the diabatic couplings, these methods also incorrectly capture the rate of decay of the couplings with the donor–acceptor distance. In contrast, Δ-ALMO(MSDFT2) accurately captures the diabatic couplings and their distance dependence for excited state ET in the Np–TCNE complex, which forms a prototypical example of donor–acceptor dyads.

FIG. 6.

Distance dependence of the LE–CT coupling (in meV) for the photoinduced ET in the naphthalene–TCNE complex obtained from TDDFT and ALMO-based calculations. The LE state corresponds to the LE2 state of naphthalene, and the CT state corresponds to the transition from naphthalene’s HOMO to TCNE’s LUMO. The calculations were performed with inter-fragment distances equal to 3.5 Å, 4.0 Å, 4.5 Å, and 5.0 Å and the reference values were obtained at the EOM-EE-CCSD/6-31+G(d) level.

FIG. 6.

Distance dependence of the LE–CT coupling (in meV) for the photoinduced ET in the naphthalene–TCNE complex obtained from TDDFT and ALMO-based calculations. The LE state corresponds to the LE2 state of naphthalene, and the CT state corresponds to the transition from naphthalene’s HOMO to TCNE’s LUMO. The calculations were performed with inter-fragment distances equal to 3.5 Å, 4.0 Å, 4.5 Å, and 5.0 Å and the reference values were obtained at the EOM-EE-CCSD/6-31+G(d) level.

Close modal

The 1-naphthol–CHCl3 complex is a prototypical system exhibiting chromophore-to-solvent electron transfer, a process that can play an important role in the non-radiative decay pathway of excited states in the condensed phase. The ultrafast photoinduced ET from 1-naphthol (NpOH) to CHCl3 was recently investigated using Marcus theory combined with TDDFT/FCD diabatization.76 Here, we examine the ability of Δ-ALMO(MSDFT2) to capture the LE–CT couplings involved in this process. Importantly, this example also allows us to demonstrate how the Δ-ALMO(MSDFT2) approach can be used in the weak coupling regime using the approaches introduced in Sec. II C.

A TDDFT calculation of the NpOH–CHCl3 complex shows that the two lowest-lying singlet excited states are dominated by local excitations on NpOH: the first is dominated by the HOMO → LUMO transition and the second by the HOMO → LUMO+1 transition. In the following discussion, we denote these locally excited states LE1 and LE2, respectively. Unlike the Np–TCNE complex investigated above, the energy of the lowest CT state in this system, which corresponds to the transition from the HOMO of NpOH to the LUMO of CHCl3, is higher than the energies of both of the LE states. Because the energy of the CT state predicted by EOM-EE-CCSD is markedly higher than those from TDDFT, which lifts the CT state up to the seventh singlet excited state of this system (with TDDFT, it is the third singlet excited state), it becomes necessary to include a relatively large number (∼10) of states in the EOM-EE-CCSD calculation. The energies of these relevant adiabatic excited states given by TDDFT and EOM-EE-CCSD are provided in Table S5 of the supplementary material.

The two photoinduced ET processes considered in this work, corresponding to the LE1 → CT and LE2 → CT pathways, are illustrated in Fig. 7. The LE2 → CT pathway involves the LE2 state that is low-lying in energy (based on EOM-EE-CCSD calculations) but is only weakly coupled to the CT state due to the destructive overlap between NpOH’s LUMO+1 and CHCl3’s LUMO. In contrast, the LE1 → CT pathway involves the higher-energy LE1 state that is more strongly coupled with the CT state since the symmetry of the LUMOs of NpOH and CHCl3 allows constructive orbital overlap (shown in Fig. 7). The reference values that we obtained using EOM-EE-CCSD/GMH for the LE1–CT and LE2–CT couplings are 72 meV and 17 meV, respectively.

FIG. 7.

The excited state electron transfer processes between 1-naphthol (NpOH) and CHCl3. The electron can be transferred from NpOH in two different locally excited states (LE1 and LE2) to CHCl3, giving rise to a CT state that corresponds to the transition from the HOMO of NpOH to the LUMO of CHCl3. The two LE states are dominated by the HOMO → LUMO and HOMO → LUMO+1 excitations, respectively. The LE2 state is weakly coupled to the CT state (indicated by the dashed arrow) due to the mismatched symmetry of NpOH’s LUMO+1 and CHCl3’s LUMO.

FIG. 7.

The excited state electron transfer processes between 1-naphthol (NpOH) and CHCl3. The electron can be transferred from NpOH in two different locally excited states (LE1 and LE2) to CHCl3, giving rise to a CT state that corresponds to the transition from the HOMO of NpOH to the LUMO of CHCl3. The two LE states are dominated by the HOMO → LUMO and HOMO → LUMO+1 excitations, respectively. The LE2 state is weakly coupled to the CT state (indicated by the dashed arrow) due to the mismatched symmetry of NpOH’s LUMO+1 and CHCl3’s LUMO.

Close modal

We now turn to the performance of TDDFT- and ALMO-based schemes in predicting the LE1–CT and LE2–CT couplings, presented in Fig. 8. The TDDFT/FCD results are in good agreement with the reference values in general, while TDDFT/GMH underestimates both couplings with a relative error of 18% for LE1–CT and 66% for LE2–CT, mirroring the good performance of TDDFT/FCD and underestimation of TDDFT/GMH that we observed for Np–TCNE. Δ-ALMO(MSDFT) overestimates the coupling between the LE1 and CT states with an error of 41%, while our scheme, Δ-ALMO(MSDFT2), again gives excellent agreement with the benchmark result, with only a relative error of 4%, even though the broken-symmetry ΔSCF treatment markedly underestimates the energy of the LE1 state (see Table S5 of the supplementary material).

FIG. 8.

Diabatic couplings (in meV) between the two LE states on 1-naphthol and the CT state that corresponds to the transition from 1-naphthol’s HOMO to CHCl3’s LUMO evaluated with TDDFT and ALMO-based diabatization schemes. The reference values obtained from EOM-EE-CCSD/GMH calculations for these two couplings are shown in gray dashed lines.

FIG. 8.

Diabatic couplings (in meV) between the two LE states on 1-naphthol and the CT state that corresponds to the transition from 1-naphthol’s HOMO to CHCl3’s LUMO evaluated with TDDFT and ALMO-based diabatization schemes. The reference values obtained from EOM-EE-CCSD/GMH calculations for these two couplings are shown in gray dashed lines.

Close modal

As discussed in Sec. II C, in the weak-coupling regime, one has to adopt an alternative form of the coupling that we denote as MSDFT2-wc. In the case of the LE2–CT coupling in the NpOH–CHCl3 complex, the weak coupling arises from the destructive overlap between NpOH’s LUMO+1 and CHCl3’s LUMO and is signaled by the smaller |Hab| than |Sab(Ea+Eb)/2|. Using the MSDFT-wc coupling prescription given in Eq. (19) to evaluate the LE2–CT coupling leads to a value of 14 meV, whose error is within ∼20% of the benchmark value (17 meV). Since Δ-ALMO(MSDFT) is less prone to issues in the weak-coupling regime, the LE2–CT diabatic coupling can be generated using the original MSDFT coupling scheme in Eq. (13), but the resulting LE2–CT coupling is much smaller than the reference value, with a relative error of 73%. These results suggest that the MSDFT2-wc approach that we proposed as a KS-DFT analogue to the generalized Slater–Condon rule for HF determinants whose interstate orbital overlap has one vanishingly small singular value in Sec. II C provides a practical and fairly accurate approach in the weak-coupling regime arising from the destructive orbital overlap in the NpOH–CHCl3 complex.

Understanding the photophysics associated with the pentacene dimer plays a central role in mechanistic studies of singlet fission,77,78,83,84 where a singlet exciton is split into two triplet excitons and which provides a promising route to efficient solar energy conversion.85–87 The initial and final states of singlet fission in the pentacene dimer can be described using fragment-based diabats, which are referred to as the LE and multiexcitonic (ME) states, respectively. Since the direct coupling between the LE and ME states has been shown to be small, it has been suggested that the transition from LE to ME state is mediated by a CT state that is more strongly coupled to both the LE and ME states via a superexchange mechanism.42,77,78,88 Here, we investigate the performance of our Δ-ALMO(MSDFT2) approach in capturing the coupling between the LE and CT diabats in the pentacene dimer as a function of intermolecular separation and monomer rotation angles. The change between the electronic configurations corresponding to these two diabats is illustrated in Fig. 9. The orientation of the two pentacene monomers (denoted as X and Y) is that used in a recent study77 and is based on the intermolecular configuration in crystalline pentacene. The LE state (|X*Y⟩) corresponds to the HOMO → LUMO transition on monomer X, while the CT state (|X+Y⟩) corresponds to the transition from the HOMO on X to the LUMO on Y.

FIG. 9.

The change in the electronic configuration for the diabatic couplings investigated for the pentacene dimer. This change is equivalent to an electron transfer from monomer X in an LE state to monomer Y, giving rise to a charge-separated pentacene dimer (the CT state). The LE state corresponds to the HOMO → LUMO transition on monomer X and the CT state corresponds to the transition from X’s HOMO to Y’s LUMO.

FIG. 9.

The change in the electronic configuration for the diabatic couplings investigated for the pentacene dimer. This change is equivalent to an electron transfer from monomer X in an LE state to monomer Y, giving rise to a charge-separated pentacene dimer (the CT state). The LE state corresponds to the HOMO → LUMO transition on monomer X and the CT state corresponds to the transition from X’s HOMO to Y’s LUMO.

Close modal

To benchmark our results, we compare them to the recent results by Zeng et al.77 who used a fourfold diabatization89,90 of adiabatic states obtained from extended multiconfigurational quasi-degenerate perturbation theory (XMCQDPT) calculations with an active space consisting of 4 electrons in 4 orbitals (4e, 4o). In particular, the two diabats we consider here are denoted as |eg⟩ (excited X with ground-state Y) and |ca⟩ (cationic X with anionic Y) in Ref. 77.

Figure 10 shows the LE–CT couplings obtained from Δ-ALMO(MSDFT) and Δ-ALMO(MSDFT2) as functions of the inter-monomer distance (dCC) and the rotation angles of monomers X and Y (ϕX and ϕY) in comparison to the reference values. For the dCC and ϕX scans, shown in Figs. 10(a) and 10(b), our Δ-ALMO(MSDFT2) approach gives excellent agreement with the XMCQDPT results with most of the relative errors below 10% and the largest being 15% at ϕX = −30°. Using our approach with the MSDFT coupling [Δ-ALMO(MSDFT)], on the other hand, leads to overestimation of the LE–CT coupling at all distances and angles in these scans, with the largest relative errors above 40%. For the scan of ϕY [Fig. 10(c)], Δ-ALMO(MSDFT2) shows reasonably good agreement with XMCQDPT at negative angles but markedly underestimates the coupling for ϕY > 0. Excluding the data point at ϕY = −30° where the magnitude of the XMCQDPT result is very close to zero (5 meV), the largest (unsigned) relative error of Δ-ALMO(MSDFT2) is 18% at ϕY = 30°. Nevertheless, Δ-ALMO(MSDFT2) still qualitatively captures the change in the LE–CT coupling strength in the ϕY > 0 regime and correctly predicts the value of ϕY (15°) corresponding to the maximum LE–CT coupling. Again, Δ-ALMO(MSDFT) systematically overestimates the LE–CT coupling in the ϕY scan, whose relative error is above 30% when ϕY < 0 while reducing to 10%–20% when ϕY > 0. Finally, we note that at ϕY = −30°, we invoked the MSDFT2-wc coupling scheme given in Eq. (19) to generate the Δ-ALMO(MSDFT2) result, since the condition |Hab|<|Sab(Ea+Eb)/2| is satisfied at this geometry. As shown in Fig. 10, this data point connects smoothly with the ones that were obtained from using the standard MSDFT2 coupling scheme.

FIG. 10.

Distance and angular dependence of the LE–CT coupling (meV) for photoinduced ET in the pentacene dimer obtained from ALMO-ΔSCF and TDDFT-based calculations. The geometric configuration of the dimer and the scanned coordinates are shown in the inset of the middle panel. (a) HLE-CT as a function of inter-fragment distance (measured by the top central C⋯ C distance, dCC, marked in the inset). [(b) and (c)] HLE-CT as functions of the rotation angle of monomers X and Y, respectively. The LE state corresponds to the HOMO → LUMO transition of monomer X and the CT state transition from the HOMO of monomer X to LUMO of monomer Y. The benchmark XMCQDPT(4e, 4o) results are from Ref. 77, for which the fourfold diabatization scheme89,90 was employed to construct the diabatic Hamiltonian and extract the diabatic couplings.

FIG. 10.

Distance and angular dependence of the LE–CT coupling (meV) for photoinduced ET in the pentacene dimer obtained from ALMO-ΔSCF and TDDFT-based calculations. The geometric configuration of the dimer and the scanned coordinates are shown in the inset of the middle panel. (a) HLE-CT as a function of inter-fragment distance (measured by the top central C⋯ C distance, dCC, marked in the inset). [(b) and (c)] HLE-CT as functions of the rotation angle of monomers X and Y, respectively. The LE state corresponds to the HOMO → LUMO transition of monomer X and the CT state transition from the HOMO of monomer X to LUMO of monomer Y. The benchmark XMCQDPT(4e, 4o) results are from Ref. 77, for which the fourfold diabatization scheme89,90 was employed to construct the diabatic Hamiltonian and extract the diabatic couplings.

Close modal

We also performed FCD diabatization of adiabatic excited states obtained from TDDFT calculations. We note that the standard GMH diabatization scheme cannot be used in this case due to the lack of a uniform charge-transfer direction. For this system, FCD diabatization requires the inclusion of at least the four lowest TDDFT excited states since the monomer HOMOs and LUMOs strongly couple with each other in the adiabatic picture, giving rise to four frontier orbitals of the full system (from HOMO−1 to LUMO +1). The details of the four-state FCD diabatization procedure are provided in Sec. S6 of the supplementary material. As shown in Fig. 10, TDDFT/FCD gives poorer agreement with the reference results than Δ-ALMO(MSDFT2) but, in general, performs similarly to Δ-ALMO(MSDFT) since it also overestimates the coupling for most of the scanned distances and rotation angles. In particular, TDDFT/FCD incorrectly predicts the turnover in the magnitude of coupling as a function of ϕY giving an early maximum at ϕY = 10° and a too rapid decay as the angle is further increased.

Here, we have introduced the Δ-ALMO(MSDFT2) method, which combines our ALMO-ΔSCF scheme to generate excited state diabats with our MSDFT2 scheme to calculate their couplings. We have shown that our method gives excellent agreement with the diabatic couplings obtained using high-level wavefunction-based schemes for a wide variety of systems relevant to DNA damage, charge separation in donor–acceptor dyads, photoinduced electron transfer to the solvent environment, and singlet fission. In addition, our ALMO approach is more generally applicable and yields better accuracy than other DFT-based methods, such as GMH and FCD diabatization of TDDFT adiabatic states, particularly for charged systems. Unlike schemes based on the adiabatic-to-diabatic transformation, Δ-ALMO(MSDFT2) gives access to variationally optimized diabatic states and thus provides easy access to nuclear gradients (i.e., forces) necessary for quantum dynamics simulations of photoinduced processes. While having a computational cost only ∼2–3 times that of the corresponding ground-state DFT calculation, our Δ-ALMO(MSDFT2) method yields LE–CT diabatic couplings typically within 10%–20% of benchmark values obtained from high-level wavefunction methods, such as GMH diabatization of EOM-CCSD adiabatic states. Importantly, although we used broken-symmetry ΔSCF to construct the excited-state diabats, our framework is flexible and therefore compatible with constructions that explicitly account for spin symmetry adaptations. This development thus opens the door to quantum dynamics simulations of photoinduced electron and hole transfer processes using diabats constructed directly from DFT calculations in the condensed phase.

See the supplementary material for procedures of GMH with three adiabatic states and FCD with a general number of adiabatic states as well as their applications to the indole–guanine complex and pentacene dimer systems, results of TDDFT and EOM-CCSD excited-state calculations required for GMH and FCD diabatization and other intermediate results, basis set sensitivity of diabatic couplings calculated with different schemes (PDF), geometries of the investigated complexes, and the original data for diabatic couplings (ZIP).

The authors greatly thank Professor Tao Zeng and Professor Nandini Ananth for providing the structures and diabatic coupling data for the pentacene dimer distance and angular scans and Professor Robert Cave for guidance on the implementation of the three-state GMH diabatization scheme. This material is based upon work supported by the National Science Foundation under Grant No. CHE-1652960. T.E.M. also acknowledges support from the Camille Dreyfus Teacher-Scholar Awards Program. This research also used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

The data that support the findings of this study are available within the article and its supplementary material.

1.
R. A.
Marcus
,
J. Chem. Phys.
24
,
966
(
1956
).
2.
N. S.
Hush
,
Trans. Faraday Soc.
57
,
557
(
1961
).
3.
R. A.
Marcus
,
Rev. Mod. Phys.
65
,
599
(
1993
).
4.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
).
5.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
6.
R. J.
Cave
and
M. D.
Newton
,
Chem. Phys. Lett.
249
,
15
(
1996
).
7.
R. J.
Cave
and
M. D.
Newton
,
J. Chem. Phys.
106
,
9213
(
1997
).
8.
A. A.
Voityuk
and
N.
Rösch
,
J. Chem. Phys.
117
,
5607
(
2002
).
9.
A. A.
Voityuk
,
J. Phys. Chem. C
117
,
2670
(
2013
).
10.
C.-P.
Hsu
,
Z.-Q.
You
, and
H.-C.
Chen
,
J. Phys. Chem. C
112
,
1204
(
2008
).
11.
C.-P.
Hsu
,
Acc. Chem. Res.
42
,
509
(
2009
).
12.
Z.-Q.
You
and
C.-P.
Hsu
,
J. Chem. Phys.
133
,
074105
(
2010
).
13.
J. E.
Subotnik
,
S.
Yeganeh
,
R. J.
Cave
, and
M. A.
Ratner
,
J. Chem. Phys.
129
,
244101
(
2008
).
14.
J. E.
Subotnik
,
R. J.
Cave
,
R. P.
Steele
, and
N.
Shenvi
,
J. Chem. Phys.
130
,
234102
(
2009
).
15.
J. E.
Subotnik
,
J.
Vura-Weis
,
A. J.
Sodt
, and
M. A.
Ratner
,
J. Phys. Chem. A
114
,
8665
(
2010
).
16.
J. E.
Subotnik
,
E. C.
Alguire
,
Q.
Ou
,
B. R.
Landry
, and
S.
Fatehi
,
Acc. Chem. Res.
48
,
1340
(
2015
).
17.
E.
Runge
and
E. K. U.
Gross
,
Phys. Rev. Lett.
52
,
997
(
1984
).
18.
A.
Dreuw
and
M.
Head-Gordon
,
Chem. Rev.
105
,
4009
(
2005
).
19.
M. E.
Casida
and
M.
Huix-Rotllant
,
Annu. Rev. Phys. Chem.
63
,
287
(
2012
).
20.
A.
Dreuw
and
M.
Head-Gordon
,
J. Am. Chem. Soc.
126
,
4007
(
2004
).
21.
Y.
Tawada
,
T.
Tsuneda
,
S.
Yanagisawa
,
T.
Yanai
, and
K.
Hirao
,
J. Chem. Phys.
120
,
8425
(
2004
).
22.
T.
Yanai
,
D. P.
Tew
, and
N. C.
Handy
,
Chem. Phys. Lett.
393
,
51
(
2004
).
23.
J.-D.
Chai
and
M.
Head-Gordon
,
J. Chem. Phys.
128
,
084106
(
2008
).
24.
M. A.
Rohrdanz
,
K. M.
Martins
, and
J. M.
Herbert
,
J. Chem. Phys.
130
,
054112
(
2009
).
25.
T.
Stein
,
L.
Kronik
, and
R.
Baer
,
J. Am. Chem. Soc.
131
,
2818
(
2009
).
26.
R.
Baer
,
E.
Livshits
, and
U.
Salzner
,
Annu. Rev. Phys. Chem.
61
,
85
(
2010
).
27.
I.
Kondov
,
M.
Čížek
,
C.
Benesch
,
H.
Wang
, and
M.
Thoss
,
J. Phys. Chem. C
111
,
11970
(
2007
).
28.
K.
Senthilkumar
,
F. C.
Grozema
,
F. M.
Bickelhaupt
, and
L. D. A.
Siebbeles
,
J. Chem. Phys.
119
,
9809
(
2003
).
29.
H.
Oberhofer
and
J.
Blumberger
,
Phys. Chem. Chem. Phys.
14
,
13846
(
2012
).
30.
C.
Schober
,
K.
Reuter
, and
H.
Oberhofer
,
J. Chem. Phys.
144
,
054103
(
2016
).
31.
Q.
Wu
and
T.
Van Voorhis
,
J. Chem. Phys.
125
,
164105
(
2006
).
32.
T.
Van Voorhis
,
T.
Kowalczyk
,
B.
Kaduk
,
L.-P.
Wang
,
C.-L.
Cheng
, and
Q.
Wu
,
Annu. Rev. Phys. Chem.
61
,
149
(
2010
).
33.
M.
Pavanello
,
T.
Van Voorhis
,
L.
Visscher
, and
J.
Neugebauer
,
J. Chem. Phys.
138
,
054101
(
2013
).
34.
P.
Ramos
and
M.
Pavanello
,
J. Chem. Theory Comput.
10
,
2546
(
2014
).
35.
A.
Cembran
,
L.
Song
,
Y.
Mo
, and
J.
Gao
,
J. Chem. Theory Comput.
5
,
2702
(
2009
).
36.
H.
Ren
,
M. R.
Provorse
,
P.
Bao
,
Z.
Qu
, and
J.
Gao
,
J. Phys. Chem. Lett.
7
,
2286
(
2016
).
37.
Y.
Mao
,
A.
Montoya-Castillo
, and
T. E.
Markland
,
J. Chem. Phys.
151
,
164114
(
2019
).
38.
S.
Difley
and
T.
Van Voorhis
,
J. Chem. Theory Comput.
7
,
594
(
2011
).
39.
Q.
Wu
and
T.
Van Voorhis
,
Phys. Rev. A
72
,
024502
(
2005
).
40.
Q.
Wu
and
T.
Van Voorhis
,
J. Chem. Theory Comput.
2
,
765
(
2006
).
41.
Q.
Wu
,
C.-L.
Cheng
, and
T.
Van Voorhis
,
J. Chem. Phys.
127
,
164119
(
2007
).
42.
W.-L.
Chan
,
T. C.
Berkelbach
,
M. R.
Provorse
,
N. R.
Monahan
,
J. R.
Tritsch
,
M. S.
Hybertsen
,
D. R.
Reichman
,
J.
Gao
, and
X.-Y.
Zhu
,
Acc. Chem. Res.
46
,
1321
(
2013
).
43.
D. J.
Thouless
,
Nucl. Phys.
21
,
225
(
1960
).
44.
Y.
Mo
,
L.
Song
, and
Y.
Lin
,
J. Phys. Chem. A
111
,
8291
(
2007
).
45.
T.
Ziegler
,
A.
Rauk
, and
E. J.
Baerends
,
Theor. Chim. Acta
43
,
261
(
1977
).
46.
T.
Kowalczyk
,
S. R.
Yost
, and
T. V.
Voorhis
,
J. Chem. Phys.
134
,
054128
(
2011
).
47.
G. M. J.
Barca
,
A. T. B.
Gilbert
, and
P. M. W.
Gill
,
J. Chem. Theory Comput.
14
,
1501
(
2018
).
48.
R. Z.
Khaliullin
,
M.
Head-Gordon
, and
A. T.
Bell
,
J. Chem. Phys.
124
,
204105
(
2006
).
49.
J. F.
Stanton
and
J.
Gauss
,
J. Chem. Phys.
101
,
8938
(
1994
).
50.
A. I.
Krylov
,
Annu. Rev. Phys. Chem.
59
,
433
(
2008
).
51.
R. S.
Mulliken
,
J. Chem. Phys.
23
,
1833
(
1955
).
52.
R. Z.
Khaliullin
,
E. A.
Cobar
,
R. C.
Lochan
,
A. T.
Bell
, and
M.
Head-Gordon
,
J. Phys. Chem. A
111
,
8753
(
2007
).
53.
P. R.
Horn
,
Y.
Mao
, and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
18
,
23067
(
2016
).
54.
A. T. B.
Gilbert
,
N. A.
Besley
, and
P. M. W.
Gill
,
J. Phys. Chem. A
112
,
13164
(
2008
).
55.
P.
Pulay
,
J. Comput. Chem.
3
,
556
(
1982
).
56.
H.
Stoll
,
G.
Wagenblast
, and
H.
Preubeta
,
Theor. Chim. Acta
57
,
169
(
1980
).
57.
P. R.
Horn
,
E. J.
Sundstrom
,
T. A.
Baker
, and
M.
Head-Gordon
,
J. Chem. Phys.
138
,
134119
(
2013
).
58.
P. O.
Löwdin
,
J. Chem. Phys.
18
,
365
(
1950
).
59.
A.
Amos
and
G.
Hall
,
Proc. R. Soc. A
263
,
483
(
1961
).
60.
A. J. W.
Thom
and
M.
Head-Gordon
,
J. Chem. Phys.
131
,
124113
(
2009
).
61.
A. F.
Morrison
,
Z.-Q.
You
, and
J. M.
Herbert
,
J. Chem. Theory Comput.
10
,
5366
(
2014
).
62.
A. F.
Morrison
and
J. M.
Herbert
,
J. Phys. Chem. Lett.
6
,
4390
(
2015
).
63.
T.
Fujita
and
Y.
Mochizuki
,
J. Phys. Chem. A
122
,
3886
(
2018
).
64.
P.
Ramos
,
M.
Papadakis
, and
M.
Pavanello
,
J. Phys. Chem. B
119
,
7541
(
2015
).
65.
Y.
Shao
,
Z.
Gan
,
E.
Epifanovsky
,
A. T. B.
Gilbert
,
M.
Wormit
,
J.
Kussmann
,
A. W.
Lange
,
A.
Behn
,
J.
Deng
,
X.
Feng
,
D.
Ghosh
,
M.
Goldey
,
P. R.
Horn
,
L. D.
Jacobson
,
I.
Kaliman
,
R. Z.
Khaliullin
,
T.
Kuś
,
A.
Landau
,
J.
Liu
,
E. I.
Proynov
,
Y. M.
Rhee
,
R. M.
Richard
,
M. A.
Rohrdanz
,
R. P.
Steele
,
E. J.
Sundstrom
,
H. L.
Woodcock
,
P. M.
Zimmerman
,
D.
Zuev
,
B.
Albrecht
,
E.
Alguire
,
B.
Austin
,
G. J. O.
Beran
,
Y. A.
Bernard
,
E.
Berquist
,
K.
Brandhorst
,
K. B.
Bravaya
,
S. T.
Brown
,
D.
Casanova
,
C.-M.
Chang
,
Y.
Chen
,
S. H.
Chien
,
K. D.
Closser
,
D. L.
Crittenden
,
M.
Diedenhofen
,
R. A.
DiStasio
,
H.
Do
,
A. D.
Dutoi
,
R. G.
Edgar
,
S.
Fatehi
,
L.
Fusti-Molnar
,
A.
Ghysels
,
A.
Golubeva-Zadorozhnaya
,
J.
Gomes
,
M. W. D.
Hanson-Heine
,
P. H. P.
Harbach
,
A. W.
Hauser
,
E. G.
Hohenstein
,
Z. C.
Holden
,
T.-C.
Jagau
,
H.
Ji
,
B.
Kaduk
,
K.
Khistyaev
,
J.
Kim
,
J.
Kim
,
R. A.
King
,
P.
Klunzinger
,
D.
Kosenkov
,
T.
Kowalczyk
,
C. M.
Krauter
,
K. U.
Lao
,
A. D.
Laurent
,
K. V.
Lawler
,
S. V.
Levchenko
,
C. Y.
Lin
,
F.
Liu
,
E.
Livshits
,
R. C.
Lochan
,
A.
Luenser
,
P.
Manohar
,
S. F.
Manzer
,
S.-P.
Mao
,
N.
Mardirossian
,
A. V.
Marenich
,
S. A.
Maurer
,
N. J.
Mayhall
,
E.
Neuscamman
,
C. M.
Oana
,
R.
Olivares-Amaya
,
D. P.
O’Neill
,
J. A.
Parkhill
,
T. M.
Perrine
,
R.
Peverati
,
A.
Prociuk
,
D. R.
Rehn
,
E.
Rosta
,
N. J.
Russ
,
S. M.
Sharada
,
S.
Sharma
,
D. W.
Small
,
A.
Sodt
,
T.
Stein
,
D.
Stück
,
Y.-C.
Su
,
A. J. W.
Thom
,
T.
Tsuchimochi
,
V.
Vanovschi
,
L.
Vogt
,
O.
Vydrov
,
T.
Wang
,
M. A.
Watson
,
J.
Wenzel
,
A.
White
,
C. F.
Williams
,
J.
Yang
,
S.
Yeganeh
,
S. R.
Yost
,
Z.-Q.
You
,
I. Y.
Zhang
,
X.
Zhang
,
Y.
Zhao
,
B. R.
Brooks
,
G. K. L.
Chan
,
D. M.
Chipman
,
C. J.
Cramer
,
W. A.
Goddard
,
M. S.
Gordon
,
W. J.
Hehre
,
A.
Klamt
,
H. F.
Schaefer
,
M. W.
Schmidt
,
C. D.
Sherrill
,
D. G.
Truhlar
,
A.
Warshel
,
X.
Xu
,
A.
Aspuru-Guzik
,
R.
Baer
,
A. T.
Bell
,
N. A.
Besley
,
J.-D.
Chai
,
A.
Dreuw
,
B. D.
Dunietz
,
T. R.
Furlani
,
S. R.
Gwaltney
,
C.-P.
Hsu
,
Y.
Jung
,
J.
Kong
,
D. S.
Lambrecht
,
W.
Liang
,
C.
Ochsenfeld
,
V. A.
Rassolov
,
L. V.
Slipchenko
,
J. E.
Subotnik
,
T.
Van Voorhis
,
J. M.
Herbert
,
A. I.
Krylov
,
P. M. W.
Gill
, and
M.
Head-Gordon
,
Mol. Phys.
113
,
184
(
2015
).
66.
J.-D.
Chai
and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
10
,
6615
(
2008
).
67.
W. J.
Hehre
,
R.
Ditchfield
, and
J. A.
Pople
,
J. Chem. Phys.
56
,
2257
(
1972
).
68.
M. J.
Frisch
,
J. A.
Pople
, and
J. S.
Binkley
,
J. Chem. Phys.
80
,
3265
(
1984
).
69.
A. D.
Becke
,
J. Chem. Phys.
88
,
2547
(
1988
).
70.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
71.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
72.
A. D.
Becke
,
J. Chem. Phys.
98
,
5648
(
1993
).
73.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
74.
S.
Hirata
and
M.
Head-Gordon
,
Chem. Phys. Lett.
314
,
291
(
1999
).
75.
J. F.
Stanton
and
R. J.
Bartlett
,
J. Chem. Phys.
98
,
7029
(
1993
).
76.
S.
Chaudhuri
,
A.
Acharya
,
E. T. J.
Nibbering
, and
V. S.
Batista
,
J. Phys. Chem. Lett.
10
,
2657
(
2019
).
77.
T.
Zeng
,
R.
Hoffmann
, and
N.
Ananth
,
J. Am. Chem. Soc.
136
,
5755
(
2014
).
78.
T. C.
Berkelbach
,
M. S.
Hybertsen
, and
D. R.
Reichman
,
J. Chem. Phys.
138
,
114103
(
2013
).
79.
C.
Butchosa
,
S.
Simon
,
L.
Blancafort
, and
A.
Voityuk
,
J. Phys. Chem. B
116
,
7815
(
2012
).
80.
M.
Rust
,
J.
Lappe
, and
R. J.
Cave
,
J. Phys. Chem. A
106
,
3930
(
2002
).
81.
M.
Head-Gordon
,
A. M.
Grana
,
D.
Maurice
, and
C. A.
White
,
J. Phys. Chem.
99
,
14261
(
1995
).
82.
J. R.
Platt
,
J. Chem. Phys.
17
,
484
(
1949
).
83.
P. M.
Zimmerman
,
F.
Bell
,
D.
Casanova
, and
M.
Head-Gordon
,
J. Am. Chem. Soc.
133
,
19944
(
2011
).
84.
X.
Feng
,
A. V.
Luzanov
, and
A. I.
Krylov
,
J. Phys. Chem. Lett.
4
,
3845
(
2013
).
85.
M. B.
Smith
and
J.
Michl
,
Chem. Rev.
110
,
6891
(
2010
).
86.
M. B.
Smith
and
J.
Michl
,
Annu. Rev. Phys. Chem.
64
,
361
(
2013
).
87.
D. N.
Congreve
,
J.
Lee
,
N. J.
Thompson
,
E.
Hontz
,
S. R.
Yost
,
P. D.
Reusswig
,
M. E.
Bahlke
,
S.
Reineke
,
T.
Van Voorhis
, and
M. A.
Baldo
,
Science
340
,
334
(
2013
).
88.
T. C.
Berkelbach
,
M. S.
Hybertsen
, and
D. R.
Reichman
,
J. Chem. Phys.
141
,
074705
(
2014
).
89.
H.
Nakamura
and
D. G.
Truhlar
,
J. Chem. Phys.
115
,
10353
(
2001
).
90.
H.
Nakamura
and
D. G.
Truhlar
,
J. Chem. Phys.
117
,
5576
(
2002
).

Supplementary Material