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 $\Delta $-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.

## I. INTRODUCTION

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 CO_{2} 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 process^{6–12} or states that are maximally localized on the donor or acceptor moieties^{13–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 states^{18,20} that are essential to photoinduced ET processes. Although these issues can somewhat be alleviated by range-separated hybrid (RSH) functionals^{21–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 method^{38} 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 work^{42} used this procedure followed by a Thouless transformation^{38,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 ΔSCF^{45,46} to target molecular excited states while maintaining a balanced accuracy for different types of electronic excited states^{46,47} as well as the ability of ALMOs^{37,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 diabatization^{6,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.

## II. METHODS

### A. LE and CT states from ALMO calculations

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 shown^{37} 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 *D*⋯*A*, where the donor (*D*) has *n*_{D} electrons and the acceptor (*A*) has *n*_{A} electrons. The ground state of this system, |*DA*⟩, can be represented as

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, *E*^{KS}[**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 (**C**_{o}) through

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 approach^{45,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

which we denote as |*D*^{*}*A*⟩_{0} 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 *i*th occupied orbital to *a*th virtual orbital.

To variationally relax the ALMOs in the LE state, we then solve the locally projected SCF (Stoll) equation^{48,56} for each fragment,

Here, *X* ∈ {*D*, *A*} is the fragment index, **F** ≡ **F**[**P**] is the AO Fock matrix constructed from the 1PDM (**P**), **P**_{X} is the projector for fragment *X*, and *ϵ*_{X} and **C**_{X} 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 **C**_{X} 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

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 *n*_{D} of the resulting matrix, and when *X* = *A*, one takes columns *n*_{D} + 1 to *n*_{D} + *n*_{A}. The matrices $[Co(\sigma oo\u22121)]X$ and **C**_{o,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(\sigma oo\u22121)]X$ and **C**_{o,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 procedure^{47} 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)$,

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 *n*_{D} orbitals that possess the largest projection onto the span of the initial set of occupied orbitals, $PD*j$,

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,

where the superscript (LE) indicates that the ALMOs are optimized within the LE state in contrast to the unrelaxed fragment orbitals in |*D*^{*}*A*⟩_{0}. 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

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.

### B. Diabatic coupling between ALMO LE and CT states

To evaluate the electronic coupling between LE and CT diabats, $\u27e8D*A|\u0124|D+A\u2212\u27e9$, one can utilize the original multistate DFT (MSDFT) approach^{35,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,

To extract the electronic coupling, $Hab=\u27e8D*A|\u0124|D+A\u2212\u27e9$, between these two ALMO diabats, one transforms **H**′ into an orthogonal basis. Using Löwdin’s symmetric orthogonalization scheme,^{58} one obtains $H=S\u22121/2H\u2032S\u22121/2$, where the interstate overlap matrix is $[S]xy\u2261\u27e8\psi x|\psi y\u27e9$ and *x*, *y* ∈ {*a*, *b*}. In the two-state case, this yields

where $Sab$ is given by

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)\u2261(Co(a))TSCo(b)$ to denote the overlap matrix between occupied orbitals from diabats |*ψ*_{a}⟩ and |*ψ*_{b}⟩.

In both MSDFT^{35,36} and MSDFT2^{37} approaches, the diagonal elements of **H**′ are the KS energies of each diabat, i.e., $Haa\u2032=EaKS[P(a)]$ and $Hbb\u2032=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\u2032$. In the original MSDFT scheme, $Hab\u2032$ is given by a KS-DFT correction to the Hartree–Fock (HF) interstate coupling^{59,60} result,

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

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,

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\u0303ab$),^{37}

where $P\u0303ab=(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 shown^{37} 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.

### C. MSDFT2 in the weak-coupling regime

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\u2032$|). 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\u2032|<|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 rule^{58,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\xafo(a)$ and $C\xafo(b)$. We then perform a singular value decomposition (SVD) on the interstate occupied orbital overlap matrix, $S\xafoo(ab)$,

and then use the **U** and **V** matrices obtained to transform the orthogonal orbitals $C\xafo(a)$ and $C\xafo(b)$ to the Löwdin-paired orbitals,^{58} $C\u0303o(a)=C\xafo(a)U$ and $C\u0303o(b)=C\xafo(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\u0303o(a)$ and $C\u0303o(b)$) corresponds to the pair of orbitals that overlap the least (associated with the smallest singular value, *s*_{1}), denoted here as $C\u0303o,1(a)$ and $C\u0303o,1(b)$. The remaining *N* − 1 vectors in each set, $C\u0303o,N\u22121(a)$ and $C\u0303o,N\u22121(b)$, correspond to the remaining occupied orbitals (where *N* denotes the total number of electrons) and **s**_{N−1} corresponds to the set of *N* − 1 singular values. This results in two re-ordered sets of Löwdin-paired orbitals: $C\u0303o(a)=[C\u0303o,1(a),C\u0303o,N\u22121(a)]$ and $C\u0303o(b)=[C\u0303o,1(b),C\u0303o,N\u22121(b)]$.

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

where

are two transition density matrix-like objects, $P\u0303N\u22121(ab)=(PN\u22121(ab)+PN\u22121(ba))/2$ is the symmetrized version of $PN\u22121(ab)$, and $S\u0303ab=\u220fi\u22601Nsi$ is the reduced interstate overlap. The KS Fock matrix (**F**_{KS}) 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 (**V**_{xc}), 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}

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\u0303o,1(a)$ and $C\u0303o,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.

## III. IMPLEMENTATION AND COMPUTATIONAL DETAILS

### A. Implementation of Δ-ALMO(MSDFT2)

We extended our previous implementation of ALMO(MSDFT2)^{37} in Q-Chem^{65} 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 method^{47} was employed in combination with the DIIS algorithm^{55} 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 scheme^{45} 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\u2032$), 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\u2032|<|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\u0303o(a)$ and $C\u0303o(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.

### B. Methods based on adiabatic-to-diabatic transformation

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

where *E*_{1} and *E*_{2} are the energies of the two adiabats, *μ*_{1} and *μ*_{2} are their dipole moments, and $\mu ^12$ denotes the projection of the transition dipole (*μ*_{12}) along the charge-transfer direction: $e^CT=(\mu 1\u2212\mu 2)/|\mu 1\u2212\mu 2|$.

The FCD method,^{8} on the other hand, defines the diabats as eigenstates of the fragment charge difference matrix (Δ*Q*_{ij} = *Q*_{ij}(*D*) − *Q*_{ij}(*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

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.

### C. Computational setup

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 functional^{66} and the 6-31+G(d) basis set^{67,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 (BLYP^{69,70} and PBE^{71}), global hybrid (B3LYP^{72} and PBE0^{73}), and two other range-separated hybrid (RSH) functionals (CAM-B3LYP^{22} and LRC-*ω*PBEh^{24}). 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-CCSD^{50,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–CHCl_{3} 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–CHCl_{3} 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 *d*_{CC} (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.

## IV. RESULTS

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–CHCl_{3} 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 fission^{42,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.

### A. Cationic indole–guanine complex

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 *H*_{GS-CT} and *H*_{LE-CT}. Since the diabatic states involved in the forward and reverse processes are identical, the diabatic couplings are the same.

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 diabatization^{80} 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 *H*_{GS-CT} = 277 meV and *H*_{LE-CT} = 60 meV.

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 densities^{81} 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 *H*_{GS-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.

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 work^{37} 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.

### B. Naphthalene–tetracyanoethylene complex

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 ^{1}L_{a} and ^{1}L_{b} 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.

The equilibrium geometry of this complex is of *C*_{2v} 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 B_{2} and B_{1} 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 B_{1} irreducible representation leading to a non-zero coupling with the CT state. Using EOM-EE-CCSD to compute the first two B_{2} 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

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 B_{1} 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 B_{1} (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.

### C. Naphthol–CHCl_{3} complex

The 1-naphthol–CHCl_{3} 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 CHCl_{3} 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–CHCl_{3} 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 CHCl_{3}, 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 CHCl_{3}’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 CHCl_{3} 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.

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

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–CHCl_{3} complex, the weak coupling arises from the destructive overlap between NpOH’s LUMO+1 and CHCl_{3}’s LUMO and is signaled by the smaller |$Hab\u2032$| 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–CHCl_{3} complex.

### D. Pentacene dimer

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 study^{77} 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.

To benchmark our results, we compare them to the recent results by Zeng *et al.*^{77} who used a fourfold diabatization^{89,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 (*d*_{CC}) and the rotation angles of monomers X and Y (*ϕ*_{X} and *ϕ*_{Y}) in comparison to the reference values. For the *d*_{CC} 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\u2032|<|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.

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.

## V. CONCLUSIONS

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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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