The accurate prediction of singlet and triplet excitation energies is an area of intense research of significant fundamental interest and critical for many applications. Most calculations of singlet and triplet energies use time-dependent density functional theory (TDDFT) in conjunction with an approximate exchange-correlation functional. In this work, we examine and critically assess an alternative method for predicting low-lying neutral excitations with similar computational cost, the ab initio Bethe-Salpeter equation (BSE) approach, and compare results against high-accuracy wavefunction-based methods. We consider singlet and triplet excitations of 27 prototypical organic molecules, including members of Thiel’s set, the acene series, and several aromatic hydrocarbons exhibiting charge-transfer-like excitations. Analogous to its impact in TDDFT, we find that the Tamm-Dancoff approximation (TDA) overcomes triplet instabilities in the BSE approach, improving both triplet and singlet energetics relative to higher level theories. Finally, we find that BSE-TDA calculations built on effective DFT starting points, such as those utilizing optimally tuned range-separated hybrid functionals, can yield accurate singlet and triplet excitation energies for gas-phase organic molecules.

The quantitative prediction and understanding of low-lying excitations in organic molecules are of significant fundamental interest and technological relevance. For example, a better understanding of multiexciton phenomena in organic molecular systems such as singlet fission (SF),1–3 a process by which a singlet exciton decays into two low-energy triplet excitations, can lead to external quantum efficiencies above 100%2,3 and is therefore desirable for next-generation solar cells and other optoelectronic applications. Such multiexciton energy conversion phenomena are dependent on a subtle balance between singlet and triplet excitation energies, and predictions of such energetics call for accurate ab initio methods.

A widely used ab initio formalism for neutral excitations is time-dependent density-functional theory (TDDFT). For gas-phase acene molecules, the performance of TDDFT with a number of exchange-correlation functionals is well-documented: overall, TDDFT with standard functionals—e.g., local, semilocal, and global hybrid exchange-correlation functionals—fails to predict triplet excitations4,5 by 0.4–1.8 eV in acenes, as well as the ordering and absolute energies of the two lowest-lying singlets,6,7 one of which has charge-transfer-like character7 (as detailed in Section III). These failures have been ascribed to (i) the so-called “low orbital overlap problem” in global hybrid functionals, in which the overlap between spatially separated molecular orbitals is usually overestimated and to (ii) triplet instabilities associated with TDDFT using standard approximate exchange-correlation functionals.5,7–10

Beyond TDDFT approaches with conventional functionals, range-separated hybrid (RSH) functionals have been shown to mitigate the low orbital overlap problem.5,7,11,12 In this class of functionals, the Coulomb potential is partitioned into short- and long-range contributions, with the important consequence that different fractions of exact exchange can be used in the short and long range.13,14 This partitioning is usually implemented as

1r=α+βerf(γr)r+1[α+βerf(γr)]r,
(1)

where the first term is treated explicitly, and the second is replaced with a semi-local functional, such as one of several generalized gradient approximations (GGAs). The α, β, and γ parameters are either fixed as in, e.g., CAM-B3LYP,14 or tuned to fulfill DFT theorems as in optimally tuned range-separated hybrid (OTRSH) functionals15 or Koopmans compliant functionals.16 Two examples of OTRSH functionals are that of Baer-Neuhauser-Livshits (BNL)17,18 and the Perdew-Burke-Ernzerhof (PBE)-based OTRSHs.15,19,20 Importantly, CAM-B3LYP, OTRSH-BNL, and other RSH functionals have proven quite successful in predicting the low-lying excitations of aromatic hydrocarbons5,7,8,11,12 and charge-transfer (CT) excitations.21–23 

An alternative approach to neutral excitations is ab initio many-body perturbation theory (MBPT), where neutral excitations are computed via two-particle Green’s functions through solution of an effective two-particle equation, the Bethe-Salpeter equation (BSE). In MBPT, solutions to the BSE build on one-particle energies and wavefunctions, usually obtained from a generalized Kohn-Sham DFT starting point within GW approximation, where G is the one-particle Green’s function and W is the screened Coulomb interaction. Referred to as the GW-BSE approach hereafter, this method has been extremely successful for solids,24 as it goes beyond DFT by including electron-hole interactions, significant for molecules and other systems. It has also been successfully applied to gas-phase molecules,25,26 and quantitative benchmark studies are beginning to appear.27–31 Yet much remains unknown about the performance of ab initio GW-BSE calculations, particularly their ability to predict acene excitations, charge-transfer-like excitations of aromatic molecules, and, more generally, the triplet excitations of organic compounds. The aim of the present work is to address these issues.

The Bethe-Salpeter equation is a formal solution to the two-particle Green’s function, giving access to excitonic wavefunctions and eigenvalues. The underlying theory and approach are explained in more detail in Refs. 24, 32, and 33. In finite systems with real wavefunctions, BSE is exactly analogous to TDDFT. As in the Casida equations of linear-response TDDFT, solutions to the BSE can take the form of an eigenvalue problem, i.e.,

(ABBA)(XsYs)=Ωs(XsYs),
(2)

where Xs and Ys are the excitonic wavefunctions, Ωs are the eigenvalues, and the A and B blocks form the resonant and coupling parts of the BSE, respectively.32 In the Tamm-Dancoff approximation (TDA),34 the B and −B blocks are neglected, resulting in the decoupling of the A and −A blocks. While the applicability and implications of the TDA in TDDFT are well documented,4,5,9,10,35,36 the quantitative impact of the TDA on GW-BSE calculations of small molecules remains underexplored.

Benchmarks against wavefunction-based theory are the norm to assess the accuracy of approximations within lower order theories or approximations to TDDFT. For example, Thiel’s set37 contains reference-quality excitation energies which were obtained with multistate multiconfigurational second-order perturbation theory such as (MS-CASPT2) and various coupled-cluster theories, such as coupled cluster with singles, doubles, and perturbative triples (CCSD(T)).37,38 Recent work28,29,31 has explored the performance of GW-BSE for Thiel’s set, reporting that GW-BSE can indeed yield quantitative singlet energetics under several approximations. However, as in TDDFT, GW-BSE-calculated triplets were found to be substantially underestimated.28,31 The limited performance of GW-BSE in this case merits further analysis.

Motivated by the success of GW-BSE in general, and its reduced computation cost relative to wavefunction methods, herein we assess the performance of different approximations to GW-BSE and determine successful approaches within this framework for the quantitative prediction of low-lying excitations of organic compounds. We evaluate GW-BSE against multireference and coupled cluster references for representative singlet and triplet excitations of 27 organic molecules, including hydrocarbons, heterocycles, aldehydes, ketones, and amides (see Figure 1). We focus on approximations to GW that enter the BSE, including hybrid functional-starting points with one-shot schemes and the effect of partially self-consistent schemes, as described in Section II. We also provide a detailed assessment of the performance of the BSE and the TDA relative to that of other two-particle Green’s function approaches and computationally less-expensive TDDFT methods.

FIG. 1.

Top: Subset of 20 organic molecules containing triplet excitations from Thiel’s set. Bottom: The general formula for an acene molecule and the three other aromatic hydrocarbons studied here: azulene, benzo[e]pyrene (BP), and dibenzo[a, c]anthracene (DBAn). H is white, C is light blue, N is dark blue, and O is red.

FIG. 1.

Top: Subset of 20 organic molecules containing triplet excitations from Thiel’s set. Bottom: The general formula for an acene molecule and the three other aromatic hydrocarbons studied here: azulene, benzo[e]pyrene (BP), and dibenzo[a, c]anthracene (DBAn). H is white, C is light blue, N is dark blue, and O is red.

Close modal

Our calculations start with a self-consistent time-independent DFT calculation using an approximate exchange-correlation functional (see below). For the molecules considered here, we minimize the total energy with respect to the density using fixed atomic coordinates for all molecules obtained from Ref. 37 (see the supplementary material for more details). Starting from the output of our DFT calculations, and using the MolGW package,39,40 we then compute one- and two-particle excitation energies with the GW and GW-BSE approaches, respectively. As is standard, our GW-BSE calculations build on single-particle states, which are coupled in the two-particle BSE equation via the electron-hole interaction kernel. With GW input, BSE is recast into an eigenvalue problem,32 the solution of which yields the energies and eigenstates of a set of neutral excitations. As detailed in prior work by us and others,41–45GW calculations are sensitive to the generalized Kohn-Sham starting point and to whether self-consistency is used. Here, we build on previous work and use three accurate GW schemes: G0W0@BHLYP (one-shot GW on top of BHLYP,46 which has 50% exact-exchange); G0W0@OTRSH-PBE;47 and eigenvalue self-consistent GW (evGW), in which the quasiparticle energies are updated (in both G and the polarizability) one or more times prior to calculating the final self-energy corrections.41 

As mentioned, our GW-BSE calculations are performed with the MolGW package, in which the frequency dependence of the GW non-local self-energy Σ(𝐫,𝐫,ω) is treated analytically, and hence is exact for a given basis set, without the need for plasmon-pole approximations. We use conventional approximations to solve the BSE: irreducible vertices are set to 1, the polarizability and other matrix elements are constructed using GW eigenvalues and DFT wavefunctions, the screened Coulomb interaction is evaluated in the random phase approximation (RPA), and a static electron-hole screening is used; see Ref. 24. We adopt the aug-cc-pVTZ basis set48 which ensures convergence better than 0.1 eV for the excitation energies shown here (see the supplementary material for details). In order to reduce the computational load, and for the purpose of parallelization, we use the resolution-of-the-identity in the Coulomb metric,49,50 as implemented in MolGW,39,45 with the well-established auxiliary basis sets of Weigend50 which are consistent with Dunning basis sets. The resolution-of-the-identity is expected to have a small effect on the GW energies, on the order of 1 meV, as we have demonstrated in the case of benzene.45 

For OTRSH-PBE, as a standard procedure for the acenes,19,20,45 we set α=0.00.2 (see the supplementary material), which fixes the amount of short-range Fock exchange to 0%–20%. Additionally, we set α+β=1 to enforce long-range asymptotic exact exchange. Then, the range-separation parameter γ is varied to achieve a minimization of the target function

J2(γ)=[IPγ(N)+EHOMOγ(N)]2+[IPγ(N+1)+EHOMOγ(N+1)]2,
(3)

where the ionization potential of the neutral species with N electrons, IPγ(N), is determined via a ΔSCF approach from total energy differences as IPγ(N)=ϵtotγ(N1)ϵtotγ(N). Here ϵtotγ(N) and ϵtotγ(N1) are the total energies of the neutral and cation species, respectively. This procedure enforces the ionization potential theorem of DFT,51–55 namely, that the energy of the Kohn-Sham highest occupied molecular orbital (HOMO) is equal to the negative of the first ionization potential. For molecules with an unbound N + 1 anionic state, only the first of these two terms is minimized, as in our previous work.45 The optimal parameters obtained within this framework for the molecules studied are listed in the supplementary material. OTRSH-BNL parameters are taken from Ref. 7.

Our TDDFT calculations are performed with QChem 4.256 with standard settings, excluding core electrons in the correlation computation and neglecting relativistic effects as usual. We use the cc-pVTZ basis set which, relative to aug-cc-pVTZ, converges the neutral-excitation energies satisfactorily: we consider TD-CAM-B3LYP, with and without the TDA, for the singlet La and Lb states for all acenes considered herein and the lowest lying triplet state for benzene, naphthalene, and anthracene. For all test cases, the difference between the augmented and unaugmented basis sets was between 0.001 and 0.087 eV, with an unsigned average difference of 0.028 eV.

In aromatic hydrocarbons, like the acenes, azulene, benzo[e]pyrene (BP), and dibenzo[a, c]- anthracene (DBAn) (see Figure 1), the two low-lying singlet excitations are labeled La1 and Lb1.57 These excitations are well-known to differ significantly in character, and these differences are widely discussed in the literature. The bright (or large oscillator strength) longitudinal ionic La1 state involves principally a transition between the highest-occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) and is often described as having charge-transfer (CT)-like character; while the dark (near-zero oscillator strength) covalent Lb1 excitation arises from a destructive interference58 of transitions that typically couple the HOMO to the LUMO+1 and the HOMO−1 to the LUMO.7,59,60

The description of these excitations as ionic or covalent comes from valence bond theory and refers to the distribution of charge in the spatial part of the component orbitals of the excited state. If in the resonance structures describing these orbitals, the density oscillates from negative to positive with respect to the carbon-atom centers, the excitation is termed “ionic.” If there is no such oscillation and the resonance structures correspond to Kekule structures, with alternating double and single bonds, the excitation is termed “covalent.”

The corresponding low-lying triplet excitations are labeled, following the same conventions as above, as La3 and Lb3, respectively. Notice that labeling of unbound molecular orbitals (e.g., the LUMO and LUMO+1 of some acenes) is somewhat arbitrary since their ordering may change depending on the choice of DFT exchange-correlation (XC) functional and basis set. In this work, we adopt a definition of the unbound LUMO of benzene and naphthalene as the first resonant state whose energy corresponds to the negative electron affinity energy measured in experiments, as detailed in Ref. 61. As this state is not stable, the resonant state has a finite lifetime. In this work, we first computed the LUMO within PBE. PBE spuriously binds the LUMO; however, the corresponding wavefunction will only serve as a basis to identify the LUMO within the more accurate approximations. We then define the LUMO calculated with a hybrid functional as the unbound state having the largest overlap with the PBE LUMO. With this simple method, we have been able to extract the LUMO states across the acene series in a consistent and reliable manner.

Predicting both La1 and Lb1 presents a challenge for TDDFT approaches. In fact, La1 excitations, with their CT-like character,7 are usually poorly predicted by standard TDDFT7,60,62 due to the known shortcoming of many standard functionals to describe such excitations. We note that this shortcoming has the potential to be ameliorated by using RSHs with asymptotic exact exchange,7 although the CT nature of La1 excitations and the ability of RSHs to overcome these shortcomings have been questioned.62 

We begin with a benchmark of GW-BSE and TDDFT against CCSD(T) for the low-lying singlet excitations of the acenes. We end with an examination of the role of the TDA within both the TDDFT and GW-BSE frameworks.

In Figure 2(a), we show the mean signed deviation (MSD=1/NiiNiEiEiref), relative to CCSD(T),8,61 of representative TDDFT-RSHs explored in this study. Our TD-OTRSH-BNL results are in excellent agreement with those reported in Refs. 7 and 8. In addition, we find that TD-OTRSH-PBE low-lying singlets are within 0.05 eV of the corresponding TD-OTRSH-BNL excitations. In fact, as previously discussed,8,62 the performance of TD-OTRSH based on BNL or PBE for La1 and Lb1 relative to CCSD(T) is very consistent: for both approaches, La1 is within 0.1–0.2 eV of the reference, but Lb1 presents larger discrepancies (0.4 eV), as does the Lb11La gap, which is within 0.6 eV.

FIG. 2.

MSD (see text for details) with respect to CCSD(T)8,61 of calculated neutral excitations of the acene molecules (n=1–6). The calculated La1 (blue bars), Lb1 (orange bars), and La3 (pink bars) excitations, and Lb11La (black bars) energy difference are shown for a few representative TDDFT and GW-BSE approaches: TD-OTRSH and TD-CAM-B3LYP in panel (a), and G0W0-BSE@BHLYP, G0W0-BSE@OTRSH-PBE, and evGW-BSE@PBE0 in panel (b).

FIG. 2.

MSD (see text for details) with respect to CCSD(T)8,61 of calculated neutral excitations of the acene molecules (n=1–6). The calculated La1 (blue bars), Lb1 (orange bars), and La3 (pink bars) excitations, and Lb11La (black bars) energy difference are shown for a few representative TDDFT and GW-BSE approaches: TD-OTRSH and TD-CAM-B3LYP in panel (a), and G0W0-BSE@BHLYP, G0W0-BSE@OTRSH-PBE, and evGW-BSE@PBE0 in panel (b).

Close modal

Having reviewed the accuracy of TDDFT-RSH for the La1 and Lb1 excitations relative to CCSD(T), we now discuss the GW-BSE results, focusing on the sensitivity to the underlying GW starting point. In Figure 2(b), we show the calculated MSD, as defined in Sec. IV A, of representative GW-BSE approaches studied here. Consistent with previously reported GW results on the charged excitations of the acenes (see, for instance, Refs. 42, 43, and 45), hybrid starting points for G0W0 or self-consistent GW approaches are required to predict accurate excitations within GW-BSE; relatively low MSDs are found within G0W0-BSE@BHLYP and evGW-BSE, in agreement with recent works.28–30 As hypothesized in Ref. 19, the OTRSH starting point is superior. In particular, we highlight that while the OTRSH-PBE starting point yields neutral excitation energies with accuracies similar to other starting-points for aromatic hydrocarbons, OTRSH-PBE leads to markedly improved triplet energetics for the molecules studied here, as discussed later. As shown in Figure 2(b), for the GW-BSE schemes considered here, the Lb1 state is predicted within 0.1–0.2 eV, whereas La1 is underestimated by at least 0.4 eV; additionally, the Lb1La1 gap is underestimated by 0.6 eV independent of GW-BSE scheme. The rather poor performance of both GW-BSE and TDDFT approaches, in the context of neutral low-lying singlet and triplet excitations of acene molecules, can be remedied by the Tamm-Dancoff approximation, as discussed next.

The fact that the TDA can improve the description of low-lying neutral excitations of the acenes has been discussed thoroughly in the TDDFT literature,4,5,9,35,63 and here we find that similar arguments apply to GW-BSE. In Figure 2, we also show the MSD of the calculated low-lying excited states using TDDFT and GW-BSE within the TDA with respect to the CCSD(T) reference. The calculated La1 and Lb1 energies of acene molecules within representative GW-BSE and TDDFT schemes are shown in Figure 3.

FIG. 3.

Low-lying singlet excitations of acenes calculated with TD-OTRSH-PBE and G0W0-BSE@OTRSH-PBE in panels (a) and (b); La1 and Lb1 excitation energies, with blue and orange lines, respectively, are compared to CCSD(T) references from Refs. 8 and 37 (dashed lines). The corresponding excitations with the TDA at the TDDFT and GW-BSE theories are shown in panels (c) and (d).

FIG. 3.

Low-lying singlet excitations of acenes calculated with TD-OTRSH-PBE and G0W0-BSE@OTRSH-PBE in panels (a) and (b); La1 and Lb1 excitation energies, with blue and orange lines, respectively, are compared to CCSD(T) references from Refs. 8 and 37 (dashed lines). The corresponding excitations with the TDA at the TDDFT and GW-BSE theories are shown in panels (c) and (d).

Close modal

1. Link between TDDFT and triplet instability

Within Hartree-Fock64 and within DFT,65 the stability of the spin-restricted solution against that of the more flexible spin-unrestricted solution requires the positive-definiteness of two matrices (one for singlet final states and one for triplet final states) that are precisely the sum of the blocks A and B [see Eq. (2)] used in time-dependent Hartree-Fock (TDHF) or in TDDFT. In other words, if either one of the A + B matrices has a negative eigenvalue, then the ground-state singlet solution is unstable against a spin-unrestricted triplet solution. This is the so-called triplet instability. Consequently, an unstable or near-unstable spin-restricted ground state implies negative or very small eigenvalues of A + B, which in turn produce non-physical or too-small neutral excitations in TDHF or in TDDFT.9,36 This is why prior work often resorts to the TDA to circumvent the spin-restricted instability situation.4,5,9,10,36,63 The TDA is thus a pragmatic way to prevent the electronic system from sampling the triplet ground state, which is spuriously too low in energy.

2. Link between the BSE and triplet instability

In the GW-BSE framework, the connection between triplet instability and the BSE matrix (Eq. (2)) is precisely analogous. However, the connection cannot be demonstrated as rigorously as for TDDFT. The BSE, evaluated in the standard fashion,24 is indeed a combination of (1) eigenvalues obtained from the dynamical GW self-energy and (2) a kernel, which is an approximate functional derivative of the static GW self-energy, namely, the static Coulomb hole plus screened exchange (COHSEX) approximation.66 Additionally, the functional derivative δW/δG is always neglected in the BSE kernel.24 Thus, following the same logic for GW-BSE as for TDDFT above, the BSE blocks A + B would then lead to stability problems (if present) in the static COHSEX spin-restricted solution. If one admits that the GW quasiparticle energies are not far from the static screened exchange energies, the connection between triplet instability and BSE can be understood, but again, not proven. In situations where the triplet instability occurs or nearly occurs in static COHSEX, the TDA to the BSE may be a good route to obtain meaningful neutral excitation energies. However, this calls for a direct numerical comparison, which we carry out below.

3. Performance of the TDA within TDDFT

As demonstrated in prior work,5,63 the TDA improves the description of the La1 singlet and the first triplet La3 states because these share a similar origin; both are covalent in the valence bond sense, and involve mainly HOMO to LUMO transitions, whereas Lb energies are virtually unmodified with the TDA. Hence, within the RSH time-dependent approaches used here, the large discrepancy (0.4 eV) between the calculated Lb1 state and CCSD(T) is not improved by the TDA (see Figure 2). On the other hand, we find that the TDA leads to an improvement in the La1 excitation energies in the asymptotic limit of longer acene molecules (see panels (a) and (c) of Figure 3). For example, for pentacene, TD-OTRSH-PBE predicts La1=2.18 eV with TDDFT and 2.48 eV with the TDA, in outstanding agreement with the CCSD(T) reference value of 2.42 eV. In brief, and in agreement with the literature,8,62 TDDFT-TDA with RSH functionals yields highly accurate CT-like La1 energetics but tends to overestimate Lb1 transition energies.

4. Performance of the TDA within GW-BSE

Having reviewed the ability of the TDA within TDDFT to predict the low-lying excitations of the acenes, we now discuss the accuracy of the TDA within GW-BSE for these transitions. Here we expand our discussion to a larger set of aromatic hydrocarbons, including azulene, BP, and DBAn (see Figure 1), which have well-characterized La and Lb states.62 In Table I, we show the calculated singlet and triplet excitations with G0W0-BSE@BHLYP (with and without the TDA) and mean deviations with respect to CCSD(T), as previously defined.

TABLE I.

Singlet and triplet energetics of representative aromatic hydrocarbons calculated with GW-BSE@BHLYP with the full-BSE (denoted simply BSE above) and the TDA. We consider benzene (Benz), naphthalene (Naph.), anthracene (Anth.), tetracene (Tetra.), pentacene (Penta.), hexacene (Hexa.), azulene (Azu.), benzo[e]pyrene (BP), and dibenzo[a, c]anthracene (DBAn). MSD and MAD with respect to CCSD(T) are also shown (see text). All energies are in units of eV.

𝑳𝒂𝐬𝐭𝐚𝐭𝐞𝐬
SingletsTriplets
BSETDACCSD(T)aBSETDACCSD(T)a
Benz. 5.89 6.13 6.5437  3.60 3.94 4.26±0.1137,61 
Naph. 4.34 4.59 4.81±0.028,37 2.65 2.93 3.20±0.1137,61 
Anth. 3.38 3.63 3.68±0.028,62 2.01 2.29 2.41±0.0761,62 
Tetra. 2.42 2.72 2.948  1.08 1.36 1.7661  
Penta. 1.88 2.21 2.428  0.57 0.91 1.3761  
Hexa. 1.48 1.84 2.058  <0b 0.58 1.0061  
Azu. 2.00 2.12 1.9462  1.35 1.42 2.1862  
BP 3.65 3.82 4.0962  1.95 2.41 2.8262  
DBAn 3.48 3.69 3.9162  1.90 2.30 2.7362  
MSD −0.43 −0.18  −0.74 −0.39  
MAD 0.44 0.22  0.74 0.39  
𝑳𝒂𝐬𝐭𝐚𝐭𝐞𝐬
SingletsTriplets
BSETDACCSD(T)aBSETDACCSD(T)a
Benz. 5.89 6.13 6.5437  3.60 3.94 4.26±0.1137,61 
Naph. 4.34 4.59 4.81±0.028,37 2.65 2.93 3.20±0.1137,61 
Anth. 3.38 3.63 3.68±0.028,62 2.01 2.29 2.41±0.0761,62 
Tetra. 2.42 2.72 2.948  1.08 1.36 1.7661  
Penta. 1.88 2.21 2.428  0.57 0.91 1.3761  
Hexa. 1.48 1.84 2.058  <0b 0.58 1.0061  
Azu. 2.00 2.12 1.9462  1.35 1.42 2.1862  
BP 3.65 3.82 4.0962  1.95 2.41 2.8262  
DBAn 3.48 3.69 3.9162  1.90 2.30 2.7362  
MSD −0.43 −0.18  −0.74 −0.39  
MAD 0.44 0.22  0.74 0.39  
𝑳𝒃𝐬𝐭𝐚𝐭𝐞𝐬
SingletsTriplets
BSETDACCSD(T)aBSETDACCSD(T)a
Benz. 5.10 5.15 5.0837  4.39 4.42 4.8637  
Naph. 4.30 4.32 4.19±0.068,37 3.76 3.79 4.0937  
Anth. 3.82 3.79 3.58±0.018,62 3.42 3.43 3.5262  
Tetra. 3.33 3.37 3.258  3.11 3.18  
Penta. 3.10 3.13 3.028  2.79 2.83  
Hexa. 2.91 2.98 2.868  b 2.66  
Azu. 3.34 3.49 3.6462  2.09 2.19 2.2062  
BP 3.57 3.60 3.5062  2.96 3.11 3.3462  
DBAn 3.58 3.61 3.5762  3.34 3.43 3.3562  
MSD 0.04 0.08  −0.23 −0.17  
MAD 0.11 0.12  0.23 0.20  
𝑳𝒃𝐬𝐭𝐚𝐭𝐞𝐬
SingletsTriplets
BSETDACCSD(T)aBSETDACCSD(T)a
Benz. 5.10 5.15 5.0837  4.39 4.42 4.8637  
Naph. 4.30 4.32 4.19±0.068,37 3.76 3.79 4.0937  
Anth. 3.82 3.79 3.58±0.018,62 3.42 3.43 3.5262  
Tetra. 3.33 3.37 3.258  3.11 3.18  
Penta. 3.10 3.13 3.028  2.79 2.83  
Hexa. 2.91 2.98 2.868  b 2.66  
Azu. 3.34 3.49 3.6462  2.09 2.19 2.2062  
BP 3.57 3.60 3.5062  2.96 3.11 3.3462  
DBAn 3.58 3.61 3.5762  3.34 3.43 3.3562  
MSD 0.04 0.08  −0.23 −0.17  
MAD 0.11 0.12  0.23 0.20  
a

CCSD(T) data from the literature.

b

The BSE Hamiltonian contains negative eigenvalues (read text).

Similar to TDDFT, in GW-BSE we find that, independent of GW self-energy scheme, the La1 and La3 states are improved within the TDA (by at least 0.2 eV, see Figure 2 and Table I). While triplet energies remain underestimated by 0.30.4 eV, singlet energies are accurately predicted with GW-BSE-TDA, with remaining discrepancies lower than 0.2 eV. The TDA also corrects the La1Lb1 energy ordering; as shown in Figure 3, these two states cross at naphthalene when following increasing/decreasing ring number n within the full-BSE (panel (b)), while the crossing is at anthracene within the TDA (panel (d)), in agreement to CCSD(T). In summary, GW-BSE within the TDA can predict—with excellent quantitative accuracy, an MSD better than 0.2 eV—both ionic CT-like and covalent singlet excitations (such as La1 and Lb1, respectively) of the acenes and other aromatic-hydrocarbons.

5. Thiel’s set

We further analyze the accuracy of GW-BSE and the TDA for a larger set of triplet excitations. We show in Figure 4 the MSDs (previously defined) of the calculated first-triplet energies of 20 organic molecules of Thiel’s set (all energies are tabulated in the supplementary material). Again, we consider several representative GW-BSE approaches using the BHLYP and OTRSH-PBE starting points for G0W0 and evGW. Additionally, we show the MSDs for the molecule categories described in Figure 1.

FIG. 4.

First-triplet excitation energies of organic molecules in Thiel’s set (see Figure 1) calculated with GW-BSE are benchmarked against reference data.37 The MSD (read text) corresponding to molecules in Series 1 is shown in blue bars, Series 2 in orange bars, Series 3 in black bars, and the total in pink bars. We consider several GW-BSE schemes with the full-BSE and the TDA.

FIG. 4.

First-triplet excitation energies of organic molecules in Thiel’s set (see Figure 1) calculated with GW-BSE are benchmarked against reference data.37 The MSD (read text) corresponding to molecules in Series 1 is shown in blue bars, Series 2 in orange bars, Series 3 in black bars, and the total in pink bars. We consider several GW-BSE schemes with the full-BSE and the TDA.

Close modal

With this larger set of excitations, it becomes clear from our calculations that G0W0@BHLYP and evGW@PBE0, known to perform reasonably well for singlet excitations,28,29 can present severe errors for triplets, with MSDs of 0.6 to −0.8 eV, as noticed first in Ref. 28 for the BHLYP starting point. The OTRSH starting point for G0W0-BSE has a relatively lower MSD of 0.4 to −0.6 eV, presumably due to the RSH optimal starting point for the underlying GW electronic structure.43,45 For all GW-BSE approaches studied here, the TDA improves the first-triplet energy, a fact that we discuss in greater depth below. Further, we note that in agreement with recent work,31GW-BSE-TDA approaches predict inaccurate triplet energies (with MSD of −0.4–0.5 eV) when using a global-hybrid starting point. Importantly, within the OTRSH starting-point, GW-BSE-TDA can result in relatively accurate first-triplet energies with a MSD of −0.19 eV.

In order to better understand the superior performance of the TDA within GW-BSE for low-lying triplet energies, we show in Figure 5 the ratio of the triplet energy calculated within the full GW-BSE and within the TDA (T/TTDA). When the ratio approaches zero or becomes negative (TDDFT predicts a negative or zero triplet energy), the triplet and its corresponding ground state become unstable, as explained before; hence this ratio acts as a measure of instability.9 In this work, we find that the full-BSE and the TDA predict similar triplet energies for benzene and naphthalene (ratio close to 0.9), but the triplet ratio drops to less than zero at hexacene independent of the GW approximation; note that for the PBE starting point to GW-BSE, the ratio becomes negative at pentacene (not shown). This implies that GW-BSE, in disagreement with CCSD(T),61 predicts triplet ground states for acenes larger than pentacene. In analogy to TDDFT, this may be a result of instabilities in the corresponding GW-BSE triplet and ground states; we leave the evaluation of stability conditions in the GW-BSE states to future work.

FIG. 5.

Ratio of the first triplet energy (T) calculated within GW-BSE diagonalizing the full BSE Hamiltonian and using the TDA. Several representative GW-BSE schemes are shown: G0W0-BSE@BHLYP in blue line with crosses, G0W0-BSE@OTRSH-PBE in orange line with circles, and evGW-BSE@PBE0 in black line with crosses. GW-BSE predicts a negative triplet energy (shown at zero) for hexacene for all GW schemes used in this work.

FIG. 5.

Ratio of the first triplet energy (T) calculated within GW-BSE diagonalizing the full BSE Hamiltonian and using the TDA. Several representative GW-BSE schemes are shown: G0W0-BSE@BHLYP in blue line with crosses, G0W0-BSE@OTRSH-PBE in orange line with circles, and evGW-BSE@PBE0 in black line with crosses. GW-BSE predicts a negative triplet energy (shown at zero) for hexacene for all GW schemes used in this work.

Close modal

As mentioned above, triplet instabilities are well-known and documented in Hartree Fock and TDDFT theories,4,5,9,10,35 and GW-BSE is similarly affected, as shown here and in Ref. 67. In TDDFT triplet instabilities are overcome with the TDA. In fact, TD-TDA triplet energetics result equivalent to those from configuration interaction singles (CIS) theory, which mixes only single Slater determinants and is the minimal post Hartree Fock method capable of predicting physical excited states.10 Not surprisingly, in GW-BSE, the TDA also overcomes triplet instabilities, as we document here, in a manner analogous to TDDFT for molecules.

Finally, we briefly comment on the performance of the BSE among other two-particle Green’s function methods. In this context, the algebraic diagrammatic construction [ADC(2)]68 and the second-order polarization propagator (SOPPA)69 methods are efficient Green’s function methods which give access to the neutral excitations of molecules. The accuracy of representative variants of ADC(2) has been studied for the low-lying singlet excitations of the acenes, from naphthalene to hexacene, in Ref. 70. Comparing these results to the CCSD(T) reference used in this work, we find an MSD for the ADC(2) variants of −0.46 and 0.36 eV (or lower) for La1 and Lb1, respectively, relatively higher than that of the best GW-BSE method used in this work (GW-BSE-TDA@OTRSH with an MSD of −0.13 and 0.17 eV for the same excitations).70–72 In Refs. 73 and 74, the lowest singlet excitations of small acenes have been calculated with SOPPA-based methods, among which the original SOPPA and SOPPA(CCSD) perform best. The MSD with respect to CCSD(T) is 0.3–0.8 eV for La1 and Lb1 and 0.4–0.5 eV for the first triplet energies of benzene and naphthalene.75 Moreover, for the triplet excitations in Thiel’s set, SOPPA and SOPPA(CCSD) yield MSDs with respect to the best theoretical estimate of −0.45 and −0.54 eV, respectively,74 which again is relatively high compared to G0W0-BSE-TDA@OTRSH in this work (with an MSD of only −0.2 eV). The superior performance of GW-BSE-TDA for these sets of excitations therefore situates the BSE as an efficient and accurate alternative to many traditional approximate methods in quantum chemistry.

In summary, we have benchmarked GW-BSE with CCSD(T) for neutral excitations of aromatic hydrocarbons and heterocycles, including the challenging La and Lb excitations for aromatic hydrocarbons heavily documented in prior work with TDDFT. We first explored the accuracy of approximations to GW-BSE and found that G0W0-BSE@OTRSH can yield accurate triplet and singlet excitations, sometimes outperforming other highly accurate approaches such as evGW-BSE and G0W0-BSE@BHLYP. In particular, for aromatic hydrocarbons, the above mentioned GW-BSE methods can predict accurate Lb1 energetics but generally present significant errors for the La1 states. This problem is remedied by using the TDA, which leads, as it does with TDDFT, to a better overall performance, overcoming triplet instabilities, improving triplet energetics and capturing quantitatively both the charge-transfer-like La and covalent Lb singlet excitations of aromatic cyclic compounds.

See supplementary material for the tabulation and calculation of neutral excitations of the acenes and other organic molecules of Thiel’s set.

This work was supported by the Center for Computational Study of Excited State Phenomena in Energy Materials at the Lawrence Berkeley National Laboratory, which is funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division under Contract No. DE-AC02-05CH11231, as part of the Computational Materials Sciences Program. This work is also supported by the Molecular Foundry through the U.S. Department of Energy, Office of Basic Energy Sciences under the same contract number. We acknowledge the use of computational resources at the National Energy Research Scientific Computing Center (NERSC). F. Bruneval acknowledges the Enhanced Eurotalent program and the France Berkeley Fund for supporting his sabbatical leave at UC Berkeley.

1.
J. E.
Anthony
,
Angew. Chem., Int. Ed.
47
,
452
(
2008
).
2.
M. B.
Smith
and
J.
Michl
,
Annu. Rev. Phys. Chem.
64
,
361
(
2013
).
3.
J.
Lee
,
P.
Jadhav
,
P. D.
Reusswig
,
S. R.
Yost
,
N. J.
Thompson
,
D. N.
Congreve
,
E.
Hontz
,
T.
Van Voorhis
, and
M. A.
Baldo
,
Acc. Chem. Res.
46
,
1300
(
2013
).
4.
M. J. G.
Peach
,
M. J.
Williamson
, and
D. J.
Tozer
,
J. Chem. Theory Comput.
7
,
3578
(
2011
).
5.
M. J. G.
Peach
and
D. J.
Tozer
,
J. Phys. Chem. A
116
,
9783
(
2012
).
6.
S.
Grimme
and
M.
Parac
,
ChemPhysChem
4
,
292
(
2003
).
7.
N.
Kuritz
,
T.
Stein
,
R.
Baer
, and
L.
Kronik
,
J. Chem. Theory Comput.
7
,
2408
(
2011
).
8.
K.
Lopata
,
R.
Reslan
,
M.
Kowalska
,
D.
Neuhauser
,
N.
Govind
, and
K.
Kowalski
,
J. Chem. Theory Comput.
7
,
3686
(
2011
).
9.
J. S.
Sears
,
T.
Koerzdoerfer
,
C.-R.
Zhang
, and
J.-L.
Brédas
,
J. Chem. Phys.
135
,
151103
(
2011
).
10.
M. E.
Casida
and
M.
Huix-Rotllant
,
Annu. Rev. Phys. Chem.
63
,
287
(
2012
).
11.
B. M.
Wong
and
T. H.
Hsieh
,
J. Chem. Theory Comput.
6
,
3704
(
2010
).
12.
R. M.
Richard
and
J. M.
Herbert
,
J. Chem. Theory Comput.
7
,
1296
(
2011
).
13.
T.
Leininger
,
H.
Stoll
,
H.-J.
Werner
, and
A.
Savin
,
Chem. Phys. Lett.
275
,
151
(
1997
).
14.
T.
Yanai
,
D. P.
Tew
, and
N. C.
Handy
,
Chem. Phys. Lett.
393
,
51
(
2004
).
15.
L.
Kronik
,
T.
Stein
,
S.
Refaely-Abramson
, and
R.
Baer
,
J. Chem. Theory Comput.
8
,
1515
(
2012
).
16.
G.
Borghi
,
A.
Ferretti
,
N. L.
Nguyen
,
I.
Dabo
, and
N.
Marzari
,
Phys. Rev. B
90
,
075135
(
2014
).
17.
R.
Baer
and
D.
Neuhauser
,
Phys. Rev. Lett.
94
,
043002
(
2005
).
18.
E.
Livshits
and
R.
Baer
,
Phys. Chem. Chem. Phys.
9
,
2932
(
2007
).
19.
S.
Refaely-Abramson
,
S.
Sharifzadeh
,
N.
Govind
,
J.
Autschbach
,
J. B.
Neaton
,
R.
Baer
, and
L.
Kronik
,
Phys. Rev. Lett.
109
,
226405
(
2012
).
20.
S.
Refaely-Abramson
,
S.
Sharifzadeh
,
M.
Jain
,
R.
Baer
,
J. B.
Neaton
, and
L.
Kronik
,
Phys. Rev. B
88
,
081204
(
2013
).
21.
T.
Stein
,
L.
Kronik
, and
R.
Baer
,
J. Chem. Phys.
131
,
244119
(
2009
).
22.
23.
T.
Stein
,
L.
Kronik
, and
R.
Baer
,
J. Am. Chem. Soc.
131
,
2818
(
2009
).
24.
G.
Onida
,
L.
Reining
, and
A.
Rubio
,
Rev. Mod. Phys.
74
,
601
(
2002
).
25.
G.
Onida
,
L.
Reining
,
R. W.
Godby
,
R.
Del Sole
, and
W.
Andreoni
,
Phys. Rev. Lett.
75
,
818
(
1995
).
26.
J. C.
Grossman
,
M.
Rohlfing
,
L.
Mitas
,
S. G.
Louie
, and
M. L.
Cohen
,
Phys. Rev. Lett.
86
,
472
(
2001
).
27.
P.
Boulanger
,
D.
Jacquemin
,
I.
Duchemin
, and
X.
Blase
,
J. Chem. Theory Comput.
10
,
1212
(
2014
).
28.
F.
Bruneval
,
S. M.
Hamed
, and
J. B.
Neaton
,
J. Chem. Phys.
142
,
244101
(
2015
).
29.
D.
Jacquemin
,
I.
Duchemin
, and
X.
Blase
,
J. Chem. Theory Comput.
11
,
3290
(
2015
).
30.
D.
Jacquemin
,
I.
Duchemin
, and
X.
Blase
,
J. Chem. Theory Comput.
11
,
5340
(
2015
).
31.
D.
Jacquemin
,
I.
Duchemin
,
A.
Blondel
, and
X.
Blase
,
J. Chem. Theory Comput.
13
,
767
(
2017
).
32.
M.
Rohlfing
and
S. G.
Louie
,
Phys. Rev. Lett.
81
,
2312
(
1998
).
33.
G.
Strinati
,
Riv. Nuovo Cimento
11
,
1
(
2008
).
34.
A. L.
Fetter
and
J. D. Q.
Walecka
,
Quantum Theory of Many-Particle Systems
(
Dover Publications
,
Mineola, NY
,
2003
).
35.
S.
Hirata
and
M.
Head-Gordon
,
Chem. Phys. Lett.
314
,
291
(
1999
).
36.
M. E.
Casida
,
F.
Gutierrez
,
F.-X.
Guan
,
J.
Gadea
,
D.
Salahub
, and
J.-P.
Daudey
,
J. Chem. Phys.
113
,
7062
(
2000
).
37.
M. R.
Silva-Junior
,
M.
Schreiber
,
S. P. A.
Sauer
, and
W.
Thiel
,
J. Chem. Phys.
133
,
174318
(
2010
).
38.
M. R.
Silva-Junior
,
S. P. A.
Sauer
,
M.
Schreiber
, and
W.
Thiel
,
Mol. Phys.
108
,
453
(
2010
).
39.
F.
Bruneval
,
T.
Rangel
,
S. M.
Hamed
,
M.
Shao
,
C.
Yang
, and
J. B.
Neaton
,
Comput. Phys. Commun.
208
,
149
(
2016
).
40.
F.
Bruneval
, “
Molgw: A slow but accurate many-body perturbation theory code
,” https://github.com/bruneval/molgw (
2015
).
41.
X.
Blase
,
C.
Attaccalite
, and
V.
Olevano
,
Phys. Rev. B
83
,
115103
(
2011
).
42.
F.
Bruneval
and
M. A. L.
Marques
,
J. Chem. Theory Comput.
9
,
324
(
2013
).
43.
L.
Gallandi
,
N.
Marom
,
P.
Rinke
, and
T.
Körzdörfer
,
J. Chem. Theory Comput.
12
,
605
(
2016
).
44.
J. W.
Knight
,
X.
Wang
,
L.
Gallandi
,
O.
Dolgounitcheva
,
X.
Ren
,
J. V.
Ortiz
,
P.
Rinke
,
T.
Körzdörfer
, and
N.
Marom
,
J. Chem. Theory Comput.
12
,
615
(
2016
).
45.
T.
Rangel
,
S. M.
Hamed
,
F.
Bruneval
, and
J. B.
Neaton
,
J. Chem. Theory Comput.
12
,
2834
(
2016
).
46.
J. B.
Foresman
and
A.
Frisch
,
Exploring Chemistry with Electronic Structure Methods: A Guide to Using Gaussian
, 2nd ed. (
Gaussian
,
Pittsburgh, PA
,
1996
).
47.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
48.
T. H.
Dunning
,
J. Chem. Phys.
90
,
1007
(
1989
).
49.
O.
Vahtras
,
J.
Almlöf
, and
M.
Feyereisen
,
Chem. Phys. Lett.
213
,
514
(
1993
).
50.
F.
Weigend
,
Phys. Chem. Chem. Phys.
4
,
4285
(
2002
).
51.
J. P.
Perdew
,
R. G.
Parr
,
M.
Levy
, and
J. L.
Balduz
,
Phys. Rev. Lett.
49
,
1691
(
1982
).
52.
U.
Salzner
and
R.
Baer
,
J. Chem. Phys.
131
,
231101
(
2009
).
53.
M.
Levy
,
J. P.
Perdew
, and
V.
Sahni
,
Phys. Rev. A
30
,
2745
(
1984
).
54.
J. P.
Perdew
and
M.
Levy
,
Phys. Rev. B
56
,
16021
(
1997
).
55.
C.-O.
Almbladh
and
U.
von Barth
,
Phys. Rev. B
31
,
3231
(
1985
).
56.
Y.
Shao
 et al,
Phys. Chem. Chem. Phys.
8
,
3172
(
2006
).
57.
J. R.
Platt
,
J. Chem. Phys.
17
,
484
(
1949
).
58.
E. B.
Guidez
and
C. M.
Aikens
,
J. Phys. Chem. C
117
,
21466
(
2013
).
59.
K.
Hirao
,
H.
Nakano
, and
K.
Nakayama
,
J. Chem. Phys.
107
,
9966
(
1997
).
60.
M.
Parac
and
S.
Grimme
,
Chem. Phys.
292
,
11
(
2003
).
61.
B.
Hajgató
,
D.
Szieberth
,
P.
Geerlings
,
F.
De Proft
, and
M. S.
Deleuze
,
J. Chem. Phys.
131
,
224321
(
2009
).
62.
B.
Moore
,
H.
Sun
,
N.
Govind
,
K.
Kowalski
, and
J.
Autschbach
,
J. Chem. Theory Comput.
11
,
3305
(
2015
).
63.
Y.-L.
Wang
and
G.-S.
Wu
,
Int. J. Quantum Chem.
108
,
430
(
2008
).
64.
R.
Seeger
and
J. A.
Pople
,
J. Chem. Phys.
66
,
3045
(
1977
).
65.
R.
Bauernschmitt
and
R.
Ahlrichs
,
J. Chem. Phys.
104
,
9047
(
1996
).
67.
R.
Zimmermann
,
Phys. Status Solidi
41
,
23
(
1970
).
68.
A. B.
Trofimov
and
J.
Schirmer
,
J. Phys. B
28
,
2299
(
1995
).
69.
E. S.
Nielsen
,
P.
Jørgensen
, and
J.
Oddershede
,
J. Chem. Phys.
73
,
6238
(
1980
).
70.
S.
Knippenberg
,
J. H.
Starcke
,
M.
Wormit
, and
A.
Dreuw
,
Mol. Phys.
108
,
2801
(
2010
).
71.

The calculated singlet energies of the acenes with strict ADC(2) [ADC(s)-s] in Ref. 70 are 4.57, 4.00, 3.61, 3.36, and 3.18 eV for Lb1, from naphthalene to hexacene, respectively, and 5.09, 3.87, 3.04, 2.46, and 2.05 eV for La1. The extended variant [ADC(2)-x] yields poorer results for these excitations.70 Similar values for the singlet excitations of naphthalene are found in Ref. 72 with a smaller TZVP basis. We use these results to calculate the MSD of ADC(2) with respect to the CCSD(T) references (shown in the supplementary material).

72.
B.
Helmich
and
C.
Hättig
,
Comput. Theor. Chem.
1040–1041
,
35
(
2014
).
73.
M. J.
Packer
,
E. K.
Dalskov
,
T.
Enevoldsen
,
H. J. A.
Jensen
, and
J.
Oddershede
,
J. Chem. Phys.
105
,
5886
(
1996
).
74.
S. P. A.
Sauer
,
H. F.
Pitzner-Frydendahl
,
M.
Buse
,
H. J. A.
Jensen
, and
W.
Thiel
,
Mol. Phys.
113
,
2026
(
2015
).
75.

We have compiled the calculated energies with SOPPA-based methods from Ref. 74 with the aug-cc-TZVP basis, which are also found in Ref. 73 with a smaller basis set. The La1 energies of benzene and naphthalene calculated with SOPPA are 5.91 and 4.19, respectively, and 5.77 and 3.97 with SOPPA(CCSD). The corresponding Lb1 energies are 4.63 and 3.78 eV with SOPPA and 4.43 and 3.49 with SOPPA(CCSD). The first triplet energies within SOPPA and the cc-pVTZ basis are 3.73 and 2.68 eV for benzene and naphthalene, respectively; the corresponding values within SOPPA(CCSD) are 3.56 and 2.44 eV.

Supplementary Material