Intermolecular Singlet and Triplet Exciton Transfer Integrals from Many-Body Green's Functions Theory.

A general approach to determine orientation and distance-dependent effective intermolecular exciton transfer integrals from many-body Green's functions theory is presented. On the basis of the GW approximation and the Bethe-Salpeter equation (BSE), a projection technique is employed to obtain the excitonic coupling by forming the expectation value of a supramolecular BSE Hamiltonian with electron-hole wave functions for excitations localized on two separated chromophores. Within this approach, accounting for the effects of coupling mediated by intermolecular charge transfer (CT) excitations is possible via perturbation theory or a reduction technique. Application to model configurations of pyrene dimers shows an accurate description of short-range exchange and long-range Coulomb interactions for the coupling of singlet and triplet excitons. Computational parameters, such as the choice of the exchange-correlation functional in the density-functional theory (DFT) calculations that underly the GW-BSE steps and the convergence with the number of included CT excitations, are scrutinized. Finally, an optimal strategy is derived for simulations of full large-scale morphologies by benchmarking various approximations using pairs of dicyanovinyl end-capped oligothiophenes (DCV5T), which are used as donor material in state-of-the-art organic solar cells.

Introduction

Energy transfer in molecular materials is a process in which excitation energy migrates among chromophores, e.g., a single molecule or a fragment of a larger macromolecule or molecular assembly. The control and optimization of this transfer is key to the functionality of many photoactive materials.[1,2] Natural light harvesting architectures reach near 100% efficiency converting light into a chemical form (electrical energy, molecular synthesis).[3,4] In contrast, the performance of many synthetic molecular devices is currently limited to significantly lower values. For instance, record power-conversion efficiencies of single-junction organic photovoltaic cells of around 11%[5] still trail the estimated thermodynamic maximum,[6] in part due to excitonic losses. Explicit simulations of the dynamics of excitons in complex molecular systems can be a powerful asset in revealing the interplay of electronic structure and molecular morphology on the micro- and mesoscale. Applied to a realistic representation of device materials, this knowledge will further new optimization strategies for the energy transfer processes involved. As the size of functional parts of these systems is typically on the scale of 10–100 nm, explicit ab initio treatment of the exciton dynamics is computationally intractable. Because of the disorder inherent to most organic materials as a result of material processing and thermal vibrations at room temperature, excitons tend to localize on one to a few chromophoric fragments. Localized exciton states |A⟩ can be used to define a supramolecular, tight binding like Hamiltonian aswhere ϵ is the associated excitation energy and J is the coupling between the two exciton states. J is sensitive to the molecular arrangement, i.e, the distance and relative orientation between chromophores, and constitutes a key ingredient in the calculation of exciton transfer rates. Depending on the type of excitation and length-scale on which the transfer occurs, different interactions and mechanisms mediate the coupling, as depicted schematically in Figure .

Figure 1

Illustration of the different pathways for exciton transfer between chromophores A and B. (i) Förster type energy transfer via exchange of a virtual photon. (ii) Dexter (charge) transfer via simultaneous hop of the electron–hole pair. (iii) CT mediated Dexter transfer via sequential hop of electron and hole.

In this paper, a dimer projection (DIPRO) method for exciton coupling based on many-body Green’s functions theory within the GW approximation and Bethe–Salpeter equation (GW-BSE)[7−9] is presented. Monomer electron–hole wave functions are used as pseudo diabatic states, and a projection method is employed to express these functions in the basis of products of single-particle functions used to determine the GW-BSE Hamiltonian in a supramolecular calculation. This dimer projection procedure allows for an efficient evaluation of the direct excitonic transfer integral J using linear algebra methods.

The GW-BSE approach, which is traditionally more rooted in the solid-state community, has recently received increasing attention from several groups for the treatment of electronically excited states of molecular systems.[10−18] One of the earlier key results is that it allows accurate description of localized (Frenkel) and bimolecular charge transfer (CT) excitons on an equal footing[12,19,20] owing to its accurate inclusion of long-range Coulomb and short-range exchange interactions. A similar accuracy, in particular of the CT excitations, can typically only be achieved by time-dependent density-functional theory (TDDFT) with functionals tuned for the specific system or by higher-order wave function based approaches at significant computational cost. As has been shown by several authors, using efficient localized orbital based implementations, GW-BSE can be readily applied to molecules or clusters of molecules of technological relevance.[16,17,21] Similarly, it has been realized that the use of simple self-consistent quasi-particle updates at GW level leads to improved agreement with reference energies (from experimental and/or high-level theories) as well as to a significant reduction of the notorious dependency on the DFT starting point.[22] On the basis of these promising findings, a variety of development directions are currently being pursued, aiming at increased efficiency, e.g., scaling with system size via improved models for microscopic dielectric screening,[23] overcoming current limitations, such as the treatment of transitions with multiple-excitation character,[22] or the inclusion in hybrid quantum-mechanical molecular-mechanics (QM/MM) calculations.[16−18] Specifically, here the focus is on developing a technique that allows to extend the accuracy of GW-BSE in molecule/cluster calculations to study excited state dynamics and coarsened models for supramolecular aggregates via eq .

Having an accurate representation of long-range Coulomb and short-range exchange interactions is of particular relevance for describing the wide range of possible excitonic coupling mechanisms, depending on the type of excitation and the involved length scales. Singlets can exhibit significant coupling even for distances exceeding 1 nm, typically estimated in the Förster picture from the interaction of the transition dipoles of the two chromophores involved in the transfer,[24] see Figure i. On shorter length scales, higher-order multipole terms in the Coulomb coupling[2] and short-range exchange effects can significantly influence the distance and orientation dependence of singlet couplings. Triplets couple exclusively via exchange interaction, which decays exponentially with chromophore distance and are therefore restricted to next neighbor transfers.[25] For the exchange based contributions to singlet and triplet coupling, two distinctly different pathways need to be considered. The electron–hole pair can transfer either as an entity (Figure ii) or via intermediate CT states[26] (see Figure iii). Inclusion of the CT mediated processes into an effective coupling is particularly important because they can, at short chromophore distances, contribute equally or even more than the direct process, depending on the details of the CT wave function and its energy relative to those of the localized excitons. Approximate CT wave functions have previously been constructed using, e.g., constrained DFT to study their coupling to localized excitations derived from TDDFT.[27]GW-BSE allows to derive both from the same Hamiltonian, and to this end, intermolecular CT excitations of types A–B+ and A+B– are constructed within GW-BSE-DIPRO as product states from the respective monomer single-particle orbitals. All the couplings between the CT excitations, as well as to the localized monomer excitons, are calculated. From these, an effective coupling Jeff is determined via a reduction technique, which maps the complete multistate system onto two effective states. Unlike previous approaches based on first-order perturbation theory, this reduction technique method is also applicable to cases with energetic resonances of CT and localized excitons.

In the following, GW-BSE-DIPRO is first applied to model configurations of pyrene dimers at various distances and orientations. As small-molecule prototype systems, these allow for a detailed analysis of the quality of the approach in application to bright and dark singlet, as well as triplet excitons. Particular emphasis is placed on the convergence of the results with respect to the number of included CT excitations as well as the differences between the perturbation and reduction techniques. After assessing the choice of the exchange-correlation functional in the density-functional theory (DFT) calculations that underly GW-BSE, optimizations of computational parameters are evaluated. For DCV5T-Me3, a dicyanovinyl end-capped oligothiophene used as donor material in state-of-the-art organic solar cells, this aims at defining a benchmark for the application of GW-BSE-DIPRO to large-scale morphologies of technologically relevant materials (Figure ).

Figure 2

Chemical structures of (a) pyrene and (b) DCV5T-Me3.

The paper is organized as follows: section gives a concise overview of the GW-BSE formalism, including computational details, followed by a full description of the GW-BSE-DIPRO methodology. In section , the results of applications to pyrene and DCV5T-Me3 are presented and discussed, devising an optimal strategy for application to large-scale morphologies. A brief summary concludes the paper.

Methodology

Essentials of GW-BSE

Here the essentials of the GW-BSE method with relevance for the subsequent derivation of excitonic coupling are briefly summarized. For more extensive discussions, the reader is referred to the reviews in refs (8,28).

In a first step quasi-particle (QP) states representing independent electron and hole excitations are constructed based on information obtained from the Kohn–Sham (KS) energy levels of DFT:In eq , Vext is the external potential (either of bare nuclei or pseudo atoms), VH the Hartree potential, and Vxc the exchange-correlation potential.

Within the GW approximation of many-body Green’s functions theory, as introduced by Hedin,[7] excitation energies are obtained as solutions of the quasi-particle equations:which contain the energy-dependent self-energy operator ∑(r,r′,E) in place of the exchange-correlation potential in eq . The self-energy operator is evaluated aswhereis the one-body Green’s function in quasi-particle approximation and W the dynamically screened Coulomb interaction. The latter is evaluated by first calculating the polarization P within the random-phase approximation (RPA):[29−31]Convolution of the polarization with the bare Coulomb interaction v yields the microscopic dielectric function ϵ = 1 – vP. After determination of the inverse dielectric function ϵ–1, the result is again convoluted with the Coulomb interaction to obtain W = ϵ–1v. Because of the time-consuming evaluation of the RPA as in eq , it is only performed for ω = 0 (static polarization) and a generalized plasmon-pole model[32] is used to extend the associated static dielectric function to the dynamic one.

In practical calculations, both G and W are determined using the ground state Kohn–Sham wave functions and energies. As DFT typically underestimates the fundamental gap EgKS, the self-energy and the resulting QP energies may typically deviate from self-consistent results by up to several 0.1 eV. Instead, an iterative procedure is employed in which a scissor shift Δ is applied to the Kohn–Sham spectrum before calculating W. The resulting quasi-particle gap EgQP, is determined, and from its difference to the Kohn–Sham gap, a new shift Δ=EgQP,–EgKS is defined. This procedure is repeated until convergence is reached. The QP energy levels are iterated in each step, implying according to eq updates to the Green’s function and thus the self-energy. Note that this (limited) self-consistency treatment does not change the QP structure of eq (due to satellite structures or other consequences of a self-consistent spectral shape of G(ω)).

The energetics of single-particle excitations are described with a high degree of accuracy by the quasi-particle energies.[9,14,21] Optical excitations, however, cannot be treated in such an effective one-particle picture. Instead, the coupled electron–hole excitation S can be expressed as a linear combination of products of occupied (v) and empty (c) quasi-particle functionsA is the electron–hole amplitude and can in case of singlet-to-singlet excitations be obtained by solving the generalized Bethe–Salpeter equation[8,33,34] within the Tamm–Dancoff approximation[28,35] (TDA)D = EvirtQP – EoccQP is defined via free interlevel transition energies, while K and Kd are the bare exchange and screened direct terms of the electron–hole interaction kernel, respectively. Finally, Ω is the transition energy of the optical excitation. For the case of triplet excitations, K vanishes and Ĥexc = D + Kd.

For the practical calculations in this paper, single-point Kohn–Sham calculations are performed using a modified[36] version of the Gaussian03 package,[37] Stuttgart/Dresden effective core potentials,[38] and the associated basis sets that are augmented by additional polarization functions[39] of d symmetry. For all steps of the actual GW-BSE calculations, the Gaussian-type orbital based implementation in the VOTCA package is employed.[12,15,40,41] It is specifically optimized for application to molecular systems by using an auxiliary basis set to express the quantities occurring in the GW self-energy operator and the electron–hole interaction in the BSE. We include orbitals of s, p, and d symmetry with the decay constants α (in a.u.) 0.20, 0.67, and 3.0 for N and S, 0.25, 0.90, 3.0 for C, and 0.4 and 1.5 for H atoms, yielding converged excitation energies. Further technical details can be found in refs (10,42).

Excitonic Coupling Elements via GW-BSE-DIPRO

In the following, diabatic states |A⟩ and |B⟩ are approximated by monomer electron–hole wave functions, |Φ⟩ and |Φ⟩, as defined in eq , because an exact diabatization is for most systems difficult or even impossible.[43] With the BSE Hamiltonian of the dimer formed by the chromophores, ĤD, one can setup an effective (2 × 2) generalized eigenvalue problemwithand is the ith eigenvector and ϵ the corresponding eigenvalue. Because the normalized monomer states are only approximations to the diabatic states, they are not necessarily orthogonal, ≠ , and J cannot be directly identified with the off-diagonal elements of . Instead, following the idea successfully established by Valeev et al. for electronic coupling,[44,45] the generalized eigenvalue problem in eq first needs to be transformed into a standard eigenvalue problemvia a Löwdin orthogonalization.[46] The choice for this technique is motivated by the fact that unlike, e.g., the Schmidt orthogonalization, it treats both wave functions on an equal footing and the resulting symmetrically orthogonalized functions are the least distant in the Hilbert space from the original functions. With the orthonormalized states |Φ̃⟩, the diagonal and off-diagonal elements of can be identified with the excitation energies ϵ and exciton coupling elements J of eq , respectively. It should be noted that this construction of approximate diabatic dimer states starting from monomer functions as above can be considered complementary to the approximate diabatization of adiabatic dimer states via the Boys or Mulliken–Hush scheme.[47−49] While such an approach is also feasible, it comes at a higher computational cost, in particular due to the fact that solutions to eq for the dimer system are required. This time-consuming step for dimers of larger molecules can be avoided using the approximation adopted here, as will be outlined in the following.

What remains is to calculate the elements of and . They can be written as ⟨Φ|Ô|Φ⟩ with i ∈ A,B and Ô ={Ĥ,1}. In the following, |v⟩, |α⟩, |β⟩ (|c⟩, |α′⟩, |β′⟩) are occupied (empty) single-particle orbitals of the dimer, monomer A, and monomer B, respectively. With this, the electron–hole wave functions for the localized monomer excitations (eq ) can be written as

By inserting the identity I = ∑|vc⟩⟨vc| twice into the expression for the off-diagonal matrix elements of eq , and using the above definitions of |Φ⟩ and |Φ⟩ from eqs and 14, we obtain

In practical calculations, Ô = ĤD is setup directly in terms of |vc⟩, so ⟨vc|ĤD|v′c′⟩ = HBSE is readily available. For Ô = 1, it holds that ⟨vc|v′c′⟩ = δδ, yielding

The quantities κ are projections of the monomer electron–hole wave functions, expressed in monomer single-particles functions, onto dimer single-particle orbitals, e.g.,

These projections are evaluated by inserting the expansion of the respective single-particle orbitals in terms of the atomic orbital basis. When dimer and monomer calculations share the same basis set of atomic functions {|χμ⟩}, it holds thatThus, the terms of type ⟨αα′|vc⟩ occurring in eq can be rewritten aswhere is the overlap matrix of the atomic orbitals.

Influence of Intermolecular CT States

While the projection technique as presented above captures the coupling mechanisms depicted in Figure i,ii), the charge transfer state mediated mechanism (iii) is not accounted for. The aim is to include these effects, which are CT inherent and whose interactions are contained in the BSE Hamiltonian but which are not represented by the subspace spanned by |Φ̃⟩.

Intermolecular CT excitons are approximated as product states of two single-particle orbitals localized on different monomers, i.e.,where αH (βH) is the highest occupied molecular orbital (HOMO) and αL (βL) is the lowest unoccupied molecular orbital (LUMO) of chromophore A(B), respectively, and i,j = 0,...,M. M is the number of additional orbitals belove (above) the HOMO (LUMO) taken into account.[50] In total, a set {|CT⟩} comprising NCT = 2(M + 1)2 CT excitations is generated according to eq .

As these approximate CT states are not orthogonal to the orthonormalized localized states after eq , each |CT⟩ is first individually orthogonalized with respect to |Φ̃⟩and then normalized via .

Eq then turns into a ([2 + NCT] × [2 + NCT]) eigenvalue problem with block structure of the augmented Hamiltonian:where FE = . In a final step, the subspace of CT states in eq is diagonalized, i.e., solvingThe eigenfunctions and energies ΩCT are used to transform eq into an ordinary eigenvalue problem:with CT = diag(ΩCT). For the special case of M = 0, i.e., construction of two CT like excitations from the respective HOMO and LUMO single-particle orbitals, this corresponds to the (2 + 2) × (2 + 2) system

Perturbation Theory

To obtain an effective excitonic coupling element Jeff between |Φ̃⟩ and |Φ̃⟩ that includes effects from coupling via intermediate CT excitations, the influence of the latter on the localized states has to be evaluated. Within first-order perturbation theory, the corrections |δΦ̃⟩ to |Φ̃⟩ due to the |CT̃⟩ can be expressed as[26,51]The modified coupling is then obtained to first-order in |δΦ̃⟩ as

All terms required to evaluate eq can be identified with elements of the Hamiltonian in eq : ⟨Φ̃|Ĥ|Φ̃⟩ is the off-diagonal element of FE (i.e., the unperturbed excitonic coupling), ⟨Φ̃|Ĥ|CT̃⟩ are elements of FE–CT, and the energies occurring in the denominator are the diagonal elements of FE and CT, respectively. For the example of the M = 0 case as in eq , the expression for the effective excitonic coupling element explicitly reads:

Reduction Method

From the structure of eq , it is apparent that the perturbative approach to account for the influence of CT excitations on excitonic coupling is not readily applicable to situations in which a CT excitation is energetically in, or close to, resonance with the localized excitations. Instead of going to even higher orders of perturbation theory, an alternative technique that starts from the augmented Hamiltonian of eq is proposed.

The main idea is to reduce the augmented (2 + NCT) × (2 + NCT) system to an effective (2 × 2) system. In spirit similar to perturbation theory, the states forming this reduced system are expected to be close to the original states |Φ̃⟩, and consequently the effects of the intermediate CT states is mapped onto a coupling between those states. To achieve this, first eq is diagonalized, yielding the eigenenergies ϵ̃ and the set of corresponding eigenvectors .

From this, two elements C̃ are chosen according to having the maxmium overlap with the states Φ̃ and Φ̃, respectively. Projecting C̃ onto the subspace spanned by Φ̃ and Φ̃, followed by a Löwdin transformation, yields new orthonormalized vectors *.

The diagonal (2 × 2) matrix ϵ* formed with the energies ϵ̃ can be transformed to its nondiagonal form using the transformation matrix = (**). Resulting is a reduced, effective systemwhich allows to read off the effective excitonic coupling Jeff as its diagonal elements.

To illustrate the differences and similarities between obtaining the effective coupling according to this reduction method (RM) and the perturbation theory (PT), it is convenient to consider a simplified model of the minimal system introduced in eq . Specifically, a symmetric system is assumed with ϵ = ϵ = ϵ, Ω1 = Ω1 = Ω, and J=JCT. Using perturbation theory, the effective coupling reads with Δϵ = Ω – ϵwith the obvious resonance for Δϵ = 0. Using the reduction method yields an analytical solution

In the limit of Δϵ → 0, Jeff,RM remains finite. Away from the energetic resonance, i.e. Δϵ ≫J,J, it holds that

For more complex systems with less symmetry or M > 0, no closed form analytical expressions for Jeff,RM can be obtained. Therefore, the method is in the following employed and assessed in practical application to realistic molecular systems.

Results

To assess the quality of the procedures outlined in the previous section, model configurations of pyrene dimers at various distances and orientations are considered. Within the TDA of GW-BSE, pyrene exhibits energetically well separated optically inactive (S1) and active (S2) singlet as well as triplet (T1) excitations. Analysis of the excitonic couplings for these different type of excitations allows scrutinization of how well the different pathways, see Figure , are accounted for. The convergence of the results with respect to the number of included CT excitations, the differences between the perturbation and reduction techniques, basisset dependence, as well as the influence of the choice of the exchange-correlation functional in the density-functional theory (DFT) calculations that underly the GW-BSE steps are evaluated. Later, optimizations of computational parameters are devised using DCV5T, defining a benchmark for the application of GW-BSE-DIPRO to large-scale morphologies of technologically relevant materials.

Model Pyrene Dimers

The geometry of a single pyrene molecule was optimized on DFT level using the PBE functional[52] with the 6-311G(d,p) basis set. From the this geometry, ideal π-stacked dimers with intermolecular distances ranging from 2.7 to 7 Å are constructed. GW-BSE calculations are performed for the monomers, as well as the dimer configuration (only setup of BSE), and coupling elements determined according to the projection method.

In this configuration in which the molecules forming the dimer are related by a symmetry transformation and the energetic states of the monomers are well separated, it is also possible to obtain the effective excitonic coupling via 2Jeff = ΔΩ, where ΔΩ is the Davydov splitting of the respective monomer excitation in the dimer. To facilitate this comparison, we also perform a full GW-BSE calculation for the dimer configuration and extract the splitting from the resulting spectrum.

Figure shows the distance dependence of |J| for (a) S1, (b) S2, and (c) T1 excitations obtained via GW-BSE-DIPRO with the reduction method. For all excitations, large deviations from the Davydov splitting are found when no intermediate CT states (NCT = 0) are included. At the typical π–π stacking distance of 3.5 Å, these deviations can be on the order of 1–2 orders of magnitude. For the two singlet states, the observed underestimation decreases exponentially with distance, typical for an exchange based coupling mechanism as the one mediated by charge transfer. Inclusion of CT excitations ameliorates this situation. For S2 and T1, taking only two CT excitations between the respective HOMO and LUMO states (M = 0, NCT = 2) into account practically recovers the split results. In contrast, no effect can be registered for S1. Here, an agreement with the Davydov split estimate can only be achieved for M = 1, i.e., by construction of additional CT excitations based on HOMO–1 and LUMO+1, respectively. This is due to the fact that, unlike S2 and T1, which have the main contribution from a HOMO → LUMO transition, the first singlet excitation in single pyrene is formed by a linear combination of HOMO–1 → LUMO and HOMO → LUMO+1, and the choice of M in the reduction method needs to reflect the composition of the various localized excited states. In all cases, converged exciton couplings are achieved including NCT = 50 CT excitations (M = 4). Note that for distances ≥6 Å, the coupling in S1 and T1 becomes so small that the split estimate becomes numerically inaccurate.

Figure 3

Distance dependence of excitonic couplings for (a) S1, (b) S2, and (c) T1 excitations in an ideally π-stacked pyrene dimer. Results obtained via GW-BSE-DIPRO with the reduction method for increasing numbers of included CT excitations are compared to the reference determined from the Davydov splitting in a full supermolecular calculation.

From the converged results and the Davydov splittings, it can be seen that the excitonic couplings based on GW-BSE simultaneously exhibit characteristics of short-range exchange and long-range Coulomb coupling depending on the type of excitation. For S1 (with a negligibly small transition dipole) and T1, |J| decays proportional to exp(−αd) with α(S1) ≈ 3.2 Å–1 and α(T1) ≈ 3.4 Å–1, respectively. In contrast, the optically active S2 shows an exponential decay for distances in the range of 3–4.5 Å (α(S2) ≈ 1.7 Å–1), before the effective coupling is dominated by slowly decaying Coulomb contributions.

To ascertain the quality of excitonic coupling elements obtained from GW-BSE, they are in the following compared to ones obtained from standard methods of similar complexity: time-dependent Hartree–Fock (TDHF), TDDFT/B3LYP, and configuration-interaction singles (CIS). Because the projection technique as used in GW-BSE-DIPRO is not available for those, the comparison is performed for the Davydov splittings in the ideally π-stacked dimer configurations, using the same ECP and basis set. The respective results are listed in Table for dimer separations of d = 3.0–6.0 Å. At all separations, GW-BSE, TDHF, and CIS agree well for the dark S1 state. Same holds for the optically active S2. In fact, here, the data from all four methods practically agree at short and long distances, while in the intermediate region TDDFT appears to decay slightly faster. All in all, a good agreement between all four methods can be noted for singlets. For the T1 state all methods show approximately the same decay of the coupling with separation. Note, however, that TDHF calculations suffered from triplet instabilities, yielding negative excitation energies, which led to the results being discarded.[53,54]

Table 1

Excitonic Coupling Elements |J| from Davydov Splitting (in eV) for S1, S2, and T1 Excitations in an Ideally π-Stacked Pyrene Dimer with Varying Distance d, as Obtained by GW-BSE, TDHF, TDDFT/B3LYP, and CIS, Respectivelya

type	d = 3.0 Å	d = 4.0 Å	d = 5.0 Å	d = 6.0 Å
Davydov Split S1
GW-BSE	6.13 × 10–1	3.70 × 10–2	1.53 × 10–3
TDHF	6.24 × 10–1	4.94 × 10–2	2.50 × 10–3
TDDFT	6.81 × 10–1	1.08 × 10–1	8.35 × 10–3
CIS	6.41 × 10–1	4.77 × 10–2	2.50 × 10–3
Davydov Split S2
GW-BSE	9.57 × 10–1	1.71 × 10–1	4.80 × 10–2	2.18 × 10–2
TDHF	8.38 × 10–1	1.40 × 10–1	5.11 × 10–2	3.22 × 10–2
TDDFT	6.74 × 10–1	7.36 × 10–2	3.45 × 10–2	2.87 × 10–2
CIS	8.87 × 10–1	1.55 × 10–1	5.60 × 10–2	3.56 × 10–2
Davydov Split T1
GW-BSE	5.12 × 10–1	2.54 × 10–2	8.51 × 10–4
TDDFT	5.26 × 10–1	3.83 × 10–2	1.55 × 10–3
CIS	3.97 × 10–1	2.33 × 10–2	1.70 × 10–3

For d = 6.0 Å, the coupling in S1 and T1 becomes so small that the split estimate becomes numerically inaccurate and is therefore omitted.

In many realistic molecular aggregates, chromophores do not arrange in an ideal π-stack as assumed in the previous section. Instead, they assume relative positions and orientations characterized by shifts and rotations which are not compatible with the basic symmetry operations. Because of the asymmetry in the geometry, an estimation of coupling elements from Davydov splits is inaccurate and the use of techniques such as GW-BSE-DIPRO is indispensable. To assess the procedure for such a case, starting from the ideal π-stacking configuration with an intermolecular distance of 6.5 Å, one molecule is rotated along its long axis from 0° to 90° (see, cf., inset in Figure b). As for the distance dependence, the convergence of the excitonic coupling elements with respect to the number of included CT states is investigated for S1, S2, and T1. The results shown in Figure generally exhibit more structure compared to the distance dependence in Figure as a consequence of intricate interactions between the two π systems upon rotation. It is also evident that for rotation angles of up to approximately 50°, converged results are obtained for NCT = 50 for all excitations. At larger rotations, strong couplings are found for S1 and T1 in particular. In this region, the convergence is much slower and up to 2450 (M = 34) intermediate CT states are required. This is probably a result of the stronger and asymmetric polarization of the dimer states with respect to the monomer calculations. It should be emphasized, however, that the two chromophores approach each other very closely for those angles. At the perpendicular configuration, the minimal distance is reduced to only 3.1 Å and concomitantly strong effects and mixing of single-particle functions can be expected.

Figure 4

Rotational dependence of excitonic couplings for (a) S1, (b) S2, and (c) T1 excitations in a pyrene dimer. The configuration at ϕ = 0° corresponds to and ideal π-stacking at a distance of 6.5 Å. Results obtained via GW-BSE-DIPRO with the reduction method for increasing numbers of included CT excitations. Molecules have been rendered with VMD.[55]

To highlight the differences and similarities, distance and rotational dependence of singlet and triplet excitonic transfer integrals as obtained by the reduction method presented in this work and the perturbation theory are compared for the model pyrene configurations. In all cases, 50 CT states are included in the distance dependence and 2450 in the rotational dependence. As is apparent from Figure a, the distance dependence of the coupling constants for S1 and T1 shows good agreement between the approaches. Only at distances smaller than 3.5 Å perturbation theory slightly overestimates the RM results, which agree with the Davydov split estimates, cf. Figure . For the optically active S2 state, however, significant deviations can be observed. For an intermolecular separation of 3.5 Å, a characteristic resonance structure is found, representing a massive overestimation of the transfer integral by nearly 2 orders of magnitude. Deviations are noticeable around this typical π-stacking distance in many molecular semiconductors up to a distance of 4.5 Å. For the rotated systems, see Figure b; both approaches yield qualitatively similar behavior with some quantitative deviations up to a factor of 2 for the close contact structures at large rotation angles. All in all, the reduction method compares favorably with the perturbation theory approach.

Figure 5

Comparison of the (a) distance and (b) rotation dependence of effective excitonic couplings in the pyrene model dimers, obtained with reduction method (RM) and via perturbation theory (PT), respectively. The distance dependence of both approaches as seen in (a) is nearly identical for S1 and T1 with the exception of short intermolecular distances. For the S2 couplings, large deviations are observed due to energetic resonance of localized and intermediate CT excitons. For the rotated structures (b) both approaches show similar qualitative behavior.

As mentioned in section , the basis set used in this work is the based on the one optimized for the use with the effective core potentials,[38] which will be referred to as bsECP in the following. All results presented to this point have been obtained from calculations, in which (d,p) polarization functions of the 6-311G basis set[39] have been added to form the basis bsECP(d,p). To further gauge the convergence of the reported values for the singlet and triplet exciton coupling elements with respect to the basis set choice, two basis sets with diffuse functions (bsECP(d,p)+, s shells, decay constant 0.0438 au for C, and 0.102741 au for H) and (bsECP(d,p)++, additional p shell, decay constant 0.0691 au for C) have been prepared as well. Three representative configurations of the pyrene dimers have been chosen: ideal π-stack geometries with distances 3.5 and 7.0 Å, respectively, and one at d = 5.5 Å with an additional rotation by 45°. For the latter, no results could be obtained using the bsECP(d,p)++ because of convergence issues in the underlying DFT calculation. The obtained couplings for S1, S2, and T1 excitons for these configurations are listed in Table . It is clear that the results are fairly independent of the choice of the basis set. Even the use of the standard bsECP without additional polarization or diffuse functions yields excitonic transfer integrals in reasonable agreement with the ones obtained by extended basis sets. This is of great significance in terms of the computational costs listed at the bottom of Table , from which it is obvious that the addition of diffuse functions in particular increases the calculation times dramatically. It should also be noted that this scaling of the computation time depending on the choice of basis set is exacerbated for molecules bigger than pyrene and is an important aspect to control in application to large-scale morphologies, as will be further discussed below.

Table 2

Basis Set Dependence of the Calculated Transfer Integrals |J| (in eV) for S1, S2, and T1 States in Representative Configurations of Pyrene Dimersa

type	bsECP	bsECP(d,p)	bsECP(d,p)+	bsECP(d,p)++
Ideal π-Stack, d = 3.5 Å
S1	1.90 × 10–1	1.85 × 10–1	1.87 × 10–1	1.95 × 10–1
S2	4.59 × 10–1	4.20 × 10–1	4.31 × 10–1	4.53 × 10–1
T1	1.42 × 10–1	1.37 × 10–1	1.39 × 10–1	1.46 × 10–1
Ideal π-Stack, d = 7.0 Å
S1	1.86 × 10–5	1.35 × 10–5	1.33 × 10–5	1.32 × 10–5
S2	2.52 × 10–2	2.16 × 10–2	2.19 × 10–2	2.29 × 10–2
T1	4.18 × 10–7	3.08 × 10–7	2.47 × 10–7	1.62 × 10–7
d = 5.5 Å, Rotation ϕ = 45°
S1	3.44 × 10–4	4.18 × 10–4	4.31 × 10–4
S2	5.91 × 10–2	4.88 × 10–2	5.06 × 10–2
T1	3.97 × 10–5	5.04 × 10–5	6.63 × 10–5
timings	1.00	1.26	1.50	2.20

The average timings for the different calculations relative to the one with the smallest bsECP set are given at the bottom.

Until this point, the discussion of the GW-BSE-DIPRO approach has been limited to results obtained via GW-BSE based on DFT calculations using the semilocal PBE functional. It is known from literature that this technique can be marred by a starting point dependence, i.e., that the computed excited states depend on the quality of the underlying ground-state calculation.[56] In particular, one-shot G0W0 techniques are affected by this. As a result of the iterative procedures regarding the refinement of quasi-particle energies as described in section for the GW-BSE implementation used here, this problem is alleviated. To confirm this, the distance dependence of the exciton transfer integrals has been recalculated using DFT with the B3LYP hybrid functional[57,58] for the pyrene dimer. The comparison to the results based on PBE as shown in Figure a reveals practical independence on the DFT starting point for all three types of excitations considered.

Figure 6

Comparison of the distance dependence of exciton transfer integrals in ideally π-stacked dimers of (a) pyrene and (b) DCV5T, obtained starting from DFT calculations using PBE and B3LYP functionals, respectively. Fifty CT states have been taken into account using the reduction method.

Optimizations for Application to Large Scale Morphologies

In the following, the focus shifts to the evaluation of GW-BSE-DIPRO to chromophores of sizes typical and relevant for application in organic devices. To this end, exciton transfer integrals in an ideal stack of DCV5T molecules are investigated for the lowest energy singlet (S1, optically active) and triplet excitons (T1), respectively.

The obtained distance dependence is shown in Figure b, with the GW-BSE calculations being based on PBE and B3LYP DFT ground states. As for the pyrene dimer, different behavior is observed for the two excitations. The coupling of the optically active singlet states is dominated by long-range Coulomb interactions, representing Förster type coupling. In contrast, T1 shows an exponential distance dependence, as expected for a Dexter type exchange coupling. Results for the two functionals are practically identical, with a shift to slightly larger S1 couplings noticeable for PBE compared to B3LYP. This is a consequence of an approximately 9% bigger transition dipole moment resulting from the former.

The above confirms that the GW-BSE-DIPRO method used in this work is well applicable to complex molecular systems of relevant size. However, the investigation of, e.g., exciton diffusion in realistic large-scale morphologies requires the calculation of tens of thousands excitonic transfer integrals. It is therefore highly desirable to devise optimizations of the involved computational procedures. To this end, approximations on DFT, GW-BSE, and GW-BSE-DIPRO levels to decrease computation times are evaluated.

On the DFT level, it has been shown for the calculation of charge transfer integrals that it is computationally advantageous to use an initial guess for the dimer calculation formed by merging the involved densities of the monomer fragments.[45] This also allows to perform only a single SCF step on the dimer instead of obtaining a fully self-consistent solution.

A possible simplification of the GW-BSE run concerns the iterative procedure used to scissor shift the Kohn–Sham spectrum before calculating W. Instead of iterating this shift for each dimer configuration, a fixed value can be predetermined, e.g., from a single representative configuration or averaging over a couple.

In Table , the run times and exciton transfer integrals for S1 and T1 in a DCV5T dimer separated by 3.5 Å are compared for different combinations of computational procedures. Using a fixed shift in GW-BSE saves on average 3–4 iterations, which for this system translates to reducing the runtime by approximately 4 min. The effect on the exciton couplings is negligible. Similarly, using the noSCF procedure in the DFT part has a slightly bigger effect decreasing the singlet coupling by about 1.5%. Combining both approximations cuts computation time by about 25%.

Table 3

Effect of Different Computational Parameters for DFT and GW-BSE Calculations on Run Times and Exciton Transfer Integrals for S1 and T1 in a DCV5T Dimer Separated by 3.5 Åa

DFT@PBE	GW-BSE	time [min:s]	|J|S1 [eV]	|J|T1 [eV]
SCF	iterate	44:50	0.5354	0.1709
SCF	fixed	40:28	0.5355	0.1710
noSCF	iterate	37:51	0.5277	0.1721
noSCF	fixed	33:40	0.5278	0.1721

Calculations were performed using four threads on a i5-4690 CPU @3.50 GHz. The value of 0.22 Ryd for the fixed shift runs was taken from the result of the iterative procedure at 3 Å.

Additionally, the distance dependence of triplet and singlet couplings in DCV5T has been calculated for all four options given in Table . The results shown in Figure underline that the choice of noSCF and fixed shift does not lead to larger deviations even at greater intermolecular separations. Even though this choice reduces the computational time by 25% compared to using fully self-consistent procedures on all levels, the remaining absolute time is still substantial. In particular, the slow distance decay of couplings for optically active singlets implies the necessity to determine transfer integrals for a number of chromophore dimers that is intractable for any first-principles based technique such as GW-BSE-DIPRO.

Figure 7

Comparison of the effect of different approximations on the distance dependence of triplet (a) and singlet (b) couplings in DCV5T. Inset in (b) shows additionally estimates from transition partial charges and transition dipole interactions for very long distances.

The setup of the dimer GW-BSE Hamiltonian matrix takes up a significant amount of the total computer time when applied to larger molecules. As stated in section , we have so far considered all nocc occupied and as many virtual single-particle levels in the construction of the product basis for the dimer. With nocc = 186 for DCV5T, this corresponds to a matrix of dimension 34569. Naturally, calculations can be accelerated by reducing the number of single particle states taken into account in this step. Using several values of the factor f to set the range single particle states considered to [nocc(1 – f):nocc(1 + f)], the relative reduction in computational time and relative deviation of singlet couplings is evaluated. As depicted in Figure , reducing f to 0.5 approximately halves the computational cost, while only a 8% decrease in accuracy is observed.

Figure 8

Relative reduction of computation time of the GW-BSE calculation (excluding the DFT part) for the DCV5T singlet coupling. Couplings were calculated for five distances and the relative errors with respect to the high dimension BSE result averaged.

Furthermore, it can be assumed that beyond some distance, the long-ranged couplings are given as classical Coulomb interactions of transitions densities of the constituent chromophores. It is convenient to map the full transition density to a set of atomic partial charges that reproduces its electrostatic potential.[59,60] Such a classical model allows the computation of many thousands excitonic transfer integrals per minute. As can be seen in Figure b, the results from this TrESP approach agree well with the GW-BSE-DIPRO couplings beyond a separation of 5 Å, a typical distance at which intermolecular exchange effects can be expected to be negligible relative to the Coulomb interactions of large transition densities. The inset shows a comparison between classical interactions of transition charges with those of transition dipoles (TrDip), a coarser and often employed representation of transition densities. At a distance of 5 Å, the dipole coupling overestimates the TrESP and GW-BSE-DIPRO results by 1 order of magnitude. Practical agreement can only be observed for separations larger than 60 Å, casting doubts about the use of the TrDip approximation for intermediate distances.

Summary

In this paper, a general approach to determine orientation and distance-dependent effective intermolecular exciton transfer integrals from many-body Green’s functions theory within the GW approximation and the Bethe–Salpeter equation (BSE) has been derived. A projection technique is employed to obtain the excitonic coupling by forming the expectation value of a supramolecular BSE Hamiltonian with electron–hole wave functions for excitations localized on two separated chromophores. Within this approach effects of coupling mediated by intermolecular charge transfer (CT) excitations are accounted for via a reduction technique that proves to be applicable to situation in which conventional perturbative approaches break down.

Application to model dimers reveals an accurate description of short-range exchange and long-range Coulomb interactions for the coupling of singlet and triplet excitons by this GW-BSE-DIPRO technique. An optimal strategy for simulations of full large-scale morphologies includes a combination of loosening of self-consistency parameters, reduction of the active space in GW-BSE, and, for optically active singlets, a change to classical transition density-based interaction models at larger distances. All developed methods are available within VOTCA software package.[41]