Geometry dependence of excitonic couplings and the consequences for configuration‐space sampling

Excitonic coupling plays a key role for the understanding of excitonic energy transport (EET) in, for example, organic photovoltaics. However, the calculation of realistic systems is often beyond the applicability range of accurate wavefunction methods so that lower‐scaling semi‐empirical methods are used to model EET events. In the present work, the distance and angle dependence of excitonic couplings of dimers of selected organic molecules are evaluated for the semi‐empirical long‐range corrected density functional based tight binding (LC‐DFTB) method and spin opposite scaled second order approximate coupled cluster singles and doubles (SOS‐CC2). While semi‐empirically scaled methods can lead to slightly increased deviations for excitation energies, the excitonic couplings and their dependence on the dimer geometry are reproduced. LC‐DFTB yields a similar accuracy range as density‐functional theory (DFT) employing the ωB97X functional while the computation time is reduced by several orders of magnitude. The dependence of the exchange contributions to the excitonic couplings on the dimer geometry is analyzed assessing the calculation of Coulombic excitonic couplings from monomer local excited states only, which reduces the computational effort significantly. The present work is a necessary first step toward the simulation of excitonic energy transport using semi‐empirical methods.

sent work, the distance and angle dependence of excitonic couplings of dimers of selected organic molecules are evaluated for the semi-empirical long-range corrected density functional based tight binding (LC-DFTB) method and spin opposite scaled second order approximate coupled cluster singles and doubles (SOS-CC2). While semi-empirically scaled methods can lead to slightly increased deviations for excitation energies, the excitonic couplings and their dependence on the dimer geometry are reproduced. LC-DFTB yields a similar accuracy range as density-functional theory (DFT) employing the ωB97X functional while the computation time is reduced by several orders of magnitude. The dependence of the exchange contributions to the excitonic couplings on the dimer geometry is analyzed assessing the calculation of Coulombic excitonic couplings from monomer local excited states only, which reduces the computational effort significantly. The present work is a necessary first step toward the simulation of excitonic energy transport using semi-empirical methods.

| INTRODUCTION
In recent years, organic semiconductors have been widely investigated and applied in electronic and photonic applications. 1,2 Organic materials are of particular interest for such applications since they combine different desirable properties. They have electronic properties such as low band-gaps that are important for these applications, 1,3 they are cheap to produce and can be molded in various forms or used as thin films due to their flexibility. 3,4 For example, systems with conjugated π-systems and aromatic rings have been investigated for their use in organic photovoltaics. 2,4,5 To be able to use sunlight as a source for electric energy, light needs to be converted into electric current. Electron donor molecules that are electronically excited transport the absorbed energy to an interface where electric current is generated due to charge separation.
This energy transport proceeds without charge-transfer (CT) and is known as excitonic energy transport (EET). In the Frenkel exciton model, 6 excitons are described via the interactions between excited states, which is also denoted excitonic coupling. 7 Excitonic couplings have been studied using a wide range of methods [8][9][10][11][12][13][14] ranging from highly accurate to fast, approximate methods. Different classes of molecules and excitonic couplings in different phases have been studied. 2,13,[15][16][17] Recently, density functional theory (DFT) and density functional based tight binding (DFTB) methods have been benchmarked for acenes and small molecules. 18,19 In order to assess the performance of different methods, a reference method has to be chosen to meet certain needs. First, results have to be consistent and accurate. Second, the method must not be computationally too expensive, as dimers of systems including multiple aromatic rings need to be feasible. For instance, the second order approximate coupled cluster singles and doubles (CC2) method scales with O N 5 À Á , which at present means that dimers of pentacene can still be calculated. It is chosen as the reference method since it is an accurate method leading to consistent results for excitation energies compared to CC3. 20 The semi-empirically scaled SOS-CC2 is of interest for larger molecular systems since in its Laplace-transformed formalism (LT-SOS-CC2), it scales with O N 4 À Á , one order of magnitude smaller than CC2. [21][22][23] The algebraic diagrammatic construction (ADC) scheme for the polarization propagator 24,25 provides a series of ab-initio methods with which excitation energies can be calculated using perturbation theory, following the typical partitioning of the Hamiltonian by Møller and Plesset. 26 At second order, ADC(2) provides excitation energies for singles excitations consistent up to second order and their accuracy is very similar to those of CC2. 27 However, ADC(2) bares the advantage of being not only size-consistent but also Hermitian. 28  shown to reliably predict some excited-state properties, 32-35 they show large errors concerning CT processes. 33,36,37 Because of this, systems with extended conjugation [38][39][40] or systems in which CT processes are present are often not described accurately. 41 Improved functionals, such as range-separated functionals that include a longrange correction of the exchange potential, help mitigate this problem. 42 Long-range corrected (LC) functionals have been tested on their performance with respect to EET processes, where the family of ωB97 functionals has shown to outperform other functionals. 43 However, these functionals do not yield a systematic improvement as well, so individual benchmarking is always necessary.
DFTB 44,45 is an alternative for very large molecular systems, since it is about two to three orders of magnitude faster than DFT using the generalized gradient approximation (GGA) and mediumsized basis sets. It can be derived from DFT using a Taylor series expansion around a reference density, followed by approximations for the two-electron integrals. 46 Besides these approximations, the accuracy of DFTB is also limited by the applied exchange-correlation functional, which was based on GGA. Only recently, a LC functional has been implemented in the framework of DFTB for the ground state 47 and extended to compute excited states properties via the linear-response formalism, 48 abbreviated LC-DFTB in the following.
It has been benchmarked for excitation energies of a large set of organic molecules including charge-transfer excitations. 48 For the application to exciton-transfer, Coulomb couplings have been investigated. 49 where E a and E b denote excitation energies for local excited (diabatic) states on monomers A and B, respectively. Analytical diagonalization of this 2 Â 2 matrix yields the excitation energies of the (adiabatic) dimer states E m and E n , and their difference is given as: Re-arranging Equation (2) shows the influence of the energy difference ΔE nm upon the excitonic coupling J: In case of symmetric dimers, that is, ΔE ab = 0, the coupling J simplifies to: The difference ΔE nm obtained from supermolecular calculations, however, does not correspond precisely to the target coupling J(m, n) of the two states in general. Equation (4) where X and Y are the atoms of the corresponding molecules. Q 0a and Q 0b are the atomic Mulliken transition charges for the excited states a and b, respectively. The pure atomic Coulomb interaction ζ XY is defined as where in which μ(r) denote basis functions and N X is the total number of basis functions on atom X. ζ XY can be calculated analogously to the electronelectron interaction term γ XY in the LC-DFTB formalism. 47 For further details of this approach the reader is referred to Ref 49. If the coupling obtained from supermolecular excited states m and n (Equation (4)) is in sufficient agreement with the Coulomb coupling of two local states a and b according to Equation (5), that is J (m, n) ≈ J C (a, b), then exchange effects are negligible. Note, however, that the Mullikan transition charges in Equation (5) are an additional approximation so that even for increased distances small deviations may remain, vide infra. As the calculation of J C (a, b) is significantly faster compared to J(m, n), using only the former would provide a computational significant advantage in the simulation of large molecular systems.

| Geometries
The dimers required for the study of interacting local states were constructed by optimizing the geometries of the monomers with B3LYP and stacking two monomers symmetrically face to face. In bulk material, however, relative positions between monomers differ by, for example, intermolecular distance and rotation angle ( In the following, two sets of molecules are discussed separately: First set A, consisting of different acenes and then, set B, consisting of guanine and purine as well as different naphthyridines ( Figure 2) in order to assess the influence of nitrogen substitution effects. In LC-DFTB the local BNL functional 68,69 was used for the shortrange part and a conventional non-local Hartree-Fock exchange for ω ! ∞ for the long-range part. 48 The range-separation parameter ω is set to ω = 0.3/a 0 for the computation of the electronic parameters. 48 DFTB uses a minimal atomic orbital basis set, which is computed from atomic Kohn-Sham equations, and an additional harmonic potential is introduced to confine the basis function. The harmonic potential is characterized by confinement radii r 0 , which determine the range of the potential and therefore the extension of the LCAO basis functions. The radii r 0 are usually optimized for properties such as atomization energies, geometries and vibrational frequencies of molecules, resulting typically in values being a factor of 2 of the covalent radii of the corresponding atoms. 45 The original LC-DFTB confinement radii lead to accurate vertical excitation energies and Coulomb couplings employing rather compact atomic orbital basis sets. 48

| Method assessment
Due to the size of the dimer systems and the large amount of distances and angles, we have chosen to use CC2 as general benchmark method. However, as CC2 is only a second-order method, its accuracy is briefly assessed with respect to ADC(3). For the present purpose excitation energies and supermolecular couplings at selected distances and angles of pyrrol, pyridine, and pyrazine dimers are compared, see Table 1, where the excitonic coupling is given as For the investigated properties ADC(3) and CC2 yield practically identical results. For example, the couplings J obtained with LC-DFTB yield a mean error of +9 meV and standard deviation of 25 meV with respect to ADC(3) and a mean error of +7 meV and a standard deviation of 21 with respect to CC2, respectively. Based on these results and the small deviations between ADC(3) and CC2, we will evaluate SOS-CC2 and SCS-CC2 with respect to CC2.

| Sampling
For sampling, a large number of dimer geometries were taken from molecular dynamics (MD) simulations of an anthracene crystal containing 10 Â 40 Â 5 molecules along the respective crystallographic with a reorganization energy λ of 0.563 eV 70 is used, where σ is the standard deviation. The latter accounts for fluctuations around the average structure.
A value to quantify the degree of influence of the dynamics is the coherence parameter, 75,76 reaching values near 1 or 0 when the coupling is defined by the average structure or by non-equilibrium conformations, respectively. The mean-square displacement (MSD) of the exciton averaged over 10,000 trajectories (N traj ) is calculated as: 77 where x l A t ð Þ is the center of mass of molecule A along trajectory l, P l A t ð Þ is the corresponding diabatic population and x l (0) is the center of the exciton at the start of the simulation (t = 0). The diffusion constant D was finally calculated as follows: In case of the model dimer systems the coupling can be obtained as half the energy difference of the coupled states, c.f. Figure 3. The underlying condition is, however, that the two monomers have identical geometries and thus degenerate excited states, so that Equation (4) can be used. In case of the MD simulation the two monomers are no longer identical and do not possess degenerate excitation energies so that Equation (3) has to be used. To obtain the individual monomer energies, for each snapshot not only the dimers have to be calculated but also the individual monomers using the parameter set 2 optimized for coupling. It must be pointed out, however, that despite using Equation (3) still numerical problems can occur. If the excitation energy gap of the dimer is smaller than the difference of monomer excitation energies, the square-root term becomes negative and the coupling becomes imaginary. This problem may be rooted in the approximate nature of the overall approach, which assumes that the dimer states are a linear combination of the two corresponding monomer states only. However, only a small amount of about 0.5% of the couplings turned out to be imaginary and the corresponding snapshots were neglected for the analysis.

| RESULTS
The results are ordered as follows. We begin by assessing the perfor- Individual results for the molecules are given in the Supporting Information. In Tables 2 and 3  N that can be reduced by one order of magnitude, that is, O N 4 À Á instead of O N 5 À Á when the Laplace transformation is used. [21][22][23] TDDFT employing the range-separated ωB97X functional yields excitation energies which are 0.3-0.4 eV higher than those obtained T A B L E 2 Errors in eV of the excitation energies E 1 and E 2 belonging to states S 1 and S 2 , respectively, compared to results obtained using CC2 of set A summed over all investigated distances r and angle φ = 0 T A B L E 3 Errors in eV of the excitation energies E 1 and E 2 belonging to states S 1 and S 2 , respectively, compared to results obtained using CC2 of set A at r = 4 Å and summed over all investigated angles  The LC-DFTB values in Tables 2 and 3 were obtained using parameter set 1, that is, with confined AOs, vide supra. The mean deviations seem to be small but the standard deviation shows that excitation energies are both overestimated and underestimated to a similar extent, the overall accuracy is comparable to that of LC-DFT.

| Excitonic couplings: Distance dependence
Having addressed the accuracy of the excitation energies, we now turn to excitonic couplings and, in particular, their dependence on intermolecular distances and angles. For example, in Figure  The couplings as obtained in this manner are discussed in the following.
Results for set A are collected in  Figure 4. The table reveals that ADC(2) has a mean error of about +0.08 eV in case of absolute excitation energies, while the mean error of the couplings J is less than 1 meV. The performance of a method with respect to excitation energies does therefore not necessarily directly translate to its performance for predicting excitonic couplings as long as the electronic structure of the state is sufficiently accurately described.
In case of the scaled CC2 variants, SCS-CC2 and SOS-CC2, the mean errors of the coupling are slightly increased compared to ADC (2)/def2-SVPD. In particular the standard deviation and the maximum errors are increased by one order of magnitude for the absolute (and relative) errors. In case of DFT and LC-DFTB the errors are again somewhat increased. LC-DFTB shows a better accuracy than DFT employing the ωB97X functional. Note, however, that in case of DFTB the couplings must be computed with a different parameter set than the absolute excitation energies. The observed deviations are overall small compared to the absolute coupling strength, so that different decay rates are within the margin of error. For example, plotting the distance dependence of total couplings will exhibit no visible effect, cf. the analysis for naphthyridines in Figure 9. As LC-DFTB and ωB97X show deviations from CC2 in a similar range which we expect for other DFT-functionals as well, we decided not to further interpret these differences.

| Excitonic couplings: Rotation-angle dependence
Results for the dependence of the couplings on the rotation angle are also collected in Table 4, and displayed for the reference method CC2 in Figure 5. While in most cases the coupling decreases from 0 to 90 steadily, naphthalene and pyrene exhibit a minimum in the coupling at 30 and a maximum at 60 , cf. Section 3.1.4. The coupling strengths at 60 decrease with increasing length of the linear acenes, that is, from naphthalene to pentacene. With increasing chain length, however, the couplings at 30 and 90 remain almost unchanged, so that only naphthalene shows a local maximum while for anthracene to pentacene a steady decrease from 0 to 90 is observed.

| Coulomb couplings
The calculation of Coulomb couplings is computationally very efficient because they are available from monomer calculations and no wavefunction overlap has to be computed. To avoid supermolecular calculations it might thus be tempting to approximate excitonic couplings as Coulomb couplings. In the following, it is addressed in which cases this approximation is numerically accurate and thus justified at the example of set A.  Figure 7 reveals that the angle has a significant impact not only on the coupling strength, c.f. Table 4, but apparently also on the wavefunction overlap, that is, the exchange contribution. A similar behavior might also be expected for displacements and tilting, as occurring in Section 3.3.
For naphthalene and pyrene, however, the supermolecular couplings exhibit local minima at 30 and local maxima at 60 , which are also observed for the reference methods, cf. Figure 5, implying that the exchange strength depends on the relative geometries in a non-trivial manner. Therefore, Coulomb couplings can be a good approximation to excitonic couplings even at small intermolecular distances, since exchange effects depend sensitively on the geometric configuration.
Indeed, this seems to be the case in the anthracene crystal, vide infra.

| Model dimers: Set B
In this section, the objective is to examine how the different position of nitrogen substitutions can influence the excitonic couplings. For this, purine and guanine as well as five different naphthyridines as shown in Figure 2 were investigated.

| Excitation energies
Analogously to set A, we start with the assessment of the accuracy of the original excitation energies obtained with the different methods.
In Tables 5 and 6

| Excitonic couplings
In  T A B L E 5 Errors in eV of the excitation energies E 1 and E 2 belonging to states S 1 and S 2 , respectively, compared to results obtained using CC2 of set B summed over all investigated distances r and angle φ = 0 to a narrow distribution. Couplings are underestimated in all cases, which is in agreement with the behavior observed in set A. SOS-CC2 results in a slightly higher ME of À15 meV and maximum deviation of 54 meV with respect to CC2 couplings.
DFT calculations employing the ωB97X functional result in a ME of +9 meV in case of distances and +4 meV in case of rotation angles.
The STD is in the same range as the semi-empirically scaled variants of CC2. For set B, the LC-DFTB method shows slightly increased deviations with respect to CC2. For example, the mean error, STD and maximum error are found to be +28 meV, 29 meV and 90 meV, respectively.
In Figure 9 the behavior of the excitonic couplings is shown for the naphthyridines using CC2 and LC-DFTB. LC-DFTB reproduces the same behavior of the couplings with respect to the distance r as CC2 for nearly all naphthyridines. When comparing the couplings calculated using CC2 in Figure 9, it becomes evident that the strength of the couplings is not strongly dependent on the position of the nitrogen atoms. Both methods show the same behavior qualitatively and quantitatively for nearly all naphthyridines.

| Dimers in a crystal: Anthracene
As pointed out, the model systems addressed so far only take into account distances and rotation angles. In real crystals, however, there are also varying relative angles of the molecular planes and horizontal translation. These influences are addressed in the following.

| Static crystal dimers
In the anthracene crystal, each molecule couples in three directions with different neighbors. The directions are denoted a, b and c, see Figure 10. In a first step, we compute supermolecular excitonic and Coulomb couplings for one selected dimer along the crystal axes a, b and c, respectively, as extracted from a crystal structure. 70 Table 8 shows Therefore it looks like, although quite surprising, that the Coulomb approximation can be a quite reasonable approach also for such densely packed molecules.
T A B L E 6 Errors in eV of the excitation energies E 1 and E 2 belonging to states S 1 and S 2 , respectively, compared to results obtained using CC2 of set B at r = 4 Å and summed over all investigated angles  In Figure 11 the coupling distributions are shown and the statistical measures can be found in Table 9. The results of the

| Influence of excitonic couplings upon diffusion coefficients
To estimate the influence of couplings on the transfer, diffusion con-

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.