Theoretical description of electronically excited states of molecular aggregates at an ab initio level is computationally demanding. To reduce the computational cost, we propose a model Hamiltonian approach that approximates the electronically excited state wavefunction of the molecular aggregate. We benchmark our approach on a thiophene hexamer, as well as calculate the absorption spectra of several crystalline non-fullerene acceptors, including Y6 and ITIC, which are known for their high power conversion efficiency in organic solar cells. The method qualitatively predicts the experimentally measured spectral shape, which can be further linked to the molecular arrangement in the unit cell.

In molecular aggregates of large π-conjugated molecules, the excitons can diffuse over considerable distances because of large intermolecular excitonic couplings. Therefore, predicting the spectral properties of such aggregates is challenging because a large number of excited electronic states need to be considered, in addition to the influence of molecular packing1–4 and vibronic effects.5–7 High-level quantum chemistry methods such as coupled-cluster or configuration-interaction are often required for accurate predictions of excited state energies of single molecules. Because of their scaling as N5−7 with the number of electrons N, these methods are difficult to apply to large molecules or molecular aggregates.8 Better scaling is offered by Green’s function approach (N3−4)9 and the time-dependent density-functional theory (TDDFT) (N2);10,11 however, these methods are still too demanding for solid-state calculations of molecular films.

To overcome these difficulties, we adopt a tailored for excited states model Hamiltonian (MH) approach.12 The coarse-grained electronic Hamiltonian is constructed using several approximations: First, interactions of diabatic states on different molecules are considered only for molecular dimers. Second, the diabatic states, and the corresponding couplings, are used as a basis set to construct the excited state wavefunctions of the aggregate.13,14 Finally, the exciton-phonon coupling is approximated by Gaussian-broadening of the discrete energy spectrum.

In our approach, the interaction of diabatic states is evaluated at the ab initio level, and we can obtain all states of interest at a reasonable computational cost. Moreover, higher-order terms, i.e., couplings involving three or four monomers, can be included to systematically improve the accuracy of the MH.

We use the proposed MH approach to predict the absorption spectra of large molecular aggregates of non-fullerene acceptors (NFAs) and to link the spectral shape to the molecular arrangement in the unit cell. Recent developments of novel NFAs use the width of the NFA absorption spectra as a design criterion: broader absorption leads to larger short-circuit currents.15–19 One such example is Y6 with its narrow optical gap.20,21 This practical example shows that the use of MH is essential for material pre-screening, where the computational cost is often a bottleneck.22–24 

To construct the MH, we first expand the excited states ΨJ as a linear combination of a set of quasi-diabatic states ΞI. The choice of these diabatic states is the key to obtaining a computationally efficient and accurate MH. In partially ordered organic semiconducting films, the low-energy excited states are primarily localized on single molecules or immediate neighbors. Therefore, we choose diabatic states that reflect this spatial localization of excitons and charges: local excitations (LE) of isolated molecules and charge transfer (CT) states of molecular pairs.13 

To construct such diabatic states, we localize occupied (i, j, k, …) and virtual (a, b, c, …) molecular orbitals on the spatial domains (molecules) αi, αj, … and αa, αb, …. The configuration subspace, X = (αi, αa), in which the corresponding excitation configuration, ΦiaX, indicates that an electron is excited from occupied molecular orbitals (MOs) i of domain αi to virtual MO a of domain αa, where it can be an LE configuration if αi = αa or a CT configuration if αiαa.

The adiabatic excited state can be written as a linear combination of MO-based excitation configurations,

ΨJ=XiaXCJ,iaXΦiaX,
(1)

where the coefficients CJ,iaX are obtained by the character analysis method.13 

The adiabatic and diabatic states are related via a unitary transformation UJI,

ΞI=JΨJUJI.
(2)

We choose the transformation matrix U such that the diabatic states ΞI are as localized as possible onto the configuration subspace X. To do this, we introduce weights for the diabatic states,

P̃IX=iaJCJ,iaXUJI2=JKUJIPJKXUKI,
(3)

where

PJKX=iaXCJ,iaXCK,iaX
(4)

is the weight of a specific diabatic state. P̃IX1 indicates localization onto a single configuration space.

The transformation matrix UJI should then maximize the function

D=XI(P̃IX)2,
(5)

where D reflects the overall quality of the diabatization scheme. It is always smaller or equal to the number of diabatic states, and the closer its value is to the number of diabatic states, the higher the quality of the diabatization scheme.

With the transformation matrix U at hand, and the selected adiabatic states excitation energy EIJ = ϵIδIJ, where δIJ is the Kronecker delta function, we can back-transform the energy matrix, EIJ, into a coupling matrix,

EIJ=ϵI,ifI=J,EIJ=0,ifIJ,
(6)
CIJ=JKUIITEIJUJJ,
(7)

where the diagonal elements and off-diagonal elements of matrix C are the diabatic excitation energies and the coupling parameters, respectively. Correspondingly, we can get the transition dipole DI of diabatic states by transforming the transition dipoles of adiabatic states as follows:

DI=UIITDI.
(8)

Figure 1 illustrates the construction of the MH using the HOMO to LUMO transitions in a molecular dimer, with four different types of excitations.

FIG. 1.

Scheme of locally excited (LE, Φh1l1 and Φh2l2) and charge transfer (CT, Φh1l2 and Φh2l1) states of a molecular dimer.

FIG. 1.

Scheme of locally excited (LE, Φh1l1 and Φh2l2) and charge transfer (CT, Φh1l2 and Φh2l1) states of a molecular dimer.

Close modal

The coupling matrix, EMH=Φ|Ĥ|Φ, and the corresponding MH matrix read as follows:

formula
(9)

where the subscripts and superscripts, 1 and 2, denote the index of the domain that the electron is excited from and to, respectively. The diagonal elements E11 and E22 are the excitation energies of the two LE excited states Φh1l1 and Φh2l2, respectively; and E12 and E21 are the excitation energies of the two CT excited states, Φh1l2 and Φh2l1, respectively. The off-diagonal elements are the coupling parameters of the corresponding diabatic states. These coupling parameters can take four different forms: Vec, De, Dh, and W. Because these terms only involve two domains (monomers), we call them two-domain terms. Referring to the description of the four different diabatic states at the simplest excitation picture, i.e., HOMO to LUMO excitation, in Fig. 1, we can directly tell the physical meanings of these four coupling parameters: Vec is the coupling between the two LE states, Φh1l1 and Φh2l2, i.e., the excitonic coupling parameter; De(De) is the coupling between Φh2l2(Φh1l1) and Φh2l1(Φh1l2), where the two coupled states have the same occupation on the occupied MOs but different occupation on the virtual MOs, which describes the “electron” transport strength between the two molecules; similar to De(De), Dh(Dh) describes the “hole” transport strength between the two molecules; and W is the coupling between the two CT states. Explicit expressions of these coupling elements, in terms of monomer orbitals, are available in the supplementary material, to show how these can be approximated. In this work, these coupling elements are directly evaluated through the diabatization scheme. As an illustration, the complete coupling matrix of a tetramer system, including three- and four-domain terms, is also described in the supplementary material. All the high order terms have at least one density function that includes different monomers. For systems where the CT state excitation energies are larger than those of the LE states, which is often the case for mono-component organic crystals, three- and four-body terms can be neglected.

Given all this information, we can now construct the MH matrix M of an aggregate. The aggregate itself is a supercell containing several sets of equivalent dimers, so the diabatization calculations only need to be performed once for each set. In addition, Förster’s equation25 is used to screen far apart dimers with negligible coupling. This is done to optimize the coupling cutoff distance and save on computational costs. In the MH matrix M of the aggregate, each diagonal element is the excitation energy of an LE or a CT state, whereas the off-diagonal elements correspond to the coupling elements between the diabatic states. Because the matrix elements are obtained from the dimer diabatization, the diagonal and off-diagonal elements corresponding to CT excitations are unique. However, this is not always the case for LE states since the LE state localized on monomer I can result from the diabatization of any dimer that contains monomer I. In cases where these LE states are qualitatively close to single molecule excitations, the corresponding excitation energies, and excitation character, of the LE states resulting from different dimers will be similar. In the MH method, we use the average excitation energy of the LE states from all equivalent dimers.

The eigenvalues of M are then the excitation energies of the aggregate,

EMH=UTMU.
(10)

The transition dipole moment of each MH state is then obtained by transforming the transition dipole moment matrix of diabatic states with the matrix U′,

DMH=UTDdiabatic,
(11)

which can be used to calculate the corresponding oscillator strength.

To ensure that the phase of the localized exciton is preserved for the same molecule belonging to different dimers, we set up a unitary matrix that connects the geometry and the transition dipole tensors of the molecule in the dimer. The unitary matrix is then used for transforming the dimers’ coupling matrix to ensure the signs of the coupling parameter are consistent in the MH matrix.

The quality of the results obtained from the MH method is highly dependent on the diabatization scheme and the chosen diabatization basis. In our scheme, the MH is constructed by using a basis of LE and CT states that are localized on not more than two molecules. The MH should, in principle, be accurate for aggregates of weakly coupled organic molecules, where the lowest-lying CT state is always above the lowest-lying LE state. This choice of the diabatic basis results in high quality MH wavefunctions but also has its limitations. For example, for some closely spaced dimers, it can be difficult to get reasonable diabatization results, where a larger number of converged adiabatic excited states is required. On the other hand, these closely spaced dimers are of the most importance for constructing the MH wavefunctions. For the same reason, the scheme can also be inefficient for strongly coupled monomers and open shell systems.

The crystal structures of five organic NFAs (Y6, ITIC, ITIC-4F, IEICO-4F, and IDTBR) were first optimized using the Vienna Ab initio Simulation Package (VASP) package.26–30 The PBE functional31 and the projector-augmented wave (PAW)27,30 method were applied to treat the electron–electron and the valence–core interactions, respectively. In addition, the empirical D3 dispersion correction term32 was added. A kinetic-energy cutoff of 520 eV was used for the plane wave basis set. The Brillouin zone was sampled using a 4 × 4 × 4 Monkhorst–Pack k-point scheme. The total energy convergence criterion was 10−4 eV. After the optimization, the Hellmann–Feynman forces were smaller than 0.01 eV/Å.

The optimized unit cells were then used to construct molecular aggregates of 11 × 11 × 11 unit cells. For every dimer in the aggregate, we estimated the transition dipole–dipole coupling Vec by using the Förster equation,25 

V=DaDbRab33(DaRab)(RabDb)Rab5,
(12)

where Da is the transition dipole vector of monomer a and Rab is the vector between the two monomers’ transition dipole centers.

As Eq. (12) is known to underestimate the coupling strength between the neighboring pairs of molecules, it is only used to select the pairs for which coupling calculations are performed. The actual coupling strengths are then calculated at the QM level for all dimers within a 3 × 3 × 3 supercell. Since dimers in distinct sets are identical to each other in terms of electronic and excitonic couplings, the diabatization procedure was performed only for one dimer in the set. For each set, we make sure that at least two pairs of the lowest energy diabatic LE states are sufficiently obtained, and for some sets, the lowest lying CT states pair is also included. As mentioned in Eq. (3), the character weights of the diabatic states can be used as an indicator of the diabatization quality. In our calculations, we normally calculated the ten lowest energy excited states, and the character weights of the above-mentioned diabatic states are always larger than 98%.

The MH matrix for a large supercell, up to 11 × 11 × 11, was then constructed by using the couplings of either the diabatization calculations if the dimer is calculated at the QM level, or with an estimated Vec[Eq. (12)] otherwise. The effect of the choice of cutoff is discussed further in Sec. IV.

All the dimers ab initio calculations were performed at the cam-b3lyp/TZVP33,34 level of theory using the TURBOMOLE 7.535 package. The excited state energies of the NFAs were corrected using the GW+BSE method.36 

To benchmark the method, we calculated the absorption spectra of a thiophene hexamer, for which both the MH method and the direct calculation are possible at the same level of theory. We placed six monomers into two layers, each composed of three monomers, as shown in Fig. 2(a). In this molecular arrangement, with an interlayer distance of 3.5 Å, some CT states have relatively low energies.

FIG. 2.

Molecular arrangements used to calculate the absorption spectra.

FIG. 2.

Molecular arrangements used to calculate the absorption spectra.

Close modal

It is known that TDDFT, without using range-separated functionals, fails to predict the energy levels of CT states.37 As CT states are essential for the MH method, the use of range-separated functionals is recommended. This can be illustrated in the thiophene hexamer: The hybrid functional b3lyp underestimates the LE character of the low-lying excited states, while another hybrid functional, bhlyp, and the range-separated hybrid functionals, cam-b3lyp33 and wb97m-v,38 all produce very similar results, where the wb86m-v predicted slightly stronger LE character (Table S1).

We performed diabatization analyses for 15 dimers and the entire hexamer. Each thiophene molecule has a pair of nearly degenerate low-lying LE states, with energies of 5.97 and 6.04 eV. π-stacked dimers have two additional pairs of low-lying CT states, with energies of 6.5 and 6.8 eV. The rest of the diabatic states either have very large energies, more than 7 eV, or their couplings to the low-lying states are small (<1 meV). From these diabatization analyses, 24 diabatic states were chosen to construct the MH.

By diagonalizing the MH matrix, we obtain the excitation energies of the first 24 adiabatic states. The energies and corresponding oscillator strengths of the first 12 states are compared to the results obtained for the entire hexamer using TDDFT in Table I. As we can see, the excitation energies differ only slightly, with differences varying from 11 to 81 meV. The oscillator strengths also follow a similar trend. Therefore, we can conclude that the MH method can accurately reproduce low lying excited states of an aggregate.

TABLE I.

Excitation energies (eV) of 12 lowest-lying adiabatic states of the thiophene hexamer and the corresponding oscillator strengths calculated by TDDFT, as well as using MHs with and without CT excitations. All calculations are performed at the cam-b3lyp/TZVP level of theory.

Without CT
TDDFTWith CT(LE only)
ΔEfΔEfΔEf
5.387 0.000 03 5.407 0.000 13 5.855 0.000 51 
5.396 0.000 16 5.413 0.000 09 5.879 0.000 38 
5.410 0.000 18 5.456 0.000 78 5.923 0.001 01 
5.419 0.000 16 5.461 0.000 01 5.967 0.000 06 
5.490 0.000 16 5.511 0.000 01 5.970 0.021 47 
5.526 0.000 28 5.560 0.000 09 5.981 0.011 24 
5.993 0.119 35 6.031 0.174 64 6.058 0.146 41 
6.013 0.342 70 6.095 0.380 72 6.107 0.436 86 
6.040 0.077 71 6.021 0.030 28 6.033 0.028 60 
10 6.058 0.153 88 6.038 0.078 13 6.062 0.116 78 
11 6.075 0.009 99 6.086 0.011 11 6.097 0.049 15 
12 6.087 0.151 66 6.098 0.055 47 6.109 0.003 43 
Without CT
TDDFTWith CT(LE only)
ΔEfΔEfΔEf
5.387 0.000 03 5.407 0.000 13 5.855 0.000 51 
5.396 0.000 16 5.413 0.000 09 5.879 0.000 38 
5.410 0.000 18 5.456 0.000 78 5.923 0.001 01 
5.419 0.000 16 5.461 0.000 01 5.967 0.000 06 
5.490 0.000 16 5.511 0.000 01 5.970 0.021 47 
5.526 0.000 28 5.560 0.000 09 5.981 0.011 24 
5.993 0.119 35 6.031 0.174 64 6.058 0.146 41 
6.013 0.342 70 6.095 0.380 72 6.107 0.436 86 
6.040 0.077 71 6.021 0.030 28 6.033 0.028 60 
10 6.058 0.153 88 6.038 0.078 13 6.062 0.116 78 
11 6.075 0.009 99 6.086 0.011 11 6.097 0.049 15 
12 6.087 0.151 66 6.098 0.055 47 6.109 0.003 43 

We can also approximate the wave function of the hexamer by removing the CT states in the MH matrix. This results in the excitation energies of the LE states becoming significantly larger (>0.4 eV), as shown in Table I, suggesting that these LE states are strongly coupled to the CT states, resulting in lower excitation energies. The six higher lying adiabatic states in the LE-only MH (7th to 12th) have only slightly larger excitation energies than observed in the LE+CT MH (about 0.02 eV). This shows that the higher lying LE states are coupled to the CT states as well; however, this coupling is weaker than that observed in the six lowest lying adiabatic states in the LE-only MH. The corresponding oscillator strengths also change, which implies that the absorption spectra would change their shape.

From the point of view of computational efficiency, TDDFT calculations of a hexamer for 24 excited states took 4080 s on 40 cores, whereas the calculation of 15 dimers, including the diabatization, took only 270 s on 40 cores. One should also note that the conventional TDDFT method scales as the square of the total number of atoms in the system, whereas the MH method only scales linearly with the number of symmetry unique dimers.

The unit cell of Y6, which is shown in Fig. 2(b), contains four molecules and two distinct conformers. The lowest lying excitations of these conformers have energies of 2.056 and 2.060 eV, with oscillator strengths of 2.34 and 2.44 a.u. The second excited states are at 2.58 and 2.59 eV, with oscillator strengths of 0.50 and 0.41 a.u. The GW+BSE correction predicts a 0.245 eV reduction of the lowest excited state energy of the Y6 molecule. This energy difference will be used for shifting the simulated absorption spectra.

For all dimers, the lowest energy diabatic states have an LE character with energies of 2.05 and 2.54 eV, very close to the two lowest energy excited states of the Y6 monomer. Only six dimer types have CT states with energies below 3.00 eV. These six dimer types show a strong J-aggregation character for which the lowest excitations are red-shifted by 0.03–0.08 eV with respect to the monomer's lowest excitation energy. Two dimer types display strong H-aggregation character, with 0.01–0.02 eV blue shifts. The rest of the dimers have shifts below 0.01 eV.

The absorption spectrum of Y6 aggregates calculated using the MH method is shown in Fig. 3, together with the absorption spectrum of a Y6 monomer and the experimentally determined absorption spectrum of a Y6 thin film.39 The calculated absorption spectrum is obtained by the Gaussian convolution of oscillator strengths with a full width at half maximum (FWHM) of 0.18 eV. The MH excitation energies with the GW+BSE correction are used to determine the peak positions. The peak with the lowest energy is used to normalize the spectra.

FIG. 3.

Absorption spectra of Y6 aggregates (MH) and thin films (experiment).

FIG. 3.

Absorption spectra of Y6 aggregates (MH) and thin films (experiment).

Close modal

The calculated peak intensities differ somewhat from those measured in the thin film, in particular, around 2.17 eV, where the intensities obtained from the MH method are a bit too strong. This is predominantly attributed to the fact that only a simple Gaussian convolution was used while simulating the spectra instead of considering the corresponding vibronic degrees of freedom. The simulated spectrum also has several sub-peaks at 1.67, 1.88, and 2.17 eV. Red arrows show the corresponding peaks in the measured spectra in Fig. 3(a).

To prove the convergence with the aggregate size, we show the spectra for several aggregates in Fig. 3(b). For larger aggregates, the intensities of the low-energy peaks increase, whereas those of the high-energy peaks decrease. However, the peak positions of the major absorption peaks do not change significantly. We can conclude that the spectra are practically unaffected by the aggregate size starting from an aggregate size of ∼10 nm (7 × 7 × 7 unit cells).

For comparison, we modified the MH by removing all CT states and recalculated the absorption spectra. As shown in Fig. 3(c), this spectrum exhibits two absorption peaks, with the first peak blue-shifting by about 0.12 eV compared to the absorption spectrum of the LE+CT MH. The second peak is red-shifting, and the peaks in the 1.4–2.2 eV range are significantly less pronounced. Therefore, we can conclude that the CT states, indeed, play an important role in the Y6 aggregates.

Finally, Fig. 3(d) shows the effect of the cutoff for the coupling elements. As one can see, the value of the cutoff does not affect the absorption spectrum much; this is as most of the inner layer (3 × 3 × 3 unit cells) couplings are always included in the MH.

At this point, we would like to mention the scaling of computational costs as compared to the traditional TDDFT calculations. Table II summarizes the dimensions of the MH matrix for different aggregate sizes. Conventional ab initio excited state calculations for such aggregates scale from the square to the fourth order of the number of atoms,40 whereas the MH method scales with the square of the matrix dimension. The construction of the MH of the Y6 aggregate requires a number of TDDFT dimer calculations and diabatization calculations. Depending on the intermolecular distances, a dimer calculation takes from 5 to 30 h on a 40-core cluster, where the larger intermolecular distance corresponds to the shorter time. In our current setup, a total number of 334 Y6 dimers have been calculated. The final MH step is done on a single core, and the CPU hours are listed in Table II.

TABLE II.

Number of molecules Nm, atoms Na, and the size of the MH matrix for aggregates corresponding to supercells of size na × nb × nc.

na, nb, ncNmNaMH sizeCPU time (s)
3,3,3 108 20 196 513 4.2 
5,5,5 500 93 500 2 727 295 
7,7,7 1372 256 564 7 949 8 529 
9,9,9 2 916 545 292 17 475 53 337 
na, nb, ncNmNaMH sizeCPU time (s)
3,3,3 108 20 196 513 4.2 
5,5,5 500 93 500 2 727 295 
7,7,7 1372 256 564 7 949 8 529 
9,9,9 2 916 545 292 17 475 53 337 

We also calculated the absorption spectrum of another NFA molecule, ITIC. The first and the second excited states of ITIC molecules extracted from the optimized crystal structure have excitation energies of 2.22 and 2.69 eV, with oscillator strengths of 3.05 and 0.00 a.u., respectively. The GW+BSE correction predicted a 0.215 eV energy shift for the lowest excited state. The excitation energies of other excited states are higher than 3 eV and, therefore, are not taken into account in the diabatic calculations. For ITIC dimers, we observed weaker J- and H-aggregation as compared to Y6 dimers. Four dimers have a J-aggregation character with about 0.03 eV red shift, and one dimer has a strong H-aggregation character, with a 0.04 eV blue shift. Other dimers have weak J- or H-aggregation character, with shifts smaller than 0.01 eV.

The absorption spectrum of an 11 × 11 × 11 supercell, with a 0.2 meV cutoff, is shown in Fig. 4, together with the experimentally measured thin film absorption spectrum.17 Although the shoulder peak at around 2.17 eV of the calculated spectra is smaller than the experimental one, the simulated and experimental spectra have a fairly similar shape. The cutoff and size dependence are shown in Figs. S9 and S10 of the supplementary material.

FIG. 4.

The simulated absorption spectrum of ITIC aggregates. Experimentally measured thin-film spectrum is blue-shifted by 0.1 eV to overlap the lowest absorption peaks.

FIG. 4.

The simulated absorption spectrum of ITIC aggregates. Experimentally measured thin-film spectrum is blue-shifted by 0.1 eV to overlap the lowest absorption peaks.

Close modal

The broader spectra of Y6 as compared to ITIC can be related to the J- and H-aggregates. The ITIC crystal unit cell contains only two molecules, with the normal vectors of the backbone planes forming an angle of about 60°. This packing pattern is simpler than the one of Y6. Moreover, Y6 ππ stacking distances are shorter than those of ITIC.

We propose a model Hamiltonian for the electronically excited state wavefunction. The wave function is approximated by a linear combination of a set of diabatic excited states. These states, and the corresponding couplings between them, are evaluated by diabatizing the excited state wavefunction of molecular dimers, which greatly reduces the computational cost. With the help of the MH, we calculate the electronically excited state wavefunctions of large molecular aggregates for a model system, a thiophene hexamer, and show that the error of the predicted excitation energies is below 0.08 eV. We also calculated the absorption spectra of aggregates of several NFAs, Y6, ITIC, ITIC-4F, IDTBR, and IEICO-4F, and compared them to the experimentally measured spectra. The general shape of the spectra and energy levels are well reproduced for most of the compounds, although there are differences in describing absorption peak intensities.

See the supplementary material for the model Hamiltonian with the higher order coupling elements, benchmarks with different density functionals, diabatization procedure details, and the system size dependencies.

This publication is based on work supported by the KAUST Office of Sponsored Research (OSR) under Award Nos. OSR-2018-CARF/CCF-3079 and OSR-CRG2018-3746. W.L. acknowledges the National Natural Science Foundation of China (Project No. 21703023) for the financial support. D.A. also acknowledges the KAUST PSE Division for hosting his sabbatical in the framework of the Divisions Visiting Faculty program. D.A. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for financial support through collaborative research centers (Grant Nos. TRR 146, SPP 2196, and 460766640). We thank Kun-Han Lin, Mukunda Mandal, Andriy Zhugayevych, and Naomi Kinaret for their fruitful discussions and proofreading of the manuscript.

The authors have no conflicts to disclose.

Wenlan Liu: Conceptualization (equal); Data curation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Denis Andrienko: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Resources (equal); Writing – original draft (equal); Writing – review & editing (equal).

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

1.
E. G.
McRae
and
M.
Kasha
,
J. Chem. Phys.
28
,
721
(
1958
).
2.
N. J.
Hestand
and
F. C.
Spano
,
J. Chem. Phys.
143
,
244707
(
2015
).
3.
L.
Craciunescu
,
S.
Wirsing
,
S.
Hammer
,
K.
Broch
,
A.
Dreuw
,
F.
Fantuzzi
,
V.
Sivanesan
,
P.
Tegeder
, and
B.
Engels
,
J. Phys. Chem. Lett.
13
,
3726
(
2022
).
4.
W.
Barford
and
M.
Marcus
,
J. Chem. Phys.
146
,
130902
(
2017
).
5.
W.
Popp
,
M.
Polkehn
,
K. H.
Hughes
,
R.
Martinazzo
, and
I.
Burghardt
,
J. Chem. Phys.
150
,
244114
(
2019
).
6.
A. M.
Alvertis
,
R.
Pandya
,
C.
Quarti
,
L.
Legrand
,
T.
Barisien
,
B.
Monserrat
,
A. J.
Musser
,
A.
Rao
,
A. W.
Chin
, and
D.
Beljonne
,
J. Chem. Phys.
153
,
084103
(
2020
).
7.
V.
Stehr
,
R. F.
Fink
,
B.
Engels
,
J.
Pflaum
, and
C.
Deibel
,
J. Chem. Theory Comput.
10
,
1242
(
2014
).
8.
W.
Zhang
,
J.
Liu
,
X.
Jin
,
X.
Gu
,
X. C.
Zeng
,
X.
He
, and
H.
Li
,
Angew. Chem.
132
,
11647
(
2020
).
9.
M.
Rohlfing
and
S. G.
Louie
,
Phys. Rev. B
62
,
4927
(
2000
).
10.
G.
Onida
,
L.
Reining
, and
A.
Rubio
,
Rev. Mod. Phys.
74
,
601
(
2002
).
11.
H.
Li
,
Y.
Zhao
,
J.
Fang
,
X.
Zhu
,
B.
Xia
,
K.
Lu
,
Z.
Wang
,
J.
Zhang
,
X.
Guo
, and
Z.
Wei
,
Adv. Energy Mater.
8
,
1702377
(
2018
).
12.
H.
Fidder
,
J.
Knoester
, and
D. A.
Wiersma
,
J. Chem. Phys.
95
,
7880
(
1991
).
13.
W.
Liu
,
B.
Lunkenheimer
,
V.
Settels
,
B.
Engels
,
R. F.
Fink
, and
A.
Köhn
,
J. Chem. Phys.
143
,
084106
(
2015
).
14.
W.
Liu
,
S.
Canola
,
A.
Köhn
,
B.
Engels
,
F.
Negri
, and
R. F.
Fink
,
J. Comput. Chem.
39
,
1979
(
2018
).
15.
Z.
Zhang
,
J.
Yuan
,
Q.
Wei
, and
Y.
Zou
,
Front. Chem.
6
,
00414
(
2018
).
16.
C.
Yan
,
S.
Barlow
,
Z.
Wang
,
H.
Yan
,
A. K.-Y.
Jen
,
S. R.
Marder
, and
X.
Zhan
,
Nat. Rev. Mater.
3
,
18003
(
2018
).
17.
Y.
Wang
,
Q.
Fan
,
X.
Guo
,
W.
Li
,
B.
Guo
,
W.
Su
,
X.
Ou
, and
M.
Zhang
,
J. Mater. Chem. A
5
,
22180
(
2017
).
18.
P.
Meredith
,
W.
Li
, and
A.
Armin
,
Adv. Energy Mater.
10
,
2001788
(
2020
).
19.
K.
Vandewal
,
S.
Mertens
,
J.
Benduhn
, and
Q.
Liu
,
J. Phys. Chem. Lett.
11
,
129
(
2020
).
20.
Y.
Yang
,
ACS Nano
15
,
18679
(
2021
).
21.
L.
Perdigón-Toro
,
H.
Zhang
,
A.
Markina
,
J.
Yuan
,
S. M.
Hosseini
,
C. M.
Wolff
,
G.
Zuo
,
M.
Stolterfoht
,
Y.
Zou
,
F.
Gao
,
D.
Andrienko
,
S.
Shoaee
, and
D.
Neher
,
Adv. Mater.
32
,
1906763
(
2020
).
22.
K.-H.
Lin
,
L.
Paterson
,
F.
May
, and
D.
Andrienko
,
npj Comput. Mater.
7
,
179
(
2021
).
23.
K.-H.
Lin
,
G.-J. A. H.
Wetzelaer
,
P. W. M.
Blom
, and
D.
Andrienko
,
Front. Chem.
9
,
800027
(
2021
).
24.
N. C.
Forero-Martinez
,
K.-H.
Lin
,
K.
Kremer
, and
D.
Andrienko
,
Adv. Sci.
9
,
2200825
(
2022
).
25.
T.
Förster
,
Ann. Phys.
437
,
55
(
1948
).
26.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
558
(
1993
).
27.
G.
Kresse
and
J.
Furthmüller
,
Comput. Mater. Sci.
6
,
15
(
1996
).
28.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
29.
G.
Kresse
and
J.
Hafner
,
J. Phys.: Condens. Matter
6
,
8245
(
1994
).
30.
G.
Kresse
and
D.
Joubert
,
Phys. Rev. B
59
,
1758
(
1999
).
31.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
32.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
,
J. Chem. Phys.
132
,
154104
(
2010
).
33.
T.
Yanai
,
D. P.
Tew
, and
N. C.
Handy
,
Chem. Phys. Lett.
393
,
51
(
2004
).
34.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
35.
TURBOMOLE V7.5 2021, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com.
36.
K.
Krause
and
W.
Klopper
,
J. Comput. Chem.
38
,
383
(
2017
).
37.
A.
Dreuw
and
M.
Head-Gordon
,
J. Am. Chem. Soc.
126
,
4007
(
2004
).
38.
N.
Mardirossian
and
M.
Head-Gordon
,
J. Chem. Phys.
144
,
214110
(
2016
).
39.
H.
Yu
,
Z.
Qi
,
J.
Zhang
,
Z.
Wang
,
R.
Sun
,
Y.
Chang
,
H.
Sun
,
W.
Zhou
,
J.
Min
,
H.
Ade
, and
H.
Yan
,
J. Mater. Chem. A
8
,
23756
(
2020
).
40.
A.
Dreuw
and
M.
Head-Gordon
,
Chem. Rev.
105
,
4009
(
2005
).

Supplementary Material