An effective potential for Frenkel excitons

Abstract

Excitation energy transfer (EET) is a ubiquitous process in life and materials sciences. Here, a new and computationally efficient method of evaluating the electronic EET couplings between interacting chromophores is introduced that is valid in a wide range of intermolecular distances. The proposed approach is based on the effective elimination of electron repulsion integrals from the excitonic Hamiltonian matrix elements via the density-fitting approach and distributed multipole approximation. The excitonic Hamiltonian represented in a basis including charge transfer (CT) states is re-cast in terms of the effective one-electron potential functions (EOPs) and adapted into the effective fragment parameter (EFP) framework. Calculations for model systems indicate that the speedup of at least three orders of magnitude, as compared to the state-of-the-art methods, can be achieved while maintaining the accuracy of the EET couplings even at short intermolecular distances.

---

1 Introduction

Excitation energy transfer (EET) is a nonadiabatic phenomenon governing the nonradiative transfer of excitation energy between two chromophores characterized by overlapping spectral features which enable their resonance coupling. It is of fundamental importance for understanding the early stages of photosynthesis or functioning of light-harvesting complexes in general. Thus EET is of central importance to many branches of materials and life sciences. Its strong dependence on the distance and orientation of the coupled chromophores led to the development of several time-resolved spectroscopic techniques, in particular the single-molecule Förster resonance energy transfer, allowing the determination of the mechanisms of biomolecular processes.^1,2

Since the EET phenomenon and the theoretical framework for its description were recently reviewed thoroughly,^3–8 here we provide only a brief outline of the recent developments. At large intermolecular distances EET coupling can be adequately described by the Förster model⁹ which assumes the purely Coulombic interaction between the transition densities of chromophores, and computationally efficient schemes have been developed over recent years that take into account the spatial extent and shape of transition densities.^10–12 However, an accurate description of EET coupling constants become more difficult at short intermolecular distances due to wavefunction overlap and exchange interactions.^13–15 A robust description of EET coupling in such cases is provided by the Frenkel exciton model¹⁶ in which the excitons are localized on individual chromophores. Such an approach is quite reliable in the weak-coupling limit¹⁷ and can be extended to account for intermolecular charge transfer states.^18–21

To accurately account for short-range effects, fragment density approaches (FDAs) were developed,^13,22,23 in which the analysis of the Hamiltonian is performed in the eigenstate basis for the initial and final EET states of interest. In particular, the fragment excitation difference (FED) approach,¹³ which is a generalization of the energy splitting in dimer (ESD) method²⁴ for nonsymmetric systems, provides reasonably accurate values of EET at all physically meaningful intermolecular distances. A different strategy was adopted within the so-called transfer-integral (TI) approach.^24,25 This technique stems from a standard two-state diabatic Landau–Zener model in which the off-diagonal Hamiltonian matrix element, called the transfer integral,²⁴ couples the initial and final states in a nonadiabatic process of electron or excitation energy transfer. This approach can be generalized to interactions among many-level systems, which introduces the so-called indirect coupling matrix elements via intermediate (bridge) states in addition to the direct interaction between initial and final states.^25,26 In the context of electronic EET this approach was proposed by Scholes and Harcourt as the superexchange coupling (SEC) method^26,27 and was recently reintroduced by Fujimoto as the transition-density-fragment interaction transfer-integral (TDFI-TI) model.^28–31

The main advantage of the TI approach is that the EET coupling can be directly expressed via the donor and acceptor transition densities by means of perturbation theory which allows time-consuming calculations in the eigenstate space of an entire system to be avoided. Unfortunately, from the perspective of ab initio electronic structure calculations all the above techniques require the evaluation of the four-center electron repulsion integrals (ERIs) and excited state properties for a given instantaneous orientation of interacting chromophores. Even though nonadiabatic dynamics simulations of interacting chromophores are nowadays possible due to growing computational power and introduction of efficient computational protocols,^18–20 the time-consuming contractions involving 4-index tensors are still a serious bottleneck which generally prohibits systematic studies of large aggregates.

In this paper we introduce a new efficient method to calculate the electronic EET couplings within the TI approach that is valid in both the long- and short-range regimes and is based on the recently introduced effective one-electron potential (EOP) operators technique to eliminate ERIs from fragment-based models.³² Although we focus here on the singlet–singlet (SS) EET, extension to the triplet–triplet (TT) EET can be readily achieved from the corresponding TI model.

To test the proposed technique we assume the configuration interaction singles³³ (CIS) wavefunction of interacting molecules. For clarity of demonstration, we consider a simple model case in which only the non-degenerate frontier orbitals of the coupled chromophores are kept active to construct the excitonic Hamiltonian.^26,27 An example of the thus formulated variant of the transfer-integral CIS (TI/CIS) method, proposed by Fujimoto,³¹ is a starting point for the introduction of the effective one-electron potential transfer-integral (EOP-TI) approach, which is the main purpose of this paper. Although in a general case more orbitals might be required to describe the bridging composite states in the SEC model, we believe that the application of the EOP method would then be straightforward. In the following sections we outline the relevant theoretical framework and discuss the results obtained for model systems, focusing on the accuracy and efficiency of the introduced approximations. Finally, we show that the proposed EOP-TI method de facto introduces a way to construct an effective fragment potential for EET couplings.

2 Theoretical methodology

2.1 Transfer-integral model of EET coupling

In the simplest version of the TI/CIS approach, the Hamiltonian of an aggregate of molecules A and B is constructed assuming the weak-coupling exciton model approximation involving the following four configurations:

|Φ₁〉 = |Ψ^(e)_A ⊗ Ψ^(g)_B〉, (1a)

|Φ₂〉 = |Ψ^(g)_A ⊗ Ψ^(e)_B〉, (1b)

|Φ₃〉 = |Ψ⁽⁺⁾_A ⊗ Ψ⁽⁻⁾_B〉, (1c)

|Φ₄〉 = |Ψ⁽⁻⁾_A ⊗ Ψ⁽⁺⁾_B〉, (1d)

where superscripts g and e denote the ground and excited states of a molecule, and ‘+’ and ‘−’ label the cationic and anionic states, respectively, whereas

|Ψ_X ⊗ Ψ_Y〉 ≡ [scr A, script letter A]{|Ψ_X〉 ⊗ |Ψ_Y〉} (2)

denotes the antisymmetrized Hartree product of the monomer wavefunctions with [scr A, script letter A] being the standard antisymmetrization operator and ‘⊗’ denoting the outer product. Instead of diagonalizing the excitonic Hamiltonian in such a basis, one may adopt the perturbation theory to obtain the approximate expression for the EET coupling constant:^26,31where E_n ≡ 〈Φ_n|[script letter H] − E₀|Φ_n〉 and E₀ is the ground state energy of the aggregate. The first two site energies are given by³¹

E₁ ≡ 〈Φ₁|[script letter H] − E₀|Φ₁〉 = E^A_e→g + ΔE^A_e→g(B), (4a)

E₂ ≡ 〈Φ₂|[script letter H] − E₀|Φ₂〉 = E^B_e→g + ΔE^B_e→g(A), (4b)

where E^X_e→g is the excitation energy of isolated molecule X and ΔE^X_e→g(Y) accounts for the mean-field environmental effect due to the other molecule, which is given bywith the effective one-electron potential given byIf basis functions 3 and 4 are approximated as the HOMO and LUMO of monomers (i.e., the frontier orbital approximation), the remaining two site energies are given by

E₃ ≡ 〈Φ₃|[script letter H] − E₀|Φ₃〉 ≈ −ε^A_H + ε^B_L − (ϕ^A_Hϕ^A_H|ϕ^B_Lϕ^B_L), (7a)

E₄ ≡ 〈Φ₄|[script letter H] − E₀|Φ₄〉 ≈ ε^A_L − ε^B_H − (ϕ^A_Lϕ^A_L|ϕ^B_Hϕ^B_H), (7b)

where ε^X_j is the energy of the jth Hartree–Fock orbital ϕ^X_j associated with molecule X, and the two-electron integral (ERI) is defined by

The first term in eqn (3) is the, so-called, direct coupling, which is composed of two contributions as follows:

〈Φ₁|[script letter H]|Φ₂〉 = V⁰_Coul + V⁰_Exch (9)

where the Förster (or Coulomb) coupling is given bywhereas the Dexter (or pure exchange) coupling isIn the above equations, P^X_g→e is the transition one-particle density matrix of molecule X in AO representation. V⁰_Coul, which is usually a leading term for SS EET (but vanishes for TT EET), can be efficiently evaluated by various multipole expansion or electrostatic potential fitting models published to date.^10–12V⁰_Exch was reported to be rather negligible for SS EET even for close intermolecular contacts,^12,31 but is known to be non-negligible for TT EET.^13,34,35 However, unlike the direct Coulombic coupling, the second- and third-order terms in eqn (3), which are collectively referred to by Fujimoto as the indirect coupling,³¹ as well as direct exchange coupling for TT EET, become increasingly important at short intermolecular distances. Unfortunately, the evaluation of such couplings is computationally quite demanding, even if the frontier orbital approximation of eqn (7) is assumed, since it requires evaluation of the ERIs.

Since the TI/CIS model was found to provide quite reliable estimates of EET coupling constants for many organic molecules,^15,31,36 we consider this model as a starting point to develop a computationally more efficient scheme that does not require time-consuming calculations of ERIs at all. First, in Section 2.2 we derive the EOP model based on the TI/CIS approach proposed earlier by Scholes and Ghiggino²⁶ and Fujimoto³¹ (hereonin referred to as the EOP-TI/CIS or, in short, EOP-TI model). Next, we validate it against the various formulations of the TI/CIS method with increasing accuracy for a few selected model complexes in Section 3. Finally, we conclude our findings in Section 4. The auxiliary derivations are outlined in the appendices.

2.2 EOP-based model

In this section, we derive the approximation of the TI/CIS model that is completely free from ERIs. This is possible due to effective elimination of ERIs according to the general prescription:³²where the capital italic letters denote subsets of orbitals associated with a particular fragment and X = A or B, [scr F, script letter F]_t is an arbitrary linear functional of ERIs, and ô^A_s and v̂^eff,A_i represent an arbitrary one-electron operator and an EOP operator associated with a molecule or fragment A, respectively. First of all, the coulombic EET coupling can be approximated by the distributed multipole expansion of the transition densities of interacting monomers in a dimer using the distributed transition-density cumulative atomic multipole moments (TrCAMM) approach, developed by some of us previously:¹²

V⁰_Coul ≈ ρ^A_e→g ⊙ ρ^B_g→e, (13)

where the symbol ‘⊙’ denotes the tensor contractions involved in the multipole expansion of the interaction energy between the generalized one-electron densities and the superscript ‘0’ indicates that the matrix elements are not affected by the overlap between molecular wavefunctions (this issue shall be addressed later in Section 2.2.4).

The diagonal elements for the first two basis states that are necessary for the overlapping corrections and also for the evaluation of the indirect EET couplings can be approximated (see the ESI†) by neglecting the environmental correction ΔE^X_e→g, i.e.,

E₁ ≈ E^A_e→g, (14a)

E₂ ≈ E^B_e→g. (14b)

Although this correction was found to be negligible in all the systems studied in this work (Tables S1 and S2, ESI†), this assumption should be tested for other chromophores and electronic transitions. All the other components of the model involve ERIs between the HOMO and LUMO of monomers A and B. These can be grouped into the following three categories:³²

1. Integrals of type (ϕ^X_iϕ^X_j|ϕ^Y_kϕ^Y_l), or Coulomb-like integrals. These involve HOMO–LUMO interactions in the expressions for E₃ and E₄ from eqn (7). These shall be treated in this work by the distributed multipole expansion of the form:

(ϕ^X_iϕ^X_j|ϕ^Y_kϕ^Y_l) ≈ ρ^X_ij ⊙ ρ^Y_kl, (15)

whereand P^X_ij is the associated effective one-particle density matrix in the AO basis given by

P^X_ij;αβ = C^X_αi(C^X_βj)*. (17)

In this work, the CAMM method is used to generate the distributed multipole moments from such density matrices, as discussed in ref. 32.

2. Integrals of type (ϕ^X_iϕ^X_j|ϕ^X_kϕ^Y_l), or overlap-like integrals. Matrix elements needed to evaluate EET coupling due to electron and hole transfer processes involve overlap-like ERIs. These shall be treated by the extended density fitting (EDF) technique,³²i.e.,

(ϕ^X_iϕ^X_j|ϕ^X_kϕ^Y_l) ≡ −(ϕ^X_k|v̂^X_ij|ϕ^Y_l) (18)

whereIn the above, the ‘DF’ label indicates the auxiliary AO basis set used for density fitting, and the effective potential operator is defined by³²with its spatial form given as(here, note that ϕ_j(r) ≡ (r|ϕ_j)).

3. Integrals of type (ϕ^X_iϕ^Y_j|ϕ^X_kϕ^Y_l), or exchange-like integrals. These cannot be treated by means of the effective one-electron operators, but nonetheless constitute the V⁰_Exch matrix element and also conribute to the indirect coupling. Such integrals can be treated by a number of approximations, which allow one to express ERIs in terms of the one-electron overlap and other ERIs.³⁷ In this work, the Mulliken approximation is assumed where only two-electron Coulomb-like integrals are necessary, i.e.,

Thus, one can find that the remaining diagonal matrix elements are

E₃ ≈ −ε^A_H + ε^B_L − ρ^A_HH ⊙ ρ^B_LL, (23a)

E₄ ≈ ε^A_L − ε^B_H − ρ^A_LL ⊙ ρ^B_HH. (23b)

In the following subsections, the off-diagonal Hamiltonian matrix elements shall be discussed in detail in order to effectively remove all ERIs from the formulation.

2.2.1 Indirect coupling via electron and hole transfer. The indirect coupling terms can be intuitively decomposed into electron, hole and charge transfer terms.^26,31 The electron transfer (ET) matrix elements arealong with the twin term, V⁰_ET2 ≡ 〈Φ₂|[script letter H]|Φ₄〉, which were introduced earlier by Fujimoto.³¹ In the above, [scr F, script letter F] is the Fock operator of the entire molecular aggregate. Note that here t^A_H→L are the CIS amplitudes of isolated monomer A, whereas reflects the approximate amplitude in a dimer with monomer B being in its ground state. Here, we also approximate the Fock operator by neglecting the interchromophore relaxation:where T̂ ≡ −½∇² is the one-electron kinetic energy operator, V̂^X_nuc is the operator associated with the electrostatic potential energy between nuclei and electrons in molecule X, and the Coulomb and exchange operators are defined byNote that the decomposition of the kinetic energy operator into particular molecular contributions is arbitrary. However, consider the following partitioning of the total Fock operator:

[scr F, script letter F]₀ ≈ [capital G, script]^A₀ + [capital G, script]^B₀, (27)

in which the auxiliary operator is defined asand is the Fock operator of isolated molecule X. Here, the contribution from the one-electron kinetic energy operator was symmetrically partitioned in between the two interacting molecules. Thus one can rewrite eqn (24) asNow, by using the EOP technique³² for overlap-like effective-potential operator matrix elements of type (AA|AB) one can collect operators associated to a particular molecule and group them into one effective one-electron operator asandIn the above transformations that involve the EDF in an auxiliary basis set space, V^A;ET_ξ;HL and V^B;ET_η;L are the EOP matrix elements associated with the above effective potentials, separately defined for molecules A and B. This allows one to further recast V⁰_ET1 intowhich has a particularly simple form as compared to the original expression in eqn (24) that involves four-center electron repulsion integrals. Note that, if the EOP matrices are considered as effective fragment parameters (EFPs), only overlap integrals between auxiliary basis functions and the LUMOs of molecules A and B are necessary, i.e.,whereand C^X_βU is the SCF LCAO-MO matrix element for isolated molecule X associated with the Uth MO. It is straightforward to show that the twin term ET2 isNote that EOP matrices for ET1 and ET2 terms share a similar structure. Therefore a joint superscript ‘ET’ is used here to label the EOP matrices and additional subscripts L and HL denote the dependence on only the LUMO or on both the HOMO and LUMO, respectively.

Similar considerations apply to the hole transfer (HT) matrix elements. The corresponding expression

can be recast aswith the effective potential matrices given byThe twin term HT2 is accordinglyAs it was shown, to compute the Hamiltonian matrix elements associated with the electron and hole transfer, in total eight different EOP matrices need to be precomputed and stored in a file. Note that each of these matrices is actually just a vector of length equal to the size of the auxiliary basis set chosen. Therefore, it can be predicted that the computational cost of the evaluation of the EOP-based Hamiltonian matrix elements is negligible as compared to the original formulation of the TI/CIS method. In Appendix A, the explicit and convenient to implement working formulae for the EOP matrix elements are given in terms of the AOs.

2.2.2 Indirect coupling via charge transfer. The Hamiltonian matrix element that is associated with the charge transfer process is

V⁰_CT = 2(ϕ^A_Hϕ^B_L|ϕ^A_Lϕ^B_H) − (ϕ^A_Hϕ^B_H|ϕ^A_Lϕ^B_L). (40)

Since it involves ERIs of type (AB|AB), EOPs cannot be readily defined here. However, the application of the Mulliken approximation from eqn (22) to eqn (40) yieldswhere

r^X_HL ≡ (ϕ^X_Hϕ^X_H|ϕ^X_Lϕ^X_L). (42)

First of all, note that r^A_HL and r^B_HL are just constant numbers associated with the isolated monomers and independent of the molecular aggregate. Also, the remaining four ERIs can be considered as Coulombic interactions between the HOMO and LUMO of unperturbed molecular wavefunctions:

(ϕ^A_Hϕ^A_H|ϕ^B_Hϕ^B_H) ≈ ρ^A_HH ⊙ ρ^B_HH, (43a)

(ϕ^A_Lϕ^A_L|ϕ^B_Lϕ^B_L) ≈ ρ^A_LL ⊙ ρ^B_LL, (43b)

(ϕ^A_Hϕ^A_H|ϕ^B_Lϕ^B_L) ≈ ρ^A_HH ⊙ ρ^B_LL, (43c)

(ϕ^A_Lϕ^A_L|ϕ^B_Hϕ^B_H) ≈ ρ^A_LL ⊙ ρ^B_HH. (43d)

In the limiting case when multipole expansion is truncated on the monopole terms, one finds thatwhere the charge centroids of the HOMOs are given by

r^X_H = (ϕ^X_H|r̂|ϕ^X_H). (45)

Note that, once occupied molecular orbitals are localized, the above approximation should yield the correct description. However, the LUMOs are usually quite delocalized. Here, we use the distributed cumulative atomic charges to represent these:³⁸and one can find that assuming this one obtainsThe two remaining ERIs can then be approximated as follows:and similarly for the twin ERIs. To summarize the derivation, the CT matrix element is approximately given byThus, the computational cost of the above expression is negligible as compared to the original formula in eqn (40).

2.2.3 Approximate pure exchange coupling. The Mulliken approximation can also be used to treat the pure exchange coupling (eqn (11)) which leads toHere, such ERIs that involve AOs in two different molecules are approximated bywhere r_μσ is the distance between the centers of the μth and σth AOs. That leads to the following approximate expression:with Q_μν ≡ (μμ|νν). Note that Q^A_μν and Q^B_λσ are just matrix elements of constant matrices. Although these are not second-rank tensors (i.e., they do not transform under Euclidean rotations as other AO matrices), they can be easily computed from a small set of proper fourth-rank tensors, which is shown in Appendix B. Therefore, the EFP methodology can still be used to calculate pure exchange Dexter contribution to the EET coupling and the computational cost of the overall procedure is negligible as compared to the evaluation of eqn (11) involving ERIs in the AO space directly.

2.2.4 Complete model of EET coupling. To summarize, the overall form of the EOP-TI/CIS model is still given as in the original TI/CIS model^26,31 by

V ≈ V_Direct + V_Indirect, (53)

where the direct and indirect coupling constants are, respectively,

V_Direct ≡ V_Coul + V_Exch + V_Ovlp, (54a)

V_Indirect ≡ V^TI(2) + V^TI(3), (54b)

withIn the above, the direct EET coupling contributions arewhereas the indirect EET coupling is evaluated from the overlap-corrected ET, HT and CT matrix elements that readwhere ‘t’ denotes one of the matrix element types (ET1, ET2, HT1, HT2, CT) and S_t is the appropriate overlap integral between basis states.

Despite the same general form, the differences between the original TI/CIS model and the newly derived EOP-based approach are quite significant. The Hamiltonian matrix elements are evaluated from the effective fragment parameters, for which only the one-electron integrals are needed. Therefore, the computational cost is significantly reduced as compared to the parent TI/CIS model. In particular, the V⁰_Coul term is evaluated by using the TrCAMM method developed previously,¹² whereas V⁰_Exch and V⁰_CT are computed using the approximate Mulliken formulae in eqn (52) and (49), respectively. In turn, the electron and hole transfer terms V⁰_ET1, V⁰_ET2, V⁰_HT1 and V⁰_HT2 are treated by the EDF with EOP effective matrices as effective parameters spanned in the auxiliary AO basis set (eqn (32), (35), (37), (39) and (58)–(63)). Note that the evaluation of eqn (56) and (57) is computationally inexpensive because it involves computing only approximate overlap integrals between basis states.³¹

On the other hand, the relaxation effects due to the intermolecular interactions are not included in the current EOP-based model. Note that in the TDFI-TI/CIS model the basis states and all the one-particle density matrices are built from the self-consistently adjusted fragment ground-state densities according to the DFI method.²⁸ However, as demonstrated in the following section, we do not expect major differences due to the relaxation effects.

2.3 Adaptation of the EOP-TI/CIS method as effective fragment potentials

In the previous sections, all interfragment ERIs were effectively eliminated from the original TI/CIS equations by using the EOP technique. Here, note that the following quantities, which are sufficient to evaluate the resulting EOP-TI/CIS EET coupling constant without calculating the interfragment ERIs, are functions of the isolated electronic chromophores:

• TrCAMMs of monomers A and B,

• excitation energies E^A_g→e and E^B_g→e, as well as the corresponding transition density matrices P^A_e→g and P^B_g→e,

• (pp|qq)^A and (pp|qq)^B tensor elements from Appendix B,

• V^A;ET_ξ;HL, V^B;ET_η;L, V^A;HT_ξ;L and V^B;HT_η;HL vector elements, each of length equal to the size of the auxiliary EOP basis set,

• CIS HOMO–LUMO amplitudes t^A_H→L and t^B_H→L,

• CAMMs of HOMO and LUMO densities from eqn (43), as well as the canonical orbital energies,

• distributed atomic charges q^A_x,L and q^B_y,L from eqn (46),

• localized HOMO centroids r^A_H and r^B_H from eqn (45), as well as intrafragment (and invariant with respect to the Euclidean rotations of chromophores) Coulomb-type ERIs r^A_HL and r^B_HL from eqn (42), and finally.

• primary and auxiliary basis sets.

Therefore, only one preparative calculation per chromophore is required to obtain all of the above quantities that constitute the EFPs, much like obtaining the EFPs for calculating the interaction energy by means of the second generation of the effective fragment potential (EFP2) method.³⁹ The rate-determining step in computing the EFPs is the calculation of the CIS amplitudes which can be nowadays performed efficiently even for large molecules.

Such EOP-TI/CIS EFPs can be subsequently used in routine calculations of EET coupling constants in conjunction with the molecular dynamics and rigid molecule algorithm as developed previously in the context of calculations of the solvatochromic vibrational frequency shifts of infrared probes in condensed phases (SolEFP).^40–45 Note that, once EFPs are computed, the evaluation of the EOP-TI/CIS equations requires only one-electron overlap integrals which are relatively inexpensive to calculate.

3 Results and discussion

All the models and approximations that were discussed so far were implemented in our inhouse plugin EOPDev⁴⁶ to the Psi4 quantum chemistry package,⁴⁷ along with the direct CIS method utilizing the Davidson–Liu simultaneous expansion approach.^48–50 In this work, the tested levels of the TI/CIS theory are characterized by the treatment of the Fock operator in the expressions for the off-diagonal ET and HT Hamiltonian matrix elements. For example, in the case of V⁽⁰⁾_ET1, the approximations denoted here as V^TI([scr F, script letter F]) and V^TI([scr F, script letter F]₀) collectively refer to using eqn (24) and (29), respectively. V^TI([scr F, script letter F]^DFI₀) indicates that the approximate Fock matrix was used as in eqn (29) but including the interfragment interactions in a self-consistent manner through the DFI process.²⁸ The EOP-based model is denoted as V^TI([scr F, script letter F]^EOP₀) and additionally utilizes the TrCAMM approximation for V⁽⁰⁾_Coul, as well as the Mulliken approximation for V⁽⁰⁾_Exch (eqn (52)) and V⁽⁰⁾_CT (two variants are tested: the CAMM-based vatiant from eqn (41) with eqn (43) or the distributed charge-based variant from eqn (49)). Environmental correction terms from eqn (5) were included in all models except the V^TI([scr F, script letter F]^EOP₀) model. In order to test the EOP-TI/CIS model we chose the EET coupling between the lowest-lying bright ¹ππ* states of ethylene (1 ¹B_1u) and 7-aminocoumarin (1 ¹A′) in stacked molecular dimers, having D_2h and C_s symmetry, respectively. To test the performance and applicability of the proposed methodology we also performed preliminary calculations on the “special pair” in the reaction center of photosystem I (PS-I).^10,14,51 To improve the convergence of the distributed multipole expansion of the Coulombic EET coupling constants, the transition densities used to compute the TrCAMMs were symmetrized as suggested by Belov et al.⁵² Occupied core orbitals were included in the CIS calculations except from the PS-I model for which the frozen-core approximation was utilized.

All the reported results were obtained assuming the 6-31G(d) basis set⁵³ and, when required, the cc-pVDZ-JKFIT density-fitting basis set.⁵⁴ To obtain EOP matrices for ET and HT Hamiltonian matrix elements, the EDF scheme³² was utilized with aug-cc-pVDZ-JKFIT⁵⁴ as the auxiliary EOP basis set. Further details of calculations and the selected geometries of the studied complexes are provided in the ESI.†

3.1 Accuracy of transfer integral approximations

The total magnitudes of EET couplings between stacked ethylene molecules in the 1 ¹B_1u states and 7-aminocoumarin in the 1 ¹A′ states, obtained for the selected geometries, are compiled in Table 1. The equilibrium geometry of the ethylene dimer was optimized using the MP2/6-31G(d) method assuming D_2h symmetry of the complex and EET couplings were then computed for the selected intermolecular distances (R_AB) by translation of the molecules along the intermolecular axis. The R_AB distance between the neighboring carbon atoms of 4.169 Å corresponds to the minimum energy structure. At a distance of 6.0 Å the V⁰_Exch term drops below 0.002 cm⁻¹ and the total coupling is virtually identical to Coulomb contribution, whereas the structure for R_AB = 3.0 Å represents the region of large overlap of molecular wavefunctions where V⁰_Exch exceeds 1700 cm⁻¹. The 7-aminocoumarin (7AC) dimer was prepared as reported in our previous study.¹² The structures of two flat 7AC molecules optimized at the MP2/6-31g(d) level were superimposed and shifted with respect to the intermolecular axis. Such a rigid potential energy (PE) scan shows a minimum at R_AB = 3.6 Å, while the regions of large and negligible overlap, identified analogously as in the ethylene dimer (i.e. based on the magnitudes of V⁰_Exch), occur at intermolecular distances of roughly 2.6 and 4.6 Å, respectively.

Table 1 EET couplings (in cm⁻¹) computed between the lowest-lying bright ¹ππ* states in the studied ethylene and 7-aminocoumarin stacked dimers at various levels of approximation (indicated in parentheses by the corresponding Fock matrix estimate) to the CIS reference values. All the results were computed for the selected intermolecular distances R_AB (given in Å) assuming the 6-31G(d) basis set

R AB	(Ethylene)2			(7AC)2
	3.0	4.169	6.0	2.6	3.6	4.6
V ^ESD	9893	1973	495	12 291	2433	868
V ^TI([scr F, script letter F])	7990	1744	495	14 576	1782	840
V ^TI([scr F, script letter F]^DFI0)	10 156	1762	495	21 241	1942	835
V ^TI([scr F, script letter F]0)	9953	1766	495	21 135	1949	842
V ^TI([scr F, script letter F]^EOP0)	10 481	1772	494	19 617	2175	837

---

The reference values of the total EET coupling V^ESD were calculated using the energy splitting in dimer (ESD) method,^24,55 in which the coupling is calculated as half of the splitting between the corresponding adiabatic excited state energies in a dimer. In the case of the studied symmetric complexes this method provides the accurate electronic coupling. The transfer-integral V^TI([scr F, script letter F]) values, computed assuming the exact Fock matrix of the dimer, are very close to the ESD reference at large and intermediate intermolecular separations and become more approximate at shorter distances. The V^TI([scr F, script letter F]) results represent the limiting accuracy of the TI/CIS frontier-orbitals model which is semi-quantitative at equilibrium geometries and becomes only qualitatively correct at shorter distances. Their calculations for the ethylene dimer took 9.2 s of CPU time on one core of the 1.2 GHz AMD EPYC 7301 16-core processor. The V^TI([scr F, script letter F]^DFI₀) values correspond to the original transition-density-fragment interaction method (TDFI-TI) of Fujimoto,³¹ because the density-fragment interaction (DFI) method^28,29 is used to take into account molecular surroundings in the Fock matrix of the chromophores in a self-consistent fashion. These are overestimated at short distances and underestimated at the equilibrium geometry by roughly 10–20% (depending on the system). As expected, all the tested models predict virtually identical results at large intermolecular separations.

The effect of mutual polarization taken into account in the TDFI-TI method via the block-localized self-consistent density-fragment interaction procedure²⁸ seems to be rather small. This may be inferred from Table 1 by a comparison of V^TI([scr F, script letter F]^DFI₀) and V^TI([scr F, script letter F]₀) magnitudes, the relative differences of which are below 2% even at short intermolecular distances. In the latter approach the dimer Fock matrix was constructed from the Fock matrices of monomers assuming the symmetric partitioning introduced in eqn (27) and (28), and thus the interacting monomer ground-state densities were not relaxed in this case. Both calculations took about 3.2 s in the same CPU (cf. Table S3 in the ESI†) which shows about 20% speedup with respect to the calculations of V^TI([scr F, script letter F]) using the exact Fock matrix of the ethylene dimer. Note, however, that the reported timings do not include the time to solve the SCF and CIS problems for the monomers – V^TI([scr F, script letter F]^DFI₀) calculations are obviously longer than V^TI([scr F, script letter F]₀) ones in terms of the elapsed time since the Hartree–Fock problem for the monomers is solved in a block-localized self-consistent fashion for the whole complex.

It is interesting to note that the corresponding EOP-TI calculations took only a few milliseconds of CPU time, which is nearly three orders of magnitude faster (cf. Table S3, ESI†), while showing virtually no signs of loss in the accuracy, as the V^TI([scr F, script letter F]^EOP₀) values obtained for the ethylene dimer reproduce consistently the V^TI([scr F, script letter F]₀) estimates with a 1–2% accuracy even at short distances. In the case of the 7AC complex the speedup is more spectacular, since V^TI([scr F, script letter F]₀) calculations took about 25 minutes while V^TI([scr F, script letter F]^EOP₀) only 79 milliseconds (nearly 20 000 times faster). The speedup is even larger for PS-I, for which the respective timings are about 137 hours vs. 6 seconds. The relative differences of total couplings are larger in this system (within about ±7%) but one could argue that at short distances the EOP-TI values are, in fact, consistently closer to the reference values. This, however, is entirely due to the cancellation of errors. The fact is that essentially all the discussed TI models show quite similar accuracy, but the EOP-TI approach is simply much more efficient.

The same conclusions arise from the analysis of the PE cuts shown in Fig. 1, where in panels (a) and (c) the total TI([scr F, script letter F])/CIS couplings and their components are plotted along with the reference ESD values for the model systems in the vicinity of their equilibrium geometries. The total TI([scr F, script letter F]) coupling (marked with gray lines and asterisks) underestimates slightly the ESD values (shown as black lines with circles) at equilibrium and shorter distances. The agreement is at least semi-quantitative unless the region of very large overlap is considered, which would however be hardly accessible in any real physical situation (note the corresponding interaction energies plotted as gray dashed lines vs. the right axis). For instance, in the 7AC dimer this occurs for ΔR_AB below −0.9 Å (cf.Fig. 1c), where the TI([scr F, script letter F]) estimates start to diverge from the ESD coupling, but the intermolecular repulsion in this region exceeds 190 kJ mol⁻¹.

Fig. 1 EET couplings and their partitioning in symmetric ethylene (a and b) and 7AC dimers (c and d) as a function of intermolecular separation. ΔR_AB is the translation with respect to the equilibrium distance. In the left panels (a and c) the reference ESD values are plotted along with the best TI/CIS estimates obtained for the exact dimeric Fock matrix. In the right panels (b and d) the TI([scr F, script letter F]₀) results plotted as solid lines are compared with their EOP-TI([scr F, script letter F]₀) estimates shown as points.

3.2 Partitioning of EET couplings

In Tables 2 and 3 the detailed partitioning of EET couplings in ethylene and 7AC dimers, respectively, is shown for the selected intermolecular separations. At large distances the total coupling V^TI is determined solely by the direct electrostatic contribution which can be accurately approximated from TrCAMM expansion. At the equilibrium geometry the exchange, overlap and third-order V^TI(3) terms are still negligible. However, the total indirect coupling is already significant and almost exclusively due to the second-order terms V^TI(2). The same applies also to the special pair of photosystem I (cf. Table S4 of the ESI†), which was chosen to test the performance of the EOP-TI scheme. The vanishingly small values of V_Ovlp terms indicate that the relevant adiabatic states of monomers in the studied model systems are virtually orthogonal. Only in the case of the ethylene dimer in the large-overlap region do they rise slightly. Even though the studied ethylene dimer is very weakly bound, having an interaction energy of less than 1 kJ mol⁻¹ in the MP2/6-31G(d) minimum energy structure, the corresponding interaction energy in the 7AC dimer amounts to about 20 kJ mol⁻¹. Thus, it appears that the above conclusions could be general for interacting systems in the vicinity of equilibrium geometry. In the region of large overlap the exchange and third-order terms become significant. In fact, the indirect coupling becomes the leading component with a dominant contribution from the second-order terms. In general the V^TI([scr F, script letter F]₀) partitioning, being a starting point for the introduction of EOP approximation, agrees quantitatively with the corresponding TDFI estimates. The differences between these two methods are essentially due to the lack of relaxation effects in the former. The largest differences are in the magnitudes of second-order terms, which still fall within a few percent margin, unless the Mulliken approximation is utilized. It is interesting to note that in the case of the special pair of PS-I (Table S4, ESI†) at the experimental geometry the agreement is virtually quantitative, for both the total coupling (171 vs. 178 cm⁻¹) and its components. The latter being mostly due to the direct electrostatic coupling V_Coul (157 vs. 161 cm⁻¹) and the second-order indirect contribution V^TI(2) (15.8 vs. 17.6 cm⁻¹).

Table 2 Partitioning of EET couplings computed between the lowest-lying bright ¹ππ* states of ethylene molecules in stacked alignment at various levels of approximation to the TI/CIS method. Indentation indicates constituents of a given term or its approximation. All the results, given in cm⁻¹, were computed for the selected intermolecular distances R_AB (given in Å) assuming the 6-31G(d) basis set. The superscript ‘M’ indicates the values estimated using the Mulliken approximation: eqn (52) for the EOP-TI model and eqn (50) for the remaining models. The V^TI(3) values in the EOP-TI results were computed assuming either the distributed multipolar expansion truncated on |R|⁻⁵ terms or on atomic monopoles (values in parentheses) in V_CT calculations

TI/CIS model	TI([scr F, script letter F])	TDFI	TI([scr F, script letter F]0)	EOP	TI([scr F, script letter F])	TDFI	TI([scr F, script letter F]0)	EOP
	R AB = 3.0				4.169
V ^TI	7990	10 156	9953	10 481 (10 500)	1744	1762	1766	1772 (1772)
VDirect	3239	3221	3239	4093	1626	1622	1626	1622
VCoul	4896	4885	4896	—	1654	1649	1654	—
V^TrCAMMCoul	5133	5120	5133	5133	1638	1634	1638	1638
VExch	−1743	−1751	−1743	—	−30	−29	−30	—
V^MExch	−1174	−1179	−1174	−1125	−19	−18	−19	−18
VOvlp	86	86	86	86	2	2	2	2
VIndirect	4751	6936	6714	6388 (6407)	118	140	141	150 (150)
V^TI(2)	5189	7708	7462	7533	118	140	141	150
V^TI(3)	−438	−773	−748	—	0	0	0	—
V^TI(3),M	−695	−1225	−1185	−1145 (−1126)	0	0	0	0 (0)

---

Table 3 Partitioning of EET couplings computed between the lowest-lying bright ¹ππ* states of two 7AC molecules in stacked alignment at various levels of approximation to the TI/CIS method. Indentation indicates constituents of a given term or its approximation. All the results, given in cm⁻¹, were computed for the selected intermolecular distances R_AB (given in Å) assuming the 6-31G(d) basis set. The superscript ‘M’ indicates the values estimated using the Mulliken approximation: eqn (52) for the EOP-TI model and eqn (50) for the remaining models. The V^TI(3) values in the EOP-TI results were computed assuming either the distributed multipolar expansion truncated on |R|⁻⁵ terms or on atomic monopoles (values in parentheses) in V_CT calculations

TI/CIS model	TI([scr F, script letter F])	TDFI	TI([scr F, script letter F]0)	EOP	TI([scr F, script letter F])	TDFI	TI([scr F, script letter F]0)	EOP
	R AB = 2.6				3.6
V ^TI	14 576	21 241	21 135	19 703 (19 617)	1782	1942	1949	2175 (2175)
VDirect	859	700	859	1987	1294	1286	1294	1352
VCoul	2803	2903	2803	—	1397	1389	1397	—
V^TrCAMMCoul	3344	3543	3344	3344	1422	1414	1422	1422
VExch	−1960	−2220	−1960	—	−104	−104	−104	—
V^MExch	−1412	−1605	−1412	−1373	−73	−73	−73	−71
VOvlp	16	18	16	15	1	1	1	1
VIndirect	13 717	20 541	20 276	17 716 (17 630)	487	656	655	823 (823)
V^TI(2)	16 427	26 271	24 776	21 545	490	659	658	827
V^TI(3)	−2710	−5730	−4500	—	−2	−3	−3	—
V^TI(3),M	−2993	−6292	−4971	−3829 (−3914)	−2	−3	−3	−4 (−4)

---

The excellent agreement of EOP-TI estimates with TI([scr F, script letter F]₀) results is demonstrated for a broader range of intermolecular distances in Fig. 1. In panels (b) and (d) the latter are plotted as solid lines while the former as points. In the vicinity of equilibrium geometry the agreement is nearly quantitative and starts deteriorating in the hardly accessible region of very large overlap. However, since the EOP-TI couplings start to underestimate the TI([scr F, script letter F]₀) values, which in turn overestimate the reference ESD values, the EOP-TI results are consistently closer to the reference. In fact, this is true also at larger intermolecular separations, since the major difference in EOP-TI components is found for the second-order terms being overestimated at large distances and underestimated at short distances.

The basis set extension effects are discussed thoroughly in the ESI† (see Fig. S1). Generally, the choice of aug-cc-pVDZ-JKFIT as the auxiliary EOP basis set is sufficient even for more flexible orbital basis sets, at least for the first-row atom molecules studied here. The choice of the orbital basis set is more critical, particularly for the calculations of indirect couplings. However, the latter constitute only about 20–30% of the total TI/CIS coupling at minimum energy geometries of the studied model systems. The direct contribution, on the other hand, is sufficiently saturated already in the split-valence polarized double-zeta basis set. Interestingly, the TI/CIS values converge quickly to the exact ESD coupling with basis set extension.

3.3 Accuracy of Mulliken and multipolar expansion approximations

The accuracy of the Mulliken approximation for the estimation of pure exchange coupling V_Exch appears to be semi-quantitative at the equilibrium geometry and becomes qualitative at shorter intermolecular distances, as the V^M_Exch term underestimates significantly the magnitude of exchange coupling (cf.Tables 2 and 3). It is interesting to note here that the simplified model of V_Exch in the EOP-TI model from eqn (52) works qualitatively well for both studied systems, as the coupling strength is systematically underestimated. In the case of indirect coupling terms, the Mulliken approximation can be utilized to estimate the magnitude of V_CT matrix elements and, hence, the third-order V^TI(3) terms. In fact, this approximation is necessary in the EOP-TI model. However, it works very well and provides the quantitative accuracy of V^M_CT terms at equilibrium geometry but becomes more approximate in the region of large overlap of molecular wavefunctions (cf.Table 4). It is also interesting to note that the magnitudes of V^M_CT terms, calculated in the EOP-TI method assuming the distributed multipolar expansion truncated either on |R|⁻⁵ terms or on atomic monopoles (values in parentheses), are virtually identical at equilibrium geometry and remain very close even at shorter distances.

Table 4 Electron, hole and charge transfer (ET, HT and CT, respectively) Hamiltonian matrix elements calculated at various levels of approximation for EET couplings between the lowest-lying bright ¹ππ* states in ethylene and 7AC dimers in stacked alignment. All the results, given in cm⁻¹, were computed for the selected intermolecular distances R_AB (given in Å) assuming the 6-31G(d) basis set. The superscript ‘M’ indicates the values estimated using the Mulliken approximation: eqn (52) for the EOP-TI model and eqn (50) for the remaining models. The V_CT values in the EOP-TI results were computed assuming the distributed multipolar expansion truncated either on |R|⁻⁵ terms or on atomic monopoles (values in parentheses)

TI/CIS model	TI([scr F, script letter F])	TDFI	TI([scr F, script letter F]0)	EOP	TI([scr F, script letter F])	TDFI	TI([scr F, script letter F]0)	EOP
(Ethylene)2	R AB = 3.0				4.169
V ET1	−4306	−4485	−4393	−4516	−1170	−1172	−1172	−1248
V HT1	6626	9422	9337	9591	1161	1379	1383	1383
V CT	−849	−853	−849	—	−15	−15	−15	—
V^MCT	−1346	−1352	−1346	−1347 (−1325)	−23	−22	−23	−23 (−22)

---

TI/CIS model	TI([scr F, script letter F])	TDFI	TI([scr F, script letter F]0)	EOP	TI([scr F, script letter F])	TDFI	TI([scr F, script letter F]0)	EOP
(7AC)2	R AB = 2.6				3.6
V ET1	−7025	−7348	−7403	−7317	−1782	−1849	−1838	−2107
V HT1	8368	12 195	11 978	11 976	1928	2510	2511	2767
V CT	−1163	−1316	−1163	—	−65	−65	−65	—
V^MCT	−1284	−1445	−1284	−1286 (−1315)	−67	−67	−67	−67 (−68)

---

In Table 4 the relevant Hamiltonian matrix elements necessary for calculating indirect coupling terms are reported for the studied model systems at selected geometries and levels of TI approximation. The magnitudes of V_ET1 and V_ET2 terms are identical due to the symmetry of the studied complexes and the same applies to the corresponding HT terms. It should be noted that the signs of these terms depend on the arbitrary phases of frontier orbitals (this is true for all the TI components). However, this does not affect the total magnitude of EET coupling. In this work we adhere to certain convention for which the total coupling is positive. Unfortunately, we have found certain discrepancies between our results and those reported by Fujimoto.³¹ A careful analysis of the formula for the CT term given in eqn (40) shows that it should always be negative. A similar inspection of the formulas for ET and HT terms (eqn (24) and (36)) shows that, even though their signs may vary, the resultant V^TI(2) and V^TI(3) terms, defined in eqn (55), should have opposite signs. Fujimoto, however, reported both terms as having the same positive sign. We were also unable to reproduce the total magnitudes of V^TI([scr F, script letter F]^DFI₀), which, as shown in Fig. 3 of Fujimoto's paper,³¹ seem to be virtually identical to V^ESD. This is quite surprising considering the crudeness of approximation (frontier orbitals and block-localized Fock matrix). For instance our TDFI-TI results, reported in Table 1, show a relative error of about 3% at a distance of 3 Å.

4 Summary and a few concluding remarks

In this work, the effective one-electron potential operators technique of electron-repulsion integral elimination was used to introduce a new efficient variant of the transfer integral approach to calculate excitation energy transfer couplings. The EOP technique was adopted to estimate the EET couplings between singlet bright states in ethylene and 7-aminocoumarine dimers assuming the TI/CIS super-exchange model of Scholes and Harcourt within the frontier orbital approximation introduced earlier by Fujimoto.

It was shown that the EOP-TI/CIS method provides semi-quantitative results even at short distances and the computational cost is reduced by at least 3 orders of magnitude, as compared to the TDFI-TI/CIS calculations. This speedup is expected to increase with the system size as in the 7-aminocoumarine dimer the EOP-TI/CIS calculations were nearly 20 000 times faster than the most efficient variant of the TDFI-TI/CIS approach.

In fact, the EOP-TI/CIS method was found to be so efficient that it is readily applicable in molecular dynamics simulations. This is due to its EOP-based formulation which defines the EET coupling in terms of effective fragment potential parameters which could be determined prior to the simulation in a fully automated fashion for a given excited state of each of the interacting molecules. We note here that the computational costs can be further reduced by applying much smaller tailored density fitting auxiliary EOP basis sets.³²

It is anticipated here that the EOP-TI/CIS method can be extended to a general multi-state model with multiple electron–hole pairs and bridging CT states,^19,20 because the mathematical form of ERIs within the TIs remains the same as that in the current EOP-TI/CIS model. In other words the extended density fitting EOP method is in principle applicable to any ET- and HT-like matrix elements, whereas the Mulliken approximation and distributed multipole expansion are applicable to any CT-like matrix element.

Conflicts of interest

There are no conflicts to declare.

Appendices

A Explicit formulae for ET and HT EOP matrix elements

As developed in ref. 32, the matrix elements of the EOP operators are computed from the EDF scheme, according to which (in general the TI/CIS case)whereIn the above equations, the auxiliary matrices are defined by

T^X;M,N_mX = S_mm⁻¹S_maT^X;M,N_aXS_XX⁻¹, (60a)

S_XX⁻¹ = ({T^X;M,N_aX}^†S_amS_mm⁻¹S_maT^X;M,N_aX)^½, (60b)

T^X;M,N_aX = S_aa⁻¹[scr Q, script letter Q]U^X;M,N_aX, (60c)

where the overlap matrix block corresponding to two AO basis sets ‘x’ and ‘y’ isand the symbol ‘m’ denotes the (generally incomplete) auxiliary EOP AO basis set, whereas the symbol ‘a’ denotes the intermediate AO basis set that approximately fulfills the resolution of identity. The similarity transformation matrix T^X;M,N_aX is obtained from the eigenvectors of the co-variance matrix, W^X;M,N_aa:The operator [scr Q, script letter Q] in eqn (60c) selects only eigenvectors (columns of) U^X;M,N_aX associated with the non-vanishing eigenvalues stored in the diagonal matrix g^X;M,N_XX. The elements of vector f^X;M_N for M = ET1, ET2, HT1 and HT2 and N = H, L and HL are given byHere, the auxiliary one-electron operator [capital G, script]^X₀ represented in space spanned by the intermediate and primary AO basis sets needs to be computed. This can be done by inverting the following matrix equation to obtain G^X_αβ:whereandIn the above, Z_x is the atomic number of the xth atom and . In this work, the auxiliary EOP basis set was assumed to be identical to the intermediate basis set, i.e., conventional density fitting was used with [R^X;M,N_ma]_ξε = δ_ξε. However, by utilizing the general EDF scheme from eqn (58), additional reduction of computational cost is possible through optimization of the auxiliary EOP basis set.

B Effective parameters for Dexter coupling

The Euclidean rotation of a shell of n_p atomic orbitals with angular momentum p is given by

|p′) = R^(p)|p), (66)

in which the shell |p) ≡ {…,φ_μ(r),…} and R^(p) is the n_p × n_p rotation matrix. Consider now the integrals of type (pp|qq). It immediately follows that

(p′p′|q′q′) = R^(p) ⊗ R^(p) ⊗ R^(q) ⊗ R^(q)::(pp|qq) (67)

and the Q matrix from eqn (52) can be readily constructed by using

Q ≡ (diag[pp]|diag[qq]), (68)

where ‘diag’ denotes taking the diagonal of a matrix. The symbol ‘::’ denotes the fourfold tensor contraction for arbitrary matrices a, b, c, d and a fourth-rank tensor B. Note that the number of (pp|qq) integrals scales quadratically with the size of the basis set which means that memory requirements to store them on a disk or in a file are relatively small and comparable to sizes of typical AO matrices. Therefore, all the (pp|qq) integrals can be treated as effective fragment parameters and the computational cost of superimposition onto a target geometry, which is described by eqn (68) and (67), is relatively small.

Acknowledgements

This project was funded by the National Science Centre, Poland (grant no. 2016/23/P/ST4/01720) within the POLONEZ 3 fellowship. The POLONEZ programme has received funding from the Horizon 2020 research and innovation programme of the European Union under the Marie Skłodowska-Curie grant agreement no. 665778.

Notes and references

1. E. Haas, Protein Folding Handbook, Wiley-VCH Verlag GmbH, Weinheim, Germany, 2005, pp. 573–633.
2. E. Lerner, T. Cordes, A. Ingargiola, Y. Alhadid, S. Chung, X. Michalet and S. Weiss, Science, 2018, 359, 6373.
3. T. Renger, Photosynth. Res., 2009, 102, 471–485.
4. Z.-Q. You and C.-P. Hsu, Int. J. Quantum Chem., 2014, 114, 102–115.
5. A. Olaya-Castro and G. D. Scholes, Int. Rev. Phys. Chem., 2011, 30, 49–77.
6. A. Chenu and G. D. Scholes, Annu. Rev. Phys. Chem., 2015, 66, 69–96.
7. C. Curutchet and B. Mennucci, Chem. Rev., 2017, 117, 294–343.
8. L. Cupellini, M. Corbella, B. Mennucci and C. Curutchet, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2019, 9, e1392.
9. T. Förster, Ann. Phys., 1948, 437, 55–75.
10. M. E. Madjet, A. Abdurahman and T. Renger, J. Phys. Chem. B, 2006, 110, 17268–17281.
11. K. J. Fujimoto, J. Chem. Phys., 2014, 141, 214105.
12. B. Błasiak, M. Maj, M. Cho and R. W. Góra, J. Chem. Theory Comput., 2015, 11, 3259–3266.
13. C.-P. Hsu, Z.-Q. You and H.-C. Chen, J. Phys. Chem. C, 2008, 112, 1204–1212.
14. M. E.-A. Madjet, F. Müh and T. Renger, J. Phys. Chem. B, 2009, 113, 12603–12614.
15. B. Shi, F. Gao and W. Liang, Chem. Phys., 2012, 394, 56–63.
16. J. Frenkel, Phys. Rev., 1931, 37, 1276–1294.
17. A. Sisto, D. R. Glowacki and T. J. Martinez, Acc. Chem. Res., 2014, 47, 2857–2866.
18. T. Fujita, S. Atahan-Evrenk, N. P. D. Sawaya and A. Aspuru-Guzik, J. Phys. Chem. Lett., 2016, 7, 1374–1380.
19. X. Li, R. M. Parrish, F. Liu, S. I. L. Kokkila Schumacher and T. J. Martínez, J. Chem. Theory Comput., 2017, 13, 3493–3504.
20. T. Fujita and Y. Mochizuki, J. Phys. Chem. A, 2018, 122, 3886–3898.
21. I. O. Glebov, M. I. Kozlov and V. V. Poddubnyy, Comput. Theor. Chem., 2019, 1153, 12–18.
22. R. J. Cave and M. D. Newton, Chem. Phys. Lett., 1996, 249, 15–19.
23. A. A. Voityuk and N. Rösch, J. Chem. Phys., 2002, 117, 5607–5616.
24. M. D. Newton, Chem. Rev., 1991, 91, 767–792.
25. M. A. Ratner, J. Phys. Chem., 1990, 94, 4877–4883.
26. G. D. Scholes and K. P. Ghiggino, J. Chem. Phys., 1994, 101, 1251–1261.
27. G. D. Scholes and R. D. Harcourt, J. Chem. Phys., 1996, 104, 5054–5061.
28. K. Fujimoto and W. Yang, J. Chem. Phys., 2008, 129, 054102.
29. K. J. Fujimoto and S. Hayashi, J. Am. Chem. Soc., 2009, 131, 14152–14153.
30. K. J. Fujimoto, J. Chem. Phys., 2010, 133, 124101.
31. K. J. Fujimoto, J. Chem. Phys., 2012, 137, 034101.
32. B. Błasiak, J. D. Bednarska, M. Chołuj, R. W. Góra and W. Bartkowiak, J. Comput. Chem., 2020 DOI:10.1002/jcc.26462.
33. J. B. Foresman, M. Head-Gordon, J. A. Pople and M. J. Frisch, J. Phys. Chem., 1992, 96, 135–149.
34. Z.-Q. You, C.-P. Hsu and G. R. Fleming, J. Chem. Phys., 2006, 124, 044506.
35. S. Yeganeh and T. V. Voorhis, J. Phys. Chem. C, 2010, 114, 20756–20763.
36. S. S. Skourtis, C. Liu, P. Antoniou, A. M. Virshup and D. N. Beratan, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 8115–8120.
37. O. Ayed, E. Bernard and B. Silvi, THEOCHEM, 1986, 135, 159–168.
38. W. Sokalski and R. Poirier, Chem. Phys. Lett., 1983, 98, 86–92.
39. M. S. Gordon, Q. A. Smith, P. Xu and L. V. Slipchenko, Annu. Rev. Phys. Chem., 2013, 64, 553–578.
40. B. Błasiak, H. Lee and M. Cho, J. Chem. Phys., 2013, 139, 044111.
41. B. Błasiak and M. Cho, J. Chem. Phys., 2014, 140, 164107.
42. B. Błasiak and M. Cho, J. Chem. Phys., 2015, 143, 164111.
43. B. Błasiak, A. W. Ritchie, L. J. Webb and M. Cho, Phys. Chem. Chem. Phys., 2016, 18, 18094–18111.
44. B. Błasiak, C. H. Londergan, L. J. Webb and M. Cho, Acc. Chem. Res., 2017, 50, 968–976.
45. R. J. Xu, B. Blasiak, M. Cho, J. P. Layfield and C. H. Londergan, J. Phys. Chem. Lett., 2018, 9, 2560–2567.
46. B. Błasiak, M. Chołuj, J. D. Bednarska, R. W. Góra and W. Bartkowiak, globulion/eopdev: EOPDev 1.0.0, 2020-11-19, 2020, DOI: http://10.5281/zenodo.4281062.
47. J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein, F. A. Evangelista, J. T. Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke, M. L. Abrams, N. J. Russ, M. L. Leininger, C. L. Janssen, E. T. Seidl, W. D. Allen, H. F. Schaefer, R. A. King, E. F. Valeev, C. D. Sherrill and T. D. Crawford, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 556–565.
48. E. R. Davidson, J. Comput. Phys., 1975, 17, 87–94.
49. B. Liu, Numerical Algorithms in Chemistry: Algebraic Methods, Lawrence Berkeley Laboratory, University of California, Berkeley, CA, USA, 1978, pp. 49–53.
50. R. M. Parrish, E. G. Hohenstein and T. J. Martínez, J. Chem. Theory Comput., 2016, 12, 3003–3007.
51. P. Jordan, P. Fromme, H. T. Witt, O. Klukas, W. Saenger and N. Krauß, Nature, 2001, 411, 909.
52. A. S. Belov, D. V. Khokhlov and V. V. Poddubnyi, Dokl. Phys. Chem., 2016, 468, 63–66.
53. W. J. Hehre, R. Ditchfield and J. A. Pople, J. Chem. Phys., 1972, 56, 2257–2261.
54. F. Weigend, Phys. Chem. Chem. Phys., 2002, 4, 4285–4291.
55. E. F. Valeev, V. Coropceanu, D. A. da Silva Filho, S. Salman and J.-L. Brédas, J. Am. Chem. Soc., 2006, 128, 9882–9886.

---

Footnote

† Electronic supplementary information (ESI) available: Diagonal Hamiltonian matrix elements, timings of various variants of TI/CIS EET coupling calculations, discussion of basis set extension effects, and Cartesian coordinates of the studied model systems. See DOI: 10.1039/d0cp04636a

---