Is the Bethe-Salpeter Formalism Accurate for Excitation Energies? Comparisons with TD-DFT, CASPT2, and EOM-CCSD.

Developing ab initio approaches able to provide accurate excited-state energies at a reasonable computational cost is one of the biggest challenges in theoretical chemistry. In that framework, the Bethe-Salpeter equation approach, combined with the GW exchange-correlation self-energy, which maintains the same scaling with system size as TD-DFT, has recently been the focus of a rapidly increasing number of applications in molecular chemistry. Using a recently proposed set encompassing excitation energies of many kinds [J. Phys. Chem. Lett. 2016, 7, 586-591], we investigate here the performances of BSE/GW. We compare these results to CASPT2, EOM-CCSD, and TD-DFT data and show that BSE/GW provides an accuracy comparable to the two wave function methods. It is particularly remarkable that the BSE/GW is equally efficient for valence, Rydberg, and charge-transfer excitations. In contrast, it provides a poor description of triplet excited states, for which EOM-CCSD and CASPT2 clearly outperform BSE/GW. This contribution therefore supports the use of the Bethe-Salpeter approach for spin-conserving transitions.

The determination of excited-state (ES) energies and properties is certainly a field in which theoretical chemistry can play a major role because ESs are both hard to capture with experimental methods and directly related to many key applications, for example, solar energy conversion, photocatalysis, and optical information storage. Today, when one wishes to model optical spectra or other ES properties for real-life molecules, one generally turns toward time-dependent density functional theory (TD-DFT) and, more precisely, to its linear-response adiabatic formulation.[1] Although this method, which benefits from a formal scaling with system size, is very successful,[2,3] it suffers from several limitations. Among these drawbacks are the impossibility to model ES with a significant multiple excitation character and the inadequacy of TD-DFT for multiconfigurational cases. In addition, TD-DFT’s excitation energies tend to significantly depend on the selected exchange-correlation functional. In practice, one can often find a functional that is suited for a given ES; for example, one can select a range-separated hybrid when studying an electronic transition presenting a significant charge-transfer (CT) character, but it remains difficult to obtain a balanced description of ES of different characters (local, Rydberg, CT, etc.) within the same compound. To achieve such description, one has to go for more refined approaches, which are unfortunately more computationally demanding. Although there have been many important developments in the field recently,[4−6] multireference methods, such as the accurate CASPT2 or MRCI theories, often remain too expensive for tackling large conjugated molecules, and, to go beyond TD-DFT, one often turns toward electron-correlated single-reference wave function schemes. In this vein, the most effective approach is probably ADC(2),[7] which delivers results similar to the ones of the more computationally involved CC2 approach[7−9] but nevertheless implies a formal scaling with system size. Although both ADC(2) and CC2 have the obvious advantage of not requiring the selection of an exchange-correlation functional, the obtained accuracy does not necessarily outperform its TD-DFT counterpart for excitation energies.[10,11] To go further, one can use EOM-CCSD or ADC(3), but these approaches are limited by their scaling. Another alternative that has been enjoying a rapidly increasing interest in molecular chemistry lately[12−18] is the Bethe–Salpeter equation (BSE) approach,[19−21] which is a specific implementation of Green’s function many-body perturbation theory.[21] BSE has been notably used to model polyacenes.[22−24] The BSE formalism relies, in particular, on accurate occupied/virtual energy levels, calculated within the GW formalism,[25] and the explicit nonlocal electron–hole interaction mediated by the screened Coulomb potential (W(, )). This latter interaction explains, in particular, that charge-transfer excitations are properly described.[26,27] Furthermore, BSE has been shown to accurately reproduce excitation energies in cyanine derivatives, for which TD-DFT calculations fail independently of the use of global or range-separated hybrids,[28] and several transitions with multiple-excitation character were observed to be better described than with TDDFT.[29] Finally, while proceeding technically as a perturbative correction based on input DFT Kohn–Sham eigenstates, which depend on the chosen exchange-correlation functional, the BSE results become (almost) independent of the selected starting point when adopting a (partially) self-consistent GW method, while conserving TD-DFT’s scaling and hence TD-DFT’s computational efficiency.[17]

To date, there have been only a few “general” benchmarks devoted to BSE/GW[9,14,29−31] and aiming at assessing the pros and cons of this approach for molecular applications. Two strategies have been principally used. On the one hand, comparisons with vertical excitation energies determined with highly correlated reference data obtained at the CASPT2 and CC3 levels on small molecules[32,33] have been carried out.[29,31] These comparisons have highlighted that performing BSE on top of eigenstates obtained through a perturbative GW correction of the DFT results, that is, using BSE/G0W0, implies a significant dependency on the selected DFT functional. Besides, such comparisons are inherently limited by the size of the considered compounds. For example, charge-transfer excited states cannot be evaluated on tiny compounds. In contrast, our groups have performed comparisons of BSE/GW 0–0 energies with experimental references.[9] While this approach allows us to consider larger molecules, it is also limited; for example, only fluorescent molecules can be treated, and an approximate model accounting for environmental effects has to be used, so that decomposing the total deviations into its different contributions is not possible. In the present contribution, we provide a complementary light by using the recently proposed Truhlar–Gagliardi (TG) set of electronic transitions.[6] This set is well balanced as it encompasses local, Rydberg, and CT states, considers both spin-conserving and singlet–triplet transitions, and does not imply the selection of a specific environmental model. It was originally designed by these authors to assess the performances of their multiconfiguration pair-density functional theory (MC-PDFT), and they concluded that this approach was both balanced (similar deviations for all families of states) and as accurate as CASPT2.[6] Herein we use exactly the same set of compounds and show that the BSE/evGW method (in which a partially self-consistent approach is used for the GW eigenvalues) is also a very competitive scheme compared with both CASPT2 and EOM-CCSD, benefiting from the facts that its computational requirements are much smaller than these two approaches and that no active space needs to be defined (in contrast with all multiconfigurational schemes).

Our main results are listed in Table and have been obtained with the aug-cc-pVTZ basis set but for the EOM-CCSD data on the two largest compounds (see Computational Details). As the majority of the reference values used in the TG set[6] are originating from experiment (see Table S-1 in the SI), we have first optimized the ground-state structures at the CCSD level using large atomic basis sets so as to improve over previously published geometries. In the SI, we compare the TD-PBE0/aug-cc-pVTZ results computed on these geometries to the TD-PBE0 values of ref (6) calculated with the same method with various basis sets and show that they are similar.

Table 1

Excitation Energies Determined for Truhlar–Gagliardi’s Set of Compoundsa

molecule	state	CCSD	CASPT2	TD-M06	TD-M06-2X	BSE	ref
acetaldehyde	1A″ (n → π*)	4.40	4.27	4.26	4.08	4.26	4.28
acetone	1A2 (n → π*)	4.58	4.44	4.43	4.26	4.28	4.43
formaldehyde	1A2 (n → π*)	4.06	3.92	3.89	3.73	3.87	4.00
pyrazine	1B3u (n → π*)	4.32	4.04	3.86	3.96	4.09	3.97
pyridazine	1B1 (n → π*)	4.02	3.67	3.44	3.65	3.78	3.60
pyridine	1B1 (n → π*)	5.19	5.06	4.71	4.87	5.03	4.74
pyrimidine	1B1 (n → π*)	4.67	4.38	4.21	4.44	4.49	4.18
s-tetrazine	1B3u (n → π*)	2.62	2.56	2.05	2.22	2.29	2.25
ethylene	1B1u (π → π*)	8.05	8.16	7.11	7.53	7.44	8.02
butadiene	1Bu (π → π*)	6.41	6.51	5.48	5.85	5.87	6.21
benzene	1B2u (π → π*)	5.22	4.83	5.28	5.58	5.21	4.90
	3B1u (π → π*)	4.02	3.98	3.70	4.33	3.57	4.12
naphthalene	1B3u (π → π*)	4.48	4.21	4.33	4.70	4.41	4.00
	3B2u (π → π*)	3.12	3.18	2.72	3.22	2.69	3.11
furan	1B2 (π → π*)	6.57	6.32	5.78	6.23	6.01	6.06
	3B2 (π → π*)	4.22	4.08	3.82	4.30	3.50	4.17
hexatriene	1Bu (π → π*)	5.55	5.32	4.64	4.95	4.98	4.93
	3Bu (π → π*)	2.78	2.67	2.24	2.67	2.22	2.69
water	Ryd. Singlet 2px → 3s	7.63	7.61	6.54	7.48	7.44	7.40
	Ryd. Triplet 2px → 3s	7.24	7.19	6.29	7.15	6.80	7.00
pNA	CT 1A1 (π → π*)	4.69	4.62	4.05	4.51	4.36	4.30
DMABN	CT 1A1 (π → π*)	5.04	4.87	4.64	4.95	4.80	4.57
B-TCNE	CT 1A1 (π → π*)	3.99	3.84	2.21	2.96	3.42	3.59

CCSD stands for EOM-CCSD, CAM for CAM-B3LYP, and BSE for BSE/evGW@M06-2X. The CASPT2 values are taken from ref (6), and the interested reader will find details regarding the impact of the size of the active space and IEPA shift in this work. The reference values are the same as in ref (6) and are taken from various sources (see the SI for details). All values are in electronvolts. pNA, DMABN, and B-TCNE stand for para-nitroaniline, dimethylaminobenonitrile, and benzene–tetracyanoethylene complex, respectively.

In a second step, we have compared the TD-DFT results obtained with two functionals, namely, M06 and M06-2X, to the corresponding BSE/evGW values determined starting with the same two functionals (see Figure and the SI for the BSE/evGW@M06 values). At the TD-DFT level, going from M06 to M06-2X tends to significantly upshift the transition energies, as expected when increasing the exact exchange ratio in the functional. It is only for the n → π* transitions of acetaldehyde, acetone, and formaldehyde that the opposite effect is found. On average, the increase attains 0.34 eV and is particularly strong for the two Rydberg transitions as well as for the B-TCNE CT excited state, a trend consistent with literature.[34,35] In contrast, the impact of the functional on the TD-DFT n → π* transition energies is more limited, which also fits the conclusions of previous reports.[36] Let us now turn toward the BSE/evGW results. In this case, an average upshift is also observed but is much smaller, 0.09 eV on average. This confirms that the partial self-consistent evGW procedure washes out most of the functional dependency, and this improvement is particularly impressive for the π → π*, Rydberg, and CT states that are much less sensitive to the starting functional than with TD-DFT. Nevertheless, small variations related to the frozen eigenvectors pertain with BSE/evGW. From Figure , one notices that these variations remain sizable for the n → π* transitions; that is, they remain on the same order of magnitude as with TD-DFT. Comparing the BSE/evGW singlet–singlet and singlet–triplet transitions in a given compound, one also notices that the latter are more dependent on the selected functional but less dependent than with TD-DFT, a conclusion that is inline with our most recent study.[37]

Figure 1

Variation of the excitation energies when changing the functional from M06 to M06-2X at the TD-DFT (top) and BSE/evGW (bottom) levels. The states are ordered as in Table . The blue, red, and green histograms correspond to valence, Rydberg, and CT transitions, respectively. The stars indicate the singlet–triplet transitions.

We now turn toward comparisons with reference values for the TG set. Table lists the mean signed error (MSE), mean unsigned error (MUE), and linear determination coefficient obtained for selected methods. Histogram plots and complementary data can be found in the SI. Let us first examine the results of the two wave function approaches, namely, EOM-CCSD and CASPT2. Both deliver very large correlations with experiment and positive MSE for all subsets of excitations; that is, they tend to overestimate the reference transition energies but provide very consistent trends (see also Figures S4 and S5 in the SI). In the literature, there have been contrasted reports regarding the accuracy of EOM-CCSD for valence and Rydberg states,[10,38] and our results indicate that the errors can be quite significant for CT transitions, a result not highlighted before to the very best of our knowledge. As expected, CASPT2 delivers accurate data with MUE as small as 0.15 eV for valence excitations, valuable results for Rydberg states, and slightly less satisfying predictions for CT states. As can be seen in Table S-6, both EOM-CCSD and CAS-PT2 are accurate for both spin-conserving and spin-changing transitions. Of the two tested functionals, M06-2X is undoubtedly the most satisfying in the TD-DFT framework, with a small MUE for valence states but quite large deviations for the CT states, as expected for a global hybrid. Nevertheless, we note that the range-separated CAM-B3LYP does not significantly improve the CT description in the present case (see Table S-5). The correlation with the reference values provided by TD-M06-2X is smaller than with wave function methods but remains large (R2 > 0.95). Clearly TD-M06-2X is a cost-effective method but for CT transitions, a conclusion in agreement with previous TD-DFT benchmarks.[39−43] In contrast, TD-M06 strongly undershoots both Rydberg and CT states, while TD-PBE0 yields, as expected,[44] rather accurate Rydberg states (Table S-5). BSE/evGW@M06-2X improves very significantly the accuracy of CT states compared with TD-M06-2X and stands as the most accurate method in Table for this class of transitions. It also provides accurate Rydberg transition energies. For the valence states, the average deviations remain rather large, about twice the one of CAS-PT2, but similar to the EOM-CCSD and MC-PDFT deviations. By analyzing more closely the BSE results (see Figure S-2 and Table S-6), we notice that these larger errors can be mainly ascribed to, on the one hand, the singlet–triplet transitions for which BSE is clearly much less accurate than most other methods,[37] and, on the other hand, ethylene, a tiny molecule for which BSE is known to be inaccurate.[29] Overall, using TG’s formula (see footnote b in Table ) to estimate the total deviation, the Bethe–Salpeter approach appears to be the method delivering both the smallest MSE and the smallest MUE. While such conclusions should be ascertained in the future for a larger set of molecules, it is certainly remarkable that no subset of transitions yields BSE MUE exceeding 0.3 eV, a success that is shared by CASPT2 and MC-PDFT only, two multiconfigurational approaches.

Table 2

Mean Signed and Unsigned Errors Given in Electronvolts for Different Subsets of Excitationsa

set	error	CCSD	CASPT2	MC-PDFTb	TD-M06	TD-M06-2X	BSE/evGW@M06-2X
Valence	MSE	0.26	0.11	–0.03	–0.21	0.05	–0.09
	MUE	0.27	0.15	0.29	0.29	0.22	0.28
Rydberg	MSE	0.24	0.20	0.07	–0.78	0.11	–0.08
	MUE	0.24	0.20	0.07	0.78	0.11	0.12
CT	MSE	0.42	0.29	–0.21	–0.52	–0.01	0.04
	MUE	0.42	0.29	0.24	0.57	0.40	0.15
Totalc	MSE	0.30	0.20	–0.06	–0.50	0.05	–0.05
	MUE	0.31	0.21	0.20	0.55	0.25	0.18
R2		0.982	0.989	0.932	0.919	0.952	0.953

CASPT2 and MC-PDFT values are from ref (6). At the bottom of the Table, we provide the linear determination coefficient, R2, obtained by comparing the full set of data to the reference values.

tPBE functional.

Computed as 1/3 Valence +1/3 Ry +1/3 CT, following ref (6).

To provide additional insights into the success of the Bethe–Salpeter approach, we have performed aug-cc-pVTZ CCSD(T) calculations of the HOMO–LUMO gap, excluding the largest DMABN and B-TCNE systems (Table S-4 in the SI), showing that the gaps predicted by the evGW@M06 and evGW@M06-2X methods are within 0.13 and 0.21 eV (MSE) of the CCSD(T) reference, respectively. This 0.08 eV MSE difference between HOMO–LUMO gaps can be compared with the 0.10 eV MSE difference between the corresponding BSE/evGW@M06 and BSE/evGW@M06-2X calculations for singlet excitations, confirming the good average correlation between the GW gap and the corresponding BSE excitation energies.[31] Consistently with ref (31), we note, however, that the GW and BSE/GW MSEs are not identical, with a small difference of 0.2 to 0.3 eV, so that it cannot be strictly concluded that an exact BSE/GW excitation energy corresponds to an exact GW HOMO–LUMO gap. Furthermore, beyond average correlations, rather large variations exist, such as in the case of ethylene, for which the error on the BSE excitation energy (0.58 eV) is much larger than that on the GW gap (0.09 eV), calling for further exploration of the pros and cons of the BSE formalism applied to organic systems.

In short, using the representative TG set of transition energies, we have shown that BSE/evGW is an effective method for studying the excitation energies of molecules. Despite its single-reference nature, its small dependency on the selected functional, and its favorable scaling, this theoretical method provides a balanced description of all families of excited-states, a success only achieved previously with multireference approaches or computationally expensive high-order wave function schemes. The obtained average unsigned deviation is comparable to the one of CASPT2, but the correlation coefficient with the reference value is smaller than that with both CASPT2 and EOM-CCSD. This effect is mainly due to the presence of singlet–triplet transitions in the TG sets, a subset for which BSE/evGW results are poor. We can therefore advocate the use of BSE/evGW for investigating spin-conserving transitions in large molecular systems. Of course, BSE/evGW also has its limitations; for example, it is it not suited for significantly multiconfigurational cases as CASPT2, and it is not systematically improvable by increasing the perturbation order as CC.

As a further consideration, we emphasize that the present BSE/evGW formalism proceeds as a perturbation theory starting from the monodeterminental Kohn–Sham ansatz ground-state wave function. Furthermore, we performed BSE calculation exploiting the statically screened Coulomb potential; namely, we adopted the analog of the adiabatic approximation used in TD-DFT. As such, the consideration that BSE calculations with dynamical screening could help in describing transitions with multiple excitation character[45] cannot be invoked here. It is therefore difficult to conclude at the present stage that the improvement over TD-DFT calculations for the present set of molecules stems potentially from a better description of cases presenting a multideterminental character. We can rather conclude that the use of the calculated microscopic screened-Coulomb potential W(, ) offers a much more flexible description of the short-, medium-, and long-range electron–hole interactions as compared with global or range-separated hybrid functionals that rely on a more “rigid” functional form. While the present study advocates the Bethe–Salpeter formalism for its remarkable accuracy and favorable scaling, one additional conclusion is that the standard adiabatic TD-DFT approach may still be improved by pursuing the search for better functionals.

Computational Details

At the exception of the B-TCNE structure that was taken in the literature,[46] the ground-state geometries of all compounds were optimized at the CCSD level with a large atomic basis set, typically def2-TZVPP (see the SI for Cartesian coordinates). The optimizations were performed with the Gaussian09 program,[47] imposing a tight convergence threshold on the forces and using 10–8 and 10–10 a.u. convergence thresholds on the CCSD and SCF energies, respectively. The EOM-CCSD and TD-DFT calculations were performed with the same program, applying Dunning’s aug-cc-pVTZ atomic basis set, but for the coupled-cluster calculations on DMABN and B-TCNE that were modeled using the aug-cc-pVDZ atomic basis set for obvious computational reasons. The CCSD(T) calculations in the SI were performed with Gaussian09 as well.[47] The GW and BSE calculations used aug-cc-pVTZ and have been achieved with the Fiesta package,[48] applying the Coulomb-fitting resolution-of-identity (RI-V) technique. It is important to stress that we used here the so-called evGW level; that is, we self-consistently converged the occupied/virtual GW energy levels while keeping the input DFT eigenfunctions frozen. We corrected at the GW level all valence occupied DFT levels and twice that amount of virtual levels. The BSE calculations were made with the “full” matrix; that is, we do not apply the Tamm–Dancoff approximation here. For the GW calculations, the starting DFT eigenstates were generated with NWChem[49] using the M06 and M06-2X hybrid functionals[50] and applying the xfine integration grid as well as 10–7 and 10–6 a.u. convergence thresholds for the total energies and densities, respectively.