GW

Article pubs.acs.org/JCTC

Assessment of the Accuracy of the Bethe−Salpeter (BSE/GW) Oscillator Strengths Denis Jacquemin,*,†,‡ Ivan Duchemin,¶,§ Aymeric Blondel,† and Xavier Blase*,∥,§ †

CEISAM Laboratory−UMRS CNR 6230, University of Nantes, 2 Rue de la Houssinière, BP 92208, 44322 Nantes Cedex 3, France Institut Universitaire de France, 1 rue Descartes, 75005 Paris Cedex 5, France ¶ Institute for Nanoscience and Cryogenics (INAC), SP2M/L_Sim, CEA/UJF Cedex 09, 38054 Grenoble, France § Grenoble Alpes University and ∥National Center for Scientific Research (CNRS), NEEL Institute, F-38042 Grenoble, France ‡

ABSTRACT: Aiming to assess the accuracy of the oscillator strengths determined at the BSE/GW level, we performed benchmark calculations using three complementary sets of molecules. In the first, we considered ∼80 states in Thiel’s set of compounds and compared the BSE/GW oscillator strengths to recently determined ADC(3/2) and CC3 reference values. The second set includes the oscillator strengths of the low-lying states of 80 medium to large dyes for which we have determined CC2/aug-cc-pVTZ values. The third set contains 30 anthraquinones for which experimental oscillator strengths are available. We find that BSE/GW accurately reproduces the trends for all series with excellent correlation coefficients to the benchmark data and generally very small errors. Indeed, for Thiel’s sets, the BSE/GW values are more accurate (using CC3 references) than both CC2 and ADC(3/2) values on both absolute and relative scales. For all three sets, BSE/GW errors also tend to be nicely spread with almost equal numbers of positive and negative deviations as compared to reference values.

1. INTRODUCTION The panel of theoretical models available for determining electronically excited-state energies and properties of molecules has become very broad. Among the available single-reference ab initio methods,1 one can distinguish three groups of techniques. The first relies on post-Hartree−Fock wave function methods and mainly encompasses variations of symmetry-adapted cluster configuration interaction (SAC−CI ),2,3 second-order polarization propagator approximation (SOPPA),4,5 coupled cluster (CC),6−12 and algebraic diagrammatic construction (ADC) approaches.13−16 These methods present the advantage of being systematically improvable as one can increase the correlation order (e.g., CC2, CCSD, CC3, CCSDT, etc.) to check the convergence of the results. However, they may also exist in several forms, e.g., linear-response (LR)6 and equationof-motion (EOM)7 formulations of CC theory have been proposed. Although the two theories yield (for a given order) the same transition energies, they yield different oscillator strengths with the EOM-CC values being approximations of their LR-CC counterparts. In addition, these wave function approaches are also computationally involved, i.e., the two leastdemanding and most popular schemes of the series, namely CC2 and ADC(2), present a formal 6(N 5) scaling with system size. The second group of models contains extension of ground-state density functional theory (DFT) to excited states. In that group, the adiabatic formulation of time-dependent DFT (TD-DFT)17,18 is certainly the most widely applied method. The formal scaling of TD-DFT is advantageous, and © XXXX American Chemical Society

very efficient implementations have been made allowing treating large systems in a variety of realistic environments.19 However, the quality of adiabatic TD-DFT’s predictions significantly depends on the selected exchange-correlation functional,20 and some families of excited-states remain difficult to capture, e.g., the states with a strong double-excitation character. As a third technique, the Bethe−Salpeter equation (BSE) formalism,21−24 a Green’s function perturbation theory, presents the same scaling as TD-DFT in its standard adiabatic implementation. In this approach, one usually first corrects the energies of input Kohn−Sham DFT eigenstates with the GW many-body perturbation technique24−27 and next determines the transition energies between these corrected levels at the BSE level. To date, the total number of applications of BSE/ GW for organic molecules remains rather limited28−57 as compared to the more than 1000 publications per year appearing in the TD-DFT field since 2011. Although it is still early to draw a definitive picture of all the pros and cons of BSE/GW compared to those of TD-DFT, recent extended benchmarks performed on the well-known Thiel’s set of molecules58,59 demonstrated an agreement with reference calculations comparable to the best TD-DFT calculations provided that the input GW energies are calculated with accuracy.54,55 Further, the BSE/GW formalism yields accurate results for systems known to be problematic within TD-DFT, such as molecules or complexes presenting charge-transfer Received: April 25, 2016

A

DOI: 10.1021/acs.jctc.6b00419 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation excitations34,36,38,39,41,45 and cyanine derivatives,46,47 which concerning the last family hints that the standard “adiabatic” BSE performs better than (adiabatic) TD-DFT for states with a double-excitation character.46,47,55,56 On the contrary, BSE seems to undergo some difficulties (too low excitation energies) for tiny molecules, probably due to a self-screening error.52,55,60 Much less is known about the accuracy of the BSE/GW oscillator strength, f, though this characteristic is also extremely important for interpreting spectra as well as for designing new compounds. In contrast, there has been several benchmarks of f for other excited-state theories, notably, for TD-DFT.20 At this level, assessments performed on compact molecules (∼2−15 atoms) concluded that global hybrid functionals with a large share (∼50%) of exact exchange or range-separated hybrids provided the most accurate results.61,62 This conclusion is in line with the results obtained by Bartlett’s group that concluded that B3LYP provides f that are too small for DNA bases.63 The two most extensive analyses of the quality of TD-DFT’s f were carried out by Caricato and co-workers64 and Silva-Junior and co-workers.59 In the former, a set of 69 excited states determined on 11 small organic molecules was designed, and EOM-CCSD’s f were used as references. An extensive screening of DFT functionals led Caricato et al. to identify CAM-B3LYP and LC-ωPBE as the two most accurate functionals. In the latter, using Thiel’s set, the authors reported mean absolute errors of 0.075, 0.055, and 0.044 for the f determined with BP86, B3LYP, and BHLYP when using CC2/TZVP values as the benchmark.65 However, the same authors have also reported, for the same set of compounds, significant discrepancies between CAS-PT2 and CC2 oscillator strengths.58 The recently proposed SOPPA analysis of oscillator strengths also show some significant deviations for the same set of compounds.66 This challenge of obtaining a highly accurate reference f was recently confirmed by Kannar and Szalay who determined CC3/TZVP f for the same set of compounds.67 Indeed, they found that the CC2 f tends to be significantly too large with an average deviation of 0.030 (or 28%) for the states for which nonzero CC2 and CC3 f were computed. Similarly, Harbach, Wormit, and Dreuw, choosing again Thiel’s set, investigated several variations of the ADC method.68 They concluded that the ADC(3/2) f generally nicely match the reference values and that ADC(2) also deliver rather accurate f but for some significant deviations w.r.t. ADC(3/2) for the states involving double-excitation contributions. Our goal here is to provide a first assessment of the accuracy of BSE/GW oscillator strengths. To this end, we determined ∼200 f for compounds divided into three molecular sets allowing us to draw relevant conclusions. In the first set, we considered Thiel’s set of compounds and used CC3/TZVP67 and ADC(3/2)/TZVP68 f values as references. This set contains high-lying excited-states, but the size of the treated molecules is limited (naphthalene and the nucleobases are the largest derivatives of the set). In the second set, we considered the 80 “real-life” dyes taken from ref 56 and compared the f of the first dipole-allowed excited states computed with BSE/GW, ADC(2), and CC2, all relying on the aug-cc-pVTZ atomic basis set. Eventually, in the third set, we considered 30 substituted anthraquinones for which experimental f values have been determined.69 Together, these three sets provide an overview of the accuracy of BSE/GW’s oscillator strengths.

2. COMPUTATIONAL DETAILS We provide in the Appendix a brief description of the Bethe− Salpeter formalism with explicit definition of the computed Scheme 1. 9,10-Anthraquinone with Numbering for Substitution Positions

oscillator strengths. More details can be found in reviews of the general formalism21−25 even though the literature on BSE/GW oscillator strength is very scarce.53 Further, the Bethe−Salpeter equation (BSE) approach also formally resembles the standard TD-DFT implementation within Casida’s formulation70 with the same difficulties when going beyond the Tamm−Dancoff approximation that the eigensystems to be inverted are not hermitian. As such, besides the specific form of the BSE resonant and nonresonant blocks involving the screened Coulomb potential, most of the algebra can also be found in standard reviews on TD-DFT.71 2.1. Overview. All GW and BSE calculations presented below have been performed with the FIESTA package,36,72,73 a Gaussian-basis implementation of the GW and BSE formalisms using the “Coulomb-fitting” resolution-of-identity (RI-V) approach and contour deformation techniques. Following previous benchmarks for the excitation energies on Thiel’s set55 and a larger set of 80 “real-life” molecules,56 our GW calculations are performed at the partially self-consistent level, namely, self-consistently updating the occupied/virtual GW energy levels while keeping the input Kohn−Sham wave functions frozen. Such a simple scheme, relying on the observation that it allows the dramatic reduction of the dependency on the starting Kohn−Sham eigenstates and provides much improved occupied/virtual energy levels72,74 as compared to those of a single-shot perturbative G0W0 approach, where the Green’s function G0 and screened Coulomb potential W0 are built from the input Kohn−Sham eigenstates, which results in significantly too small G0W0 HOMO−LUMO gaps in particular when starting from a semilocal DFT functional. As shown in the Appendix, the quality of the GW quasiparticle energies directly condition the BSE Hamiltonian through the occupied-to-virtual electronic energy level differences (Appendix eq 15) and through the screened Coulomb potential W built from the irreducible polarizability (Appendix eq 10). It was demonstrated in previous studies36,54,55 that, in the case of molecular systems, excellent (singlet) BSE excitation energies can be obtained when starting from self-consistent evGW calculations, namely, a GW scheme providing reliable quasiparticle energies. This is an important conclusion showing that the electron−hole binding energy is accurately calculated, namely, that the good performances of the BSE excitations energies do not benefit from error cancellations between the quasiparticle gap and the excitonic binding energy. Further, because the BSE Hamiltonian derives from the GW self-energy, the calculated excitation energies depend very weakly on the starting DFT Kohn−Sham eigenstates when starting from an evGW quasiparticle spectrum as compared to a calculation starting from non self-consistent G0W0 calculations. As we will see below, the same holds for the BSE oscillator strengths. B


Article

Journal of Chemical Theory and Computation

Table 1. Comparison between the Oscillator Strengths Obtained with the BSE/G0W0 and the BSE/evGW Approaches Starting with Four Different DFT Eigenstatesa BSE/G0W0@DFT molecule E-butadiene cyclopentadiene benzene s-triazine

BSE/evGW@DFT

state

M06-L

M06

M06-2X

M06-HF

M06-L

M06

M06-2X

M06-HF

CC

Bu B2 A1 E1u A″2 E′

0.548 0.070 0.438 0.465 0.009 0.200

0.618 0.081 0.508 0.549 0.011 0.338

0.654 0.086 0.555 0.580 0.017 0.357

0.704 0.096 0.625 0.642 0.019 0.470

0.681 0.084 0.505 0.540 0.015 0.329

0.704 0.087 0.556 0.592 0.016 0.371

0.692 0.090 0.575 0.597 0.017 0.404

0.698 0.095 0.619 0.633 0.014 0.401

0.725 0.092 0.563 0.630 0.016 0.386

1 1 3 1 1 1

a

All results are obtained with the TZVP atomic basis set. In the rightmost column, we give the reference CC values taken from ref 67. See the caption and footnotes of Table 2 for more details regarding these values.

aug-cc-pVTZ atomic basis sets,55 and we do not further discuss the accuracy of these transition energies here. 2.3. Real-Life Dyes. The second set of compounds is constituted of 80 “real-life” molecules. The representation of these molecules as well as all Cartesian coordinates can be found in ref 56. In this set, the ground-state geometries have been obtained at the M06-2X/6-31+G(d) level of approximation using Gaussian0981 with a tight optimization threshold (1 × 10−5 on residual rms forces) and the ultraf ine DFT integration grid. The f have been determined with the aug-ccpVTZ atomic basis set at the BSE/evGW level. These values are compared to TD-M06-2X/6-311++G(2df,2p), ADC(2)/augcc-pVTZ, and CC2/aug-cc-pVTZ results, considering the latter as reference. The TD-DFT values were determined with Gaussian09,81 whereas both ADC(2) and CC2 values have been computed with Turbomole.82 We only considered the first dipole-allowed state for each compound. The corresponding transition energies have been analyzed in ref 56. 2.4. Anthraquinones. The third set consists of 30 9,10anthraquinone dyes (see Scheme 1) for which experimental oscillator strengths have been determined by Labhart in dichloromethane.69 Though this anthraquinone set has been modeled previously with TD-DFT,83 there is, to the best of our knowledge, no previous BSE/GW calculations for this family of compounds. The optimal ground-state geometries have been determined with Gaussian0981 at the M06-2X/6-311G(d,p) level using a tight convergence threshold and the ultraf ine integration grid. The BSE/evGW calculations have been performed with the aug-cc-pVDZ atomic basis set considering 5 evGW iterations and correcting 40 DFT levels around the gap (20 occupied and 20 virtuals). These parameters have been shown previously to be suited for conjugated organic compounds.57

Concerning the RI-V approximation, our calculations have been performed using Weigend’s RI auxiliary bases.75 A complete presentation of the theories used can be found in ref 55. Here, the DFT starting eigenstates used to perform these GW calculations have been determined with the NWChem package76 selecting the M06-2X exchange-correlation functional.77 This choice is justified by previous investigations54,57 that showed that it ensures both a maximal accuracy and a rapid convergence of the transition energies. The same holds for oscillator strengths as demonstrated below. The NWChem DFT calculations relied on the xf ine integration grid with total energies and densities converged at least to 10−7 a.u. and 10−6 a.u., respectively. Further computational details specific to each set of molecules are given below. All reference CC2, CCSD, and CC3 f reported in this work have been obtained from the literature considering the linearresponse (LR) formalism.58,65,67 We recall, however, that the differences between LR f and their computationally more efficient EOM counterparts are generally small.67 As the computed f vary over several orders of magnitude, it is not straightforward to provide a relevant statistical analysis. To circumvent this problem, as already detailed by other authors,64,67 we report both deviations of the absolute and relative f in the following. With the former approach, the weight is given on the states or compounds presenting large f, whereas the latter, which considers the relative errors (in %) compared to the reference data, gives equal importance to all f irrespective of their amplitudes (but may lead to large deviations when f is negligible). In the text and tables below, we report both the mean signed error (MSE) and mean absolute error (MAE) in these two units. 2.2. Thiel’s Set. The first series of compounds is composed of the well-known Thiel’s set of molecules encompassing small organic conjugated molecules and DNA bases.58,65,78 The geometries used, obtained at the MP2/6-31G(d) level, can be found in ref 58. TZVP was selected as the atomic basis set, a choice made to ensure consistency with the reference ADC(3/ 2)68 and CC3 values,67 and also to avoid strong state mixing with the (normally high-lying) Rydberg states. Such mixing, typical of strongly diffuse atomic basis sets,79,80 makes the assignments more complex. Nevertheless, as CC2/aug-ccpVTZ f are available for this set of compounds,65 we have performed a series of test BSE/evGW calculations with aug-ccpVTZ as well. Consistent with our recent work,55 we corrected all valence DFT levels (and an equivalent number of virtual levels) with the evGW approach, and 5 evGW iterations were typically used to reach convergence. For this set of molecules, the interested reader will find elsewhere BSE/GW benchmarks of the transition energies performed with both TZVP54 and

3. RESULTS AND DISCUSSION 3.1. First Set: Thiel’s Series. For four molecules (six transitions), we have compared the results obtained with the BSE/G0W0 and BSE/evGW approaches starting with four different DFT eigenstates determined with the M06 series of exchange-correlation functionals: M06-L (0% of exact exchange), M06 (27%), M06-2X (54%), and M06-HF (100%). We recall that the G0W0 approach is a pertubative method such that a significant dependency on the starting DFT functional is expected. In contrast, in the evGW scheme, the eigenvalues are self-consistently converged such that only a minor dependency on the frozen eigenfunctions pertains. The results are given in Table 1. Paralleling the conclusions obtained previously for transition energies,55,57 we note that (i) BSE/G0W0 oscillator strengths (f) tend to be too small and are significantly sensitive C


Article


Table 2. Transition Energies (ΔE in eV) and Oscillator Strengths ( f) Determined for Thiel’s Set at the BSE/evGW@M06-2X/ TZVP Levela reference f values

BSE/evGW molecule ethene E-butadiene all-E-hexatriene all-E-octratetraene cyclopropene cyclopentadiene

norbornadiene benzene naphthalene

furan

pyrrole

imidazole

pyrdine

pyrazine

pyrimidine

pyridazine

s-triazine s-tetrazine

formaldehyde

state 1 1 1 1 1 1 1 2 3 1 2 1 1 2 2 3 1 2 3 2 1 3 2 1 3 1 1 2 2 3 4 3 4 1 1 1 2 2 1 1 2 3 2 1 2 2 1 2 3 1 1 1 1 2 1 2 2 1

B1u (π → π*) Bu (π → π*) Bu (π → π*) Bu (π → π*) B2 (π → π*) B1 (σ → π*) B2 (π → π*) A1 (π → π*) A1 (π → π*) B2 (π → π*) B2 (π → π*) E1u (π → π*) B2u (π → π*) B3u (π → π*) B2u (π → π*) B2u (π → π*) B2 (π → π*) A1 (π → π*) A1 (π → π*) A1 (π → π*) B2 (π → π*) A1 (π → π*) A′ (π → π*) A″ (n → π*) A′ (π → π*) B1 (n → π*) B2 (π → π*) A1 (π → π*) B2 (π → π*) A1 (π → π*) A1 (π → π*) B2 (π → π*) B2 (π → Ryd.) B3u (n → π*) B2u (π → π*) B1u (π → π*) B2u (π → π*) B1u (π → π*) B1 (n → π*) B2 (π → π*) A1 (π → π*) A1 (π → π*) B2 (π → π*) B1 (n → π*) A1 (π → π*) B1 (n → π*) B2 (π → π*) B2 (π → π*) A1 (π → π*) A″2 (n → π*) E′ (π → π*) B3u (n → π*) B2u (π → π*) B3u (n → π*) B1u (π → π*) B1u (π → π*) B2u (π → π*) B1 (σ → π*)

ΔE

f

CC3

ADC(3/2)

7.847 6.034 5.025 4.383 6.748 6.729 5.194 7.112 8.312 6.308 7.337 7.024 4.536 5.918 6.200 8.198 6.260 6.936 8.291 6.653 6.471 8.003 6.573 6.854 7.124 5.169 5.355 6.404 7.263 7.290 9.373 9.613 9.448 4.241 5.164 6.565 7.747 7.612 4.579 5.587 6.679 7.439 7.643 3.927 5.434 6.540 6.490 7.167 7.454 4.895 7.727 2.420 5.273 6.790 6.883 7.337 8.109 8.668

0.360 0.692 1.096 1.528 0.079 0.001 0.090 0.013 0.575 0.037 0.192 0.597 0.076 1.269 0.248 0.463 0.161 0.002 0.448 0.007 0.173 0.468 0.154 0.004 0.023 0.005 0.036 0.017 0.481 0.514 0.001 0.013 0.126 0.007 0.090 0.066 0.359 0.395 0.006 0.034 0.044 0.412 0.435 0.006 0.023 0.006 0.007 0.440 0.425 0.017 0.404 0.006 0.061 0.012 0.001 0.353 0.274 0.107

0.389b 0.725b 1.129 1.549 0.082b 0.001b 0.092b 0.004b 0.563b 0.027 0.185 0.630 0.085 1.325 0.239 0.498 0.155 0.001 0.450 0.004 0.167 0.478 0.081 0.004 0.082 0.005 0.021 0.014 0.482 0.526 0.004 0.030 0.114 0.007 0.062 0.070 0.376 0.407 0.006 0.021 0.043 0.391 0.415 0.006 0.014 0.005 0.012 0.340 0.433 0.016 0.386 0.007 0.044 0.011 0.002 0.349 0.307 0.003b

0.423 0.806 1.257 1.724 0.096 0.001 0.101 0.001 0.029c 0.036 0.221 0.707 0.097 1.524 0.299

D

0.160 0.000 0.345 0.004 0.171 0.458 0.112 0.000 0.064 0.005 0.027 0.010 0.538 0.620

0.007 0.077 0.070 0.446 0.496 0.007 0.028 0.025 0.462 0.523 0.007 0.019 0.005 0.012 0.492 0.466 0.016 0.443 0.008 0.000d 0.011 0.004 0.421 0.013 0.091


Article

Journal of Chemical Theory and Computation Table 2. continued reference f values

BSE/evGW molecule acetone p-benzoquinone formamide

acetamide

propanamide cytosine

thymine

uracil

adenine

ΔE

state 2 1 1 2 1 2 3 4 1 2 3 2 3 2 1 2 3 4 2 3 4 2 3 4 2 1 3 2 4

A1 (π → π*) A1 (π → π*) B1u (π → π*) B1u (π → π*) A″ (n → π*) A′ (π → π*) A′ (n → Ryd.) A′ (n → Ryd.) A″ (n → π*) A′ (π → π*) A′ (π → π*) A′ (π → π*) A′ (π → π*) A′ (π → π*) A″ (n → π*) A″ (n → π*) A′ (π → π*) A′ (π → π*) A′ (π → π*) A′ (π → π*) A′ (π → π*) A′ (π → π*) A′ (π → π*) A′ (π → π*) A′ (π → π*) A″ (n → π*) A′ (π → π*) A″ (n → π*) A′ (π → π*)

10.714 9.099e 5.269 7.875 5.675 7.593 8.607 8.862 5.695 7.902i 10.869 7.851j 10.465 4.795 5.234 5.922 5.769 6.636 5.294 6.579 6.858 5.412 6.532 6.998 5.326 5.388 5.231 6.011 6.501

f 0.308 0.116e 0.466 0.462 0.001 0.080 0.303 0.142 0.001 0.198i 0.269 0.160j 0.113 0.061 0.002 0.000 0.143 0.436 0.186 0.069 0.229 0.185 0.048 0.175 0.063 0.001 0.226 0.001 0.461

CC3 b

0.347 0.245 0.485 0.131f 0.001b 0.030b,g 0.388b,g 0.122b 0.001 0.207 0.263 0.170 0.136 0.046 0.001 0.001 0.130 0.520 0.172 0.072 0.197 0.174 0.046 0.152 0.002k 0.001 0.297k 0.002 0.513k

ADC(3/2) 0.096 0.134 0.590 0.001 0.135h 0.299 0.001 0.132 0.145 0.126 0.112 0.071 0.002 0.000 0.143 0.538 0.238 0.055 0.242 0.242 0.040 0.179 0.006 0.001 0.274 0.002 0.249l

a

CC3/TZVP values were taken from ref 67, and ADC(3/2)/TZVP values were taken from ref 68. bCCSDT/TZVP value. cADC(2)-s value: 0.695, inconsistency probably due to switching between 3 A1 and 4 A1 states. dADC(2)-s value: 0.055. eAnother closed state (9.830 eV) presents f of 0.169. f CCSD value: 0.506. gStrongly mixed states. hADC(2)-s value: 0.373. iMixed with a closed state at 7.540 eV (f = 0.034). jMixed with a closed state at 7.534 eV ( f = 0.039). kCCSD/TZVP value. lApparently mixed with the 5 A′ state presenting f = 0.257; ADC(2)-s f for 4 A′ and 5 A′ are 0.563 and 0.044, respectively.

Figure 1. (left) Comparison between BSE/evGW f and CC reference values for Thiel’s set of compounds (81 f considered, see text). The central black line indicates a perfect match, and the dashed red line is obtained through simple linear regression. (right) Relative BSE/evGW error in % as a function of the CC f for the same set of data.

to the selected starting eigenstates; (ii) the BSE/G0W0@M06-L values are especially far from both convergence and reference values, highlighting that such calculations can provide very inaccurate results; (iii) the BSE/evGW oscillator strengths are much less affected by the starting DFT eigenstates than their

BSE/G0W0 counterparts, as only small variations related to the frozen eigenfunctions pertain; (iv) the results obtained starting with the two standard hybrids, i.e., BSE/evGW@M06 and BSE/evGW@M06-2X, are extremely similar; and (v) M06-2X is probably the most adequate starting point as it allows a fast E


Article


Table 3. Results of Statistical Analysis Performed for the BSE/evGW f of Table 2 Considering Various References and Subgroups of Statesa absolute values

relative values 2

ref

state

MSE

MAE

R

ADC(3/2)

all (81) π → π* (63) n → π* (15) “safe states 20%”d (40) all (70)

−0.000 −0.003 0.000 −0.009 −0.018

0.018 0.020 0.000 0.014 0.014

0.991 0.992 0.981 0.997 0.970

CC

MSE (%)

MAE (%)

6.7b 8.5c −1.9 −1.5 1.1e

22.5b 23.4c 20.0 9.6 22.6e

a

The states displaying strong state-mixing were not considered during this analysis. bObtained removing two clear outliers (1 B1 of formaldehyde and 2 A′ of adenine) presenting >1000% difference (see Table 2). cObtained removing one clear outlier (2 A′ of adenine). dStates for which the absolute deviation between the reference CC and ADC(3/2) f does not exceed 20%. eObtained removing three clear outliers (2 A1 of cyclopentadiene, 2 B2u of s-tetrazine, and 2 A′ of adenine) presenting 1200, 2007, and 950% differences, respectively (see Table 2).

Table 4. Basis Set Effects for the f of Selected States in Thiel’s Seta BSE/evGW molecule ethene E-butadiene all-E-hexatriene all-E-octratetraene norbornadiene benzene furan

imidazole s-triazine p-benzoquinone a

state 1 1 1 1 1 2 1 1 2 3 2 1 1 1

B1u (π → π*) Bu (π → π*) Bu (π → π*) Bu (π → π*) B2 (π → π*) B2 (π → π*) E1u (π → π*) B2 (π → π*) A1 (π → π*) A1 (π → π*) A′ (π → π*) A″2 (n → π*) E′ (π → π*) B1u (π → π*)

CC2

TZVP

AVTZ

ΔBS

TZVP

AVTZ

ΔBS

0.360 0.692 1.096 1.528 0.037 0.192 0.597 0.161 0.002 0.448 0.154 0.017 0.404 0.462

0.318 0.646 1.051 1.475 0.034 0.123 0.567 0.153 0.000 0.405 0.155 0.013 0.267 0.440

−11.7 −6.6 −4.1 −3.5 −8.1 −35.9 −5.3 −5.0 −100.0 −9.6 +0.6 −23.5 −33.9 −4.8

0.431 0.809 1.272 1.757 0.023 0.185 0.694 0.172 0.003 0.506 0.088 0.017 0.441 0.538

0.375 0.734 1.246 1.749 0.020 0.083 0.664 0.184 0.000 0.358 0.112 0.015 0.410 0.530

−13.0 −9.3 −2.0 −0.5 −13.0 −55.1 −4.3 +6.9 −100.0 −29.2 −27.3 −11.8 −7.0 −1.5

All CC2 values are taken from ref 65. AVTZ stands for aug-cc-pVTZ. ΔBS gives the impact of going from TZVP to aug-cc-pVTZ, expressed in %.

ethene and ∼1% for octratetraene. For the other conjugated hydrocarbons, the only significant discrepancy is obtained for the 2 A1 state of cyclopentadiene that is weakly allowed. In the aromatic and heterocycles series, one also notices the good accuracy of BSE/evGW data with only three states for which rather large discrepancies could be noted (2 and 3 A′ of imidazole and 3 B2 of pyridine). The BSE/evGW values statistically spread around the CC3 values (22/45 larger, 20/45 smaller, 3/45 equal). In the aldehyde, ketone, and amide groups, it is difficult to obtain definitive conclusions due to the limited number of clear-cut cases, but the f of most of the states not suffering from strong state-mixing with any method are accurately described by BSE/evGW. For instance, for the B1 excited-state of formaldehyde, the BSE/evGW prediction (0.107) is in good agreement with the ADC(3/2) value (0.091), but both are in sharp contrast with the CC3 result (0.003). Eventually, for the nucleobases, BSE/evGW f presents the correct order of magnitude except for one clear exception (the 2 A′ state of adenine). To perform our statistical analysis, we have only considered states not evidently influenced by significant state mixing (see footnotes in Table 2 for problematic cases), which led to 81 f for CC and 70 f for ADC(3/2). A graphical comparison between BSE/evGW and reference CC f can be found in Figure 1, and the results of the statistical treatment are given in Table 3. A first striking result is the very good correlation between the BSE/evGW and reference results with linear coefficients of determination, R2, of 0.98 or 0.99 for n → π* and π → π*

convergence, as illustrated by the relatively small differences between the BSE/evGW and BSE/G0W0 oscillator strengths determined starting with this functional. For this small subset, we also notice that the TD-DFT oscillator strengths are not strongly affected by the selected functional, but this question is beyond our scope here. Below, we consider the BSE/evGW@ M06-2X data only. The oscillator strengths determined for the complete Thiel’s set of compounds can be found in Table 2 together with the reference values obtained at both ADC(3/2)68 and CC67 levels using the same atomic basis set. For the latter, we have selected the most refined f determined in ref 67, that is, mostly CC3/ TZVP results but for the smallest molecules (CCSDT/TZVP) and for the A′ states of adenine (CCSD/TZVP). We have, of course, not listed states forbidden by symmetry that present null f irrespective of the theoretical model. Likewise, dipoleallowed states presenting f < 1 × 10−3 for all methods, e.g., the lowest B3u state of p-benzoquinone, have been ignored. Eventually, six high-lying states listed in ref 67 were not reported as no unambiguous correspondence with BSE/evGW results could be made, typically due to strong state mixing.84 From the results of Table 2, one first notices a reasonable agreement between the CC3 and ADC(3/2) trends except for a few cases where strong state mixing occurs or for which ADC(2) and ADC(3/2) f do not match (see Table 2 footnotes and discussion below). In the conjugated olefin series, BSE/ evGW provides f slightly smaller than the CC reference values, though the underestimation remains rather limited, ∼7% for F


Article

Journal of Chemical Theory and Computation Table 5. Oscillator Strengths Determined for the 80 Compounds of Ref 56a molecule

TD

ADC(2)

CC2

BSE

molecule

TD

ADC(2)

CC2

BSE

I II III IV V VI VII VIII IX X XI XII XIII XIV XV XVI XVII XVIII XIX XX XXI XXII XXIII XXIV XXV XXVI XXVII XXVIII XXIX XXX XXXI XXXII XXXIII XXXIV XXXV XXXVI XXXVII XXXVIII XXXIX XLc

0.469 0.498 0.357 0.927 0.206 0.214 1.261 0.244 0.049 0.393 1.539 0.070 1.003 0.468 1.501 0.427 0.409 1.250 0.892 0.258 1.312 1.021 0.353 0.438 0.777 0.818 1.989 0.738 0.920 2.059 0.314 0.613 0.390 0.408 0.821 0.249 0.930 0.066 0.066 0.505

0.599 0.515 0.344 0.914 0.235 0.246 1.247 0.204 0.043 0.445 1.277 0.067 0.741 0.538 1.643 0.297 0.438 1.007 0.802 0.288 1.396 1.231 0.346 0.381 0.864 0.799 2.058 0.761 0.640 2.135 0.258 0.691 0.340 0.449 0.683 0.276 0.883 0.049 0.063 0.734

0.547 0.499 0.337 0.891 0.266 0.239b 1.313 0.247 0.038 0.381 1.792 0.062 0.960 0.504 1.683 0.498 0.431 1.248 0.965 0.235 1.351 1.145 0.326 0.454 0.789 0.857 2.214 0.743 0.798 2.222 0.338 0.680 0.326 0.402 0.859 0.234 0.890 0.057 0.067 0.440

0.443 0.412 0.325 0.874 0.189 0.182 1.223 0.201 0.043 0.319 1.270 0.058 1.043 0.474 1.481 0.379 0.328 1.116 0.757 0.242 1.300 1.000 0.316 0.392 0.601 0.651 1.991 0.622 0.880 2.140 0.273 0.579 0.339 0.348 0.733 0.338 0.882 0.051 0.062 0.411

XLI XLII XLIII XLIV XLV XLVI XLVII XLVIII XLIX L LI LII LIII LIV LV LVI LVII LVIII LIX LX LXI LXII LXIII LXIVc LXV LXVI LXVII LXVIII LXIX LXX LXXI LXXII LXXIII LXXIV LXXV LXXVI LXXVII LXXVIII LXXIX LXXX

0.981 0.945 1.487 0.797 0.072 0.064 0.616 0.873 0.044 1.716 1.484 0.628 0.602 0.427 0.433 0.657 0.671 0.310 0.143 1.180 0.683 0.800 0.693 0.365 0.493 1.024 0.049 0.106 0.503 0.159 0.210 1.721 1.333 0.628 1.001 0.334 0.303 1.451 0.898 0.603

1.192 1.222 1.703 0.861 0.085 0.077 0.696 0.979 0.043 1.752 1.250 0.541 0.522 0.445 0.433 0.614 0.493 0.342 0.098 1.327 0.771 0.910 0.692 0.492 0.439 1.076 0.050 0.108 0.374 0.121 0.179 1.543 1.186 0.533 0.756 0.349 0.310 1.428 0.601 0.553

1.130 1.186 1.694 0.809 0.071 0.062 0.653 0.877 0.039 1.841 1.602 0.675 0.658 0.445 0.446 0.725 0.483 0.325 0.095 1.333 0.761 0.882 0.698 0.336 0.542 1.221 0.043 0.107 0.298 0.099 0.161 1.782 1.340 0.520 0.940 0.324 0.288 1.506 0.883 0.661

0.955 0.929 1.461 0.680 0.060 0.052 0.527 0.768 0.037 1.736 1.249 0.548 0.522 0.344 0.353 0.638 0.611 0.257 0.124 1.199 0.661 0.774 0.669 0.322 0.458 1.238 0.042 0.082 0.446 0.136 0.181 1.635 1.330 0.598 1.041 0.289 0.256 1.440 0.741 0.580

a We used the same numbering as in Schemes 1−4 of that previous work and considered absorption from the ground-state geometries. TD stands for TD-M06-2X, and BSE stands for BSE/evGW@M06-2X. All calculations are with aug-cc-pVTZ but the TD-DFT relies on 6-311++G(2df,2p). bCC2 yields two nearly degenerate states (4.111 and 4.117 eV); the sum of the two f is given. cDegenerate B states with f corresponding to one state given.

deviate substantially from the CC3 values for the same set of compounds.67 Indeed, comparing the CC2 and CC3 f of ref 67, we determined an MSE of 0.030 (28.0%) and an MAE of 0.033 (32.8%). In other words, BSE/evGW results are significantly closer to the CC3 references than their CC2 counterparts that tend to be too large. Likewise, by comparing the two highly correlated methods of Table 2, CC3 and ADC(3/2), a group of 72 common f (with no strong mixing, nor states presenting a null f) could be defined. For that set, ADC(3/2) slightly overestimates the CC3 results (MSE of 0.019 or 2.2%), but the absolute deviation remains non-negligible (MAE of 0.042 or 27.5%, again larger than the BSE/evGW values), highlighting the difficulty of obtaining very accurate oscillator strengths. This is why we also considered the series of “safe states” in Table 3, that is, a group for which the CC and ADC(3/2) oscillator strengths are within ±20% of each other. For this

Table 6. Results of a Statistical Analysis Performed for the Molecules of Table 5 Considering the CC2 Values As Reference absolute values

relative values

method

MSE

MAE

R2

MSE (%)

MAE (%)

TD-M06-2X ADC(2) BSE/evGW

−0.027 −0.035 −0.070

0.060 0.077 0.092

0.983 0.958 0.969

1.6 −0.5 −7.9

10.5 12.0 14.5

states, respectively. Next, the trends obtained when using CC or ADC f as reference are similar. It is also obvious that the BSE/evGW f are very close to the reference values for the majority of states with an MSE of −4 × 10−4 (+6.7%) and an MAE of 0.018 (22.5%) with respect to CC. If the relative error percentages might appear quite large, we recall that CC2 f also G


Article


Figure 2. (left) Comparison between TD-M06-2X (top), ADC(2) (center), and BSE/evGW (bottom) f for the set of molecules of Table 5. CC2 is used as reference. See the caption of Figure 1 for more details.

subset containing 40 f, one also finds that BSE/evGW is extremely accurate with MAE of 0.014, or 9.6%, with respect to CC3 and an R2 very close to unity. A possible source of error in the previous comparison is the use of the TZVP atomic basis set, which could be too compact, especially for coupled-cluster calculations. To provide the first insights into the role of the atomic basis set, we have determined the f with aug-cc-pVTZ for a subset of states. The results are listed in Table 4. As no CC3 values are available with this atomic basis set, we have used CC2/aug-cc-pVTZ f obtained in ref 65. From Table 4, one indeed notices significant changes (mostly decreases) of the f when increasing the size of the basis set, though most tested cases display absolute variations smaller than 15%. Except for a few states (3 A1 state of furan, 2 A′ of imidazole, and 1 E′ of s-triazine), there is also a reasonable agreement between the trends observed for both methods: the states that are the most sensitive to basis set effects are similar at the CC2 and BSE/evGW levels. As one exception, going from TZVP to aug-cc-pVTZ substantially

degrades the agreement between BSE/evGW and CC2 f for the (1 E′) state of s-triazine. For the states of Table 4, the MSE (MAE) between BSE/evGW and CC2 attains −0.056 (0.069) with TZVP and −0.060 (0.080) with aug-cc-pVTZ, whereas the average change of f induced by enlarging the basis set amounts to −0.036 (BSE/evGW) and −0.033 (CC2) hinting that the conclusions obtained above should not strongly vary when extending the atomic basis set. 3.2. Second Set: Real-Life Compounds. Let us now turn toward larger molecules closer to the typical target of BSE/ evGW calculations. For the set of 80 compounds of ref 56, the f obtained with four theoretical approaches are listed in Table 5. Table 6 provides the results of the corresponding statistical analysis, and Figure 2 gives a graphical representation of the obtained results. As we considered only dipole-allowed states in quite large π-conjugated molecules, the average f are significantly larger than in the previous set. Indeed, they range from 0.037 to 2.140 with a mean of 0.640. However, the size of the molecules implies the use of CC2 reference values H


Article


of M06-2X that contains 54% of exact exchange, a ratio close to ideal for obtaining accurate oscillator strengths with TD-DFT (see Introduction). More expected is the good match between ADC(2) and CC2 f that deviate by only 12% on average. The ADC(2) values are nevertheless a bit more spread around the CC2 references than the TD-M06-2X data, explaining the smaller R2 of the former approach in Table 6. Similarly to their TD-DFT counterparts, BSE/evGW f are too large for the smallest f but become too small when the CC2 f become large (see bottom right panel in Figure 2). The general correlation between BSE/evGW and CC2 f remains excellent (R2 = 0.969), and the trend of BSE/evGW to yield smaller f than CC2 yields a rather small MSE and MAE of −7.9% and 14.5%, respectively. This error is clearly mitigated by the fact that the CC2 values are too large by ∼10−20% (see above). In short, BSE/evGW provides accurate f for this set of large molecules judged on the expected good performances of both CC2 and ADC(2). 3.3. Third Set: Substituted Anthraquinones. The results obtained for anthraquinones are listed in Table 7 and compared to the experimental data of Labhart69 in Figure 3. First, we note that the vertical transition energies are larger (by ∼0.2−0.4 eV) than the energies corresponding to the measured absorption wavelengths. This discrepancy can be partly explained by the neglect of both solvation and vibronic effects that have been modeled elsewhere for anthraquinones.85 Second, for both the 1-hydroxy and 2-hydroxy-anthraquinones, we have computed the oscillator strengths of the lowest state with aug-cc-pVTZ and obtained f = 0.135 and 0.039, respectively. These values closely match their aug-cc-pVDZ counterparts (0.133 and 0.039, respectively), confirming that the diffuse-polarized double-ζ atomic basis set is sufficient for our purposes. Third, the f measured for anthraquinones increase regularly with the number of hydrogen bonds formed with the ketone (1, 4, 5, and 8 substitution with amines and hydoxyls), a trend correctly restored by theory. Figure 3 compares the BSE/evGW f with experimental values. At least two of the f reported to be less accurate experimentally69 are indeed clear outliers. For this reason, the five “inaccurate” f have been removed from the statistical analysis. As in the other sets, BSE/evGW provides a balanced description of oscillator strengths. For the 25 accurate f reported in the experimental paper, we found an MSE of +0.005 (−2.07%) and an MAE of 0.020 (15.7%) with a similar number of dyes for which the BSE/evGW f are smaller (14/25) or larger (11/25) than their experimental counterparts. Eventually, as can be seen in the left panel of Figure 3, the correlation between experimental and theoretical values remains satisfactory (R2 = 0.819). It is less impressive than in the two other sets of compounds but more sources of error play a role here (geometry, solvent effects, etc.).

Table 7. Transition Energies (ΔE in eV) and Oscillator Strengths ( f) Determined at the BSE/evGW@M06-2X/augcc-pVDZ Level for Substituted Anthraquinones (see Scheme 1)a BSE/evGW subst. 2-F none 2-Cl 2,3-Cl 1,8-Cl 1,4-Cl 2-OMe 2-OH 1-OMe 1-OH 2-NH2 1,2-OH 1,5-OH 1,8-OH 1-NH2 2-NMe2 1,4-OH 1,5-NH2 1,8-NH2 1-NHPh 1-NO2, 4,5,8-OH 1-OH, 4-NH2 1-OH, 2,4-NH2 1-NHMe, 4-OMe 1,4-NH2 1-NHMe,4-NH2 1,5-NH2, 4,8-OH 1-NH2, 4-NHPh 1,4-NHMe 1,4-NH2, 2-NO2

ΔE 4.277 4.291 4.245 4.226 4.050 3.913 3.720 3.822 4.048 3.413 3.478 3.292 3.307 3.270 3.064 3.201 2.836 3.008 2.914 2.871 2.714 2.624 2.620 2.837 2.517 2.382 2.388 2.386 2.257 2.174

experiment

f b

0.077 0.085c 0.094d 0.075e 0.112 0.092 0.042 0.039 0.088 0.133 0.061 0.110 0.243 0.246 0.131 0.088 0.167 0.250 0.226 0.171 0.206 0.167 0.158 0.156 0.173 0.184 0.306 0.229 0.193 0.175

λmax

ΔE

f

325 327 330 330 344 350 363 365 380 405 410 416 428 430 465 470 476 480 492 508 510 520 530 540 550g 590 590 590 620 645

3.815 3.792 3.757 3.757 3.604 3.542 3.416 3.397 3.263 3.061 3.024 2.980 2.897 2.883 2.666 2.638 2.605 2.583 2.520 2.441 2.431 2.384 2.339 2.296 2.254 2.101 2.101 2.101 2.000 1.922

0.095 0.094 0.102 0.102 0.102 0.119 0.060f 0.047f 0.099 0.113 0.091 0.119 0.191 0.225 0.114 0.109 0.131 0.206f 0.173 0.112 0.208 0.081f 0.175 0.142 0.157 0.192 0.139f 0.231 0.199 0.166

a

The experimental values (taken from ref 69) are given on the righthand-side: λmax in nm, the corresponding ΔE in eV and f. These experimental values have been measured in dichloromethane. bThree nearby states: 4.099 eV ( f = 0.010), 4.277 eV ( f = 0.046), and 4.394 eV ( f = 0.020); total f reported. cTwo nearly degenerate states at 4.291 eV ( f = 0.013) and 4.292 eV ( f = 0.072); total f reported. dThree nearby states: 4.088 eV ( f = 0.029), 4.245 eV ( f = 0.049), and 4.355 eV (f = 0.017); total f reported. eThree nearby states: 4.086 eV ( f = 0.006), 4.226 eV (f = 0.053), and 4.307 eV ( f = 0.016); total f reported. fValues given as estimated in the original experimental paper. g This wavelength, reported in ref 69, could correspond to the second maxima of the vibronic progression.

that are probably slightly too large as explained in the previous section. If we take Thiel’s set and limit the CC2/CC3 comparisons to states with f > 0.037 ( f > 0.200), we obtain relative MSE and MAE of 19.6 and 20.7% (11.2 and 11.9%), respectively, indicating that if the relative amplitude of the CC2 overestimation decreases with increasing f, it remains significant even for compounds presenting large oscillator strengths. From Figure 2 and Table 6, one notes that TD-M06-2X provides an astonishing match with CC2 with very large correlation and tiny MSE and MAE on both relative and absolute scales. Nevertheless, we note that TD-M06-2X has a tendency to slightly overestimate (underestimate) CC2 oscillator strengths for small (large) f (top right panel in Figure 2). This success of TD-DFT is likely due to the selection

4. CONCLUSIONS For the first time, we have provided an assessment of the accuracy of the oscillator strengths obtained with BSE/GW for a large number of molecules and states. To this end, we have considered three sets of molecules and used various reference values from the literature. Contrary to transition energies, for which a large number of highly accurate benchmark data can be found, obtaining very accurate oscillator strengths is not an easy task, as differences between various highly correlated wave function approaches can be non-negligible. First, we demonstrated that the partially self-consistent approach consisting in reinjecting the eigenvalues in the GW process until convergence is reached but keeping frozen the DFT I


Article


Figure 3. (left) Comparison between BSE/evGW f and their experimental counterparts for substituted anthraquinones. (right) Relative error in % as a function of the experimental f for the same data. The blue squares correspond to the data reported as inaccurate in the original experimental work.

where, e.g., 1 = (r1,t1) is a space-time coordinate, ρ, Uext, and G are the charge density, a local external perturbation, and the time-ordered one-body Green’s function, respectively, by the four-point response quantity

eigenfunctions is an efficient approach allowing for reducing the impact of the selected starting DFT functional on the BSE/GW oscillator strengths, whereas BSE/G0W0 f tend to be both significantly affected by the starting DFT eigenstates and too small if a pure (nonhybrid) exchange-correlation functional is selected. These conclusions are in line with those obtained previously for transition energies. Second, for Thiel’s sets for which both ADC(3/2) and CC3 f are available, we found that BSE/evGW estimates are very accurate. Indeed, the differences between BSE/evGW f and the CC3 reference values are significantly smaller than the discrepancies between CC3 and both ADC(3/2) or CC2 results, the latter model tending to provide slightly too large f. Third, for a set of 80 dyes (20−70 atoms), BSE/evGW provides oscillator strengths smaller (∼8%) than the CC2 values, whereas TD-M06-2X provides an extremely good match with CC2, keeping in mind that, at least for Thiel’s set, CC2 values overestimate their CC3 counterparts. Eventually, for 25 anthraquinones for which experimental f have been measured, BSE/evGW again yields accurate estimates with an equal number of positive and negative deviations and a mean absolute deviation limited to ∼0.020, or 15.7%. In short, despite the disparities of the three sets of molecules, a main conclusion can therefore be reached: BSE/evGW oscillator strengths are both consistent and accurate. Finally, we stress that, even if performing partially self-consistent GW calculations comes with an increased computational cost compared to both BSE/G0W0 and TDDFT, this 6(N 4) approach delivers accurate excited-state energies and oscillator strengths that are almost independent of the starting exchange-correlation functional. More works are certainly needed to assess the accuracy of BSE/evGW for other excited-state properties (dipole and quadrupole moments, structures, vibrational energies, etc.), but we see these first results as quite encouraging.

L(1, 2; 3, 4) = −i

∂G(1, 2) L̃(1, 2; 3, 4) ≡ −i ∂Utot(3, 4) L(1, 2; 3, 4) = L̃(1, 2; 3, 4) +

(3)

∫ d(5, 6)L̃(1, 2; 5, 5)v

(5, 6)L(6, 6; 3, 4)

(4)

In the standard implementation of the Bethe−Salpeter approximation, the polarizability L̃ is linked to the four point independent quasi particle polarizability, χ0, by the following Dyson equation, χ0 (1, 2; 3, 4) = −iG(1, 3)G(4, 2)

L̃(1, 2; 3, 4) = χ0 (1, 2; 3, 4) −

(5)

∫ d(5, 6)χ0 (1, 2; 5, 6)

W (5, 6)L̃(5, 6; 3, 4)

(6)

where the screened potential W describes the screened interaction between the propagating electron and hole. Because optical excitation energies are provided by the poles of L, one does not need to solve L̃ explicitly and can directly assess L by substituting eq 6 within eq 4. In practice, the poles of L are obtained by studying the kernel of L−1(ω) expressed in the frequency domain, further assuming an adiabatic approximation for W

APPENDIX We will not review here the details of the GW and Bethe− Salpeter formalisms presented in several reviews24 and briefly overview the definition of the calculated Bethe−Salpeter eigenstates and related oscillator strengths only. The BSE formalism generalizes the standard calculation of the reducible polarizability, ∂G(1, 1+ ) ∂ρ(1) = −i ∂Uext(2) ∂Uext(2)

(2)

where Uext is now a nonlocal external perturbation. Similarly to the reducible polarizability χ, L can be described through a Dyson equation involving the bare coulomb potential v and the four point response function L̃ associated with a variation in the total potential

■

χ (1, 2) =

∂G(1, 2) ∂Uext(3, 4)

L(1, 2; 3, 4; ω)−1 = χ0 (1, 2; 3, 4; ω)−1 − δ(1, 2)δ(3, 4)v (1, 3) + δ(1, 3)δ(2, 4)W (1, 2; ω = 0)

(7)

where, e.g., 1 = r1 is now a space coordinate only. In this frequency representation, the (static) screened Coulomb potential can be simply related to the reducible polarizability

(1) J


Article

Journal of Chemical Theory and Computation χ or, alternatively, the inverse (optical) dielectric function in the low-frequency limit

⎛⎡ R − iηI C ⎤ ⎡ I 0 ⎤⎞ ⎟ ⎥ − ω⎢ L−s →1 s(ω) = − ⎜⎜⎢ ⎣ 0 −I ⎥⎦⎟⎠ R − iηI ⎦ ⎝⎣ C

W (1, 2; ω = 0) = v(1, 2)

∬

+

d3 d4v(1, 3)χ (3, 4; ω = 0)v(4, 2)

W (1, 2; ω = 0) =

∫ d3ϵ−1(1, 3; ω = 0)v(3, 2)

⎤ C ⎡−I 0 ⎤ ⎛⎡ R − iηI ⎡ I 0 ⎤⎞ ⎟ ⎥ − ω⎢ =⎢ ·⎜⎜⎢ ⎥ ⎣ 0 I ⎦ ⎝⎣ − C ⎣ 0 I ⎥⎦⎟⎠ −R + iηI ⎦

(8)

This non-Hermitian problem requires taking into account left and right eigenvectors of the right-most member of eq 18 to −1 write the spectral decomposition of Ls→s . However, right and left eigensolutions are linked through the problem symmetry and one can show that only the definition of the right eigenvectors appears in the final expression. Furthermore, the (right) eigensolutions are paired because to any positive + eigenenergy solution (|Φ+R λ ⟩, Eλ − i0 ) correspond a negative + −R energy solution (|Φλ ⟩, −Eλ + i0 ), the two being related by

(9)

As such, in contrast with TD-DFT exchange-correlation kernels associated with (semi)local functionals, it is a nonlocal operator allowing, e.g., to couple nonoverlapping hole and electron distributions of interest in charge-transfer excitations. Noting (ψσa,Eσa) and (ψσi ′,Eσi ′), the GW occupied and virtual molecular orbitals of spin indices σ and σ′, respectively, the expression of χ0 within this representation reads

Φ+λ R (1, 2) =

⎡ ψ σ(1)ψ σ ′*(2)ψ σ *(3)ψ σ ′(4) i a i χ0 (1, 2; 3, 4; ω) = ∑ ⎢ a σ′ σ ⎣ ω − (Ei − Ea ) + iη aiσσ ′ ⎢

⇒

∑ (BλiaΦia(1, 2) + AλiaΦai(1, 2)) (19)

Ls → s(1, 2; 3, 4; ω) ⎛ Φ+R (1, 2)Φ+R (3, 4) Φ−λ R (1, 2)Φ−λ R (3, 4) ⎞ λ ⎟ = ∑⎜ λ − ω − Eλ + iη+ ω + Eλ − iη+ ⎠ λ ⎝ (20)

The reducible polarizability is then obtain through contraction of Ls→s

⎧ ⎫ 1 ⎨Φkl (1, 2) = ψk(1)ψl(2)((α(1)α(2) + β(1)β(2))⎬ ⎩ ⎭ 2

χ (1, 2, ω) = Ls → s(1, 1; 2, 2; ω)

Using the fact that Φλ(1) ≡ a final expression for χ

(11)

χ0 is still diagonal under these conditions and can thus be easily inverted within this representation. In order to facilitate the expression of the singlet diagonal block Ls→s of L, one can introduce the following notations:

∫ d(1, 2)ψm(1)ψk(1)W (1, 2)ψn(2)ψl(2)

2) =

After some algebra work, one finds finally that the spectral representation of Ls→s writes

Assuming symmetric eigenstates in spin, one can show (see Appendix B of ref 53, for example) that L should be first block diagonalized in spin indices to reveal one singlet and three triplet independent subspaces. The basis set supporting the singlet components then reads

≡

Φ−λ R (1,

vc

ψi (1)ψaσ *(2)ψi σ ′*(3)ψaσ(4) ⎤ ⎥ ω − (Eaσ − Eiσ ′) − iη ⎥⎦ (10)

kl W mn

∑ (AλiaΦia(1, 2) + BλiaΦai(1, 2)) ia

σ′

−

χ (1, 2, ω) =

∑ λ

Φ+R λ (1,1)

2Eλ 2

ω − Eλ2

(21)

=

Φ−R λ (1,1),

we obtain as

Φλ(1)Φλ(2) (22)

Defining the dipolar transition moments as dμλ ≡ ⟨μ|Φλ⟩ =

(12)

2

∑ ⟨ψi|μ|ψa⟩(Aλia + Bλia) ia

kl V mn ≡

∫ d(1, 2)ψm(1)ψn(1)v(1, 2)ψk(2)ψl(2)

≡ ⟨Φia|Ls → s(ω = 0)|Φjb⟩

R iajb = δabδij(Ei − Ea) + 2Viajb − Wiajb

(13)

αμν(ω) ≡ ⟨μ|χ (ω)|ν⟩ =

Ciajb = 2Viajb − W iabj

∑ λ

2Eλ 2

ω − Eλ2

dμλdνλ (24)

Defining now the oscillator strength from the averaged polarizability

(14)

α̅ (ω) = (15)

1 Tr[αμν(ω)] = 3

∑ λ

fλ 2

ω − Eλ2

(25)

one obtains the oscillator strength provided in this manuscript 2 fλ = Eλ ∑ (dμλ)2 3 (26) μ

and: Ciajb ≡ ⟨Φia|Ls → s(ω = 0)|Φbj⟩

(23)

we obtain for expression of the polarizability tensor

and then define the resonant (R) and antiresonant (C) parts of the projected L−1 s→s operator as R iajb

(18)

(16)

As a final remark, we note that replacing the screened Coulomb potential W by its bare analogue v in eq 12 leads to the standard time-dependent Hartree−Fock formalism, allowing careful comparisons with standard quantum chemistry codes to validate the present implementation.

(17)

We can now search for the kernel of L−1 s→s as the eigensolution of the matrix pencil K


Article


■

(27) Godby, R. W.; Schlüter, M.; Sham, L. J. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 10159−10175. (28) Tiago, M. L.; Chelikowsky, J. R. Solid State Commun. 2005, 136, 333−337. (29) Tiago, M. L.; Kent, P. R. C.; Hood, R. Q.; Reboredo, F. A. J. Chem. Phys. 2008, 129, 084311. (30) Ma, Y.; Hao, R.; Shao, G.; Wang, Y. J. Phys. Chem. A 2009, 113, 5066−5072. (31) Palummo, M.; Hogan, C.; Sottile, F.; Bagalá, P.; Rubio, A. J. Chem. Phys. 2009, 131, 084102. (32) Kaczmarski, M. S.; Ma, Y.; Rohlfing, M. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 115433. (33) Ma, Y.; Rohlfing, M.; Molteni, C. J. Chem. Theory Comput. 2010, 6, 257−265. (34) Rocca, D.; Lu, D.; Galli, G. J. Chem. Phys. 2010, 133, 164109. (35) Garcia-Lastra, J. M.; Thygesen, K. S. Phys. Rev. Lett. 2011, 106, 187402. (36) Blase, X.; Attaccalite, C. Appl. Phys. Lett. 2011, 99, 171909. (37) Duchemin, I.; Deutsch, T.; Blase, X. Phys. Rev. Lett. 2012, 109, 167801. (38) Baumeier, B.; Andrienko, D.; Ma, Y.; Rohlfing, M. J. Chem. Theory Comput. 2012, 8, 997−1002. (39) Faber, C.; Duchemin, I.; Deutsch, T.; Blase, X. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 155315. (40) Hogan, C.; Palummo, M.; Gierschner, J.; Rubio, A. J. Chem. Phys. 2013, 138, 024312. (41) Faber, C.; Boulanger, P.; Duchemin, I.; Attaccalite, C.; Blase, X. J. Chem. Phys. 2013, 139, 194308. (42) Duchemin, I.; Blase, X. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 245412. (43) Varsano, D.; Coccia, E.; Pulci, O.; Mosca Conte, A.; Guidoni, L. Comput. Theor. Chem. 2014, 1040−1041, 338−346. (44) Coccia, E.; Varsano, D.; Guidoni, L. J. Chem. Theory Comput. 2014, 10, 501−506. (45) Baumeier, B.; Rohlfing, M.; Andrienko, D. J. Chem. Theory Comput. 2014, 10, 3104−3110. (46) Boulanger, P.; Jacquemin, D.; Duchemin, I.; Blase, X. J. Chem. Theory Comput. 2014, 10, 1212−1218. (47) Boulanger, P.; Chibani, S.; Le Guennic, B.; Duchemin, I.; Blase, X.; Jacquemin, D. J. Chem. Theory Comput. 2014, 10, 4548−4556. (48) Koval, P.; Foerster, D.; Sánchez-Portal, D. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 155417. (49) Körbel, S.; Boulanger, P.; Duchemin, I.; Blase, X.; Marques, M. A. L.; Botti, S. J. Chem. Theory Comput. 2014, 10, 3934−3943. (50) Hahn, T.; Geiger, J.; Blase, X.; Duchemin, I.; Niedzialek, D.; Tscheuschner, S.; Beljonne, D.; Bässler, H.; Köhler, A. Adv. Funct. Mater. 2015, 25, 1287−1295. (51) Krause, K.; Harding, M. E.; Klopper, W. Mol. Phys. 2015, 113, 1952−1960. (52) Hirose, D.; Noguchi, Y.; Sugino, O. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 205111. (53) Ljungberg, M. P.; Koval, P.; Ferrari, F.; Foerster, D.; SanchezPortal, D. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 92, 075422. (54) Bruneval, F.; Hamed, S. M.; Neaton, J. B. J. Chem. Phys. 2015, 142, 244101. (55) Jacquemin, D.; Duchemin, I.; Blase, X. J. Chem. Theory Comput. 2015, 11, 3290−3304. (56) Jacquemin, D.; Duchemin, I.; Blase, X. J. Chem. Theory Comput. 2015, 11, 5340−5359. (57) Jacquemin, D.; Duchemin, I.; Blase, X. Mol. Phys. 2016, 114, 957−967. (58) Schreiber, M.; Silva-Junior, M. R.; Sauer, S. P. A.; Thiel, W. J. Chem. Phys. 2008, 128, 134110. (59) Silva-Junior, M. R.; Schreiber, M.; Sauer, S. P. A.; Thiel, W. J. Chem. Phys. 2008, 129, 104103. (60) Rebolini, E.; Toulouse, J. J. Chem. Phys. 2016, 144, 094107. (61) Tawada, T.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. J. Chem. Phys. 2004, 120, 8425−8433.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS D.J. acknowledges K. Burke (UCI) for suggesting to go beyond the analysis of transition energies with BSE/evGW. The authors are indebted to Fabien Bruneval for benchmark comparisons with the MOLGW code. D.J. acknowledges the European Research Council (ERC) and the Région des Pays de la Loire for financial support in the framework of a Starting Grant (Marches -278845) and the LumoMat project, respectively. X.B. acknowledges partial support from the French National Research Agency under contract ANR-12-BS04 “PANELS”. This research used resources of (i) the GENCI-CINES/IDRIS, (ii) CCIPL (Centre de Calcul Intensif des Pays de Loire), (iii) a local Troy cluster, and (iv) HPC resources from ArronaxPlus (grant ANR-11-EQPX-0004 funded by the French National Agency for Research).

■

REFERENCES

(1) Dreuw, A.; Head-Gordon, M. Chem. Rev. 2005, 105, 4009−4037. (2) Nakatsuji, H. J. Chem. Phys. 1991, 94, 6716−6727. (3) Nakatsuji, H.; Ehara, M. J. Chem. Phys. 1993, 98, 7179−7184. (4) Bak, K. L.; Koch, H.; Oddershede, J.; Christiansen, O.; Sauer, S. P. A. J. Chem. Phys. 2000, 112, 4173−4185. (5) Falden, H. H.; Falster-Hansen, K. R.; Bak, K. L.; Rettrup, S.; Sauer, S. P. A. J. Phys. Chem. A 2009, 113, 11995−12012. (6) Koch, H.; Jorgensen, P. J. Chem. Phys. 1990, 93, 3333−3344. (7) Stanton, J. F.; Bartlett, R. J. J. Chem. Phys. 1993, 98, 7029−7039. (8) Christiansen, O.; Koch, H.; Jørgensen, P. Chem. Phys. Lett. 1995, 243, 409−418. (9) Christiansen, O.; Koch, H.; Jorgensen, P. J. Chem. Phys. 1995, 103, 7429−7441. (10) Koch, H.; Christiansen, O.; Jorgensen, P.; Sanchez de Merás, A. M.; Helgaker, T. J. Chem. Phys. 1997, 106, 1808−1818. (11) Kállay, M.; Gauss, J. J. Chem. Phys. 2004, 121, 9257−9269. (12) Caricato, M. J. Chem. Phys. 2013, 139, 114103. (13) Schirmer, J.; Trofimov, A. B. J. Chem. Phys. 2004, 120, 11449− 11464. (14) Hellweg, A.; Grün, S. A.; Hättig, C. Phys. Chem. Chem. Phys. 2008, 10, 4119−4127. (15) Krauter, C. M.; Pernpointner, M.; Dreuw, A. J. Chem. Phys. 2013, 138, 044107. (16) Dreuw, A.; Wormit, M. WIREs Comput. Mol. Sci. 2015, 5, 82− 95. (17) Runge, E.; Gross, E. K. U. Phys. Rev. Lett. 1984, 52, 997−1000. (18) Casida, M. E. In Time-Dependent Density-Functional Response Theory for Molecules; Chong, D. P., Ed.; Recent Advances in Density Functional Methods; World Scientific: Singapore, 1995; Vol. 1, pp 155−192. (19) Laurent, A. D.; Adamo, C.; Jacquemin, D. Phys. Chem. Chem. Phys. 2014, 16, 14334−14356. (20) Laurent, A. D.; Jacquemin, D. Int. J. Quantum Chem. 2013, 113, 2019−2039. (21) Sham, L. J.; Rice, T. M. Phys. Rev. 1966, 144, 708−714. (22) Hanke, W.; Sham, L. J. Phys. Rev. Lett. 1979, 43, 387−390. (23) Strinati, G. Phys. Rev. Lett. 1982, 49, 1519−1522. (24) Onida, G.; Reining, L.; Rubio, A. Rev. Mod. Phys. 2002, 74, 601− 659. (25) Hedin, L. Phys. Rev. 1965, 139, 796−823. (26) Hybertsen, M. S.; Louie, S. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 34, 5390−5413. L


Article

Journal of Chemical Theory and Computation (62) Miura, M.; Aoki, Y.; Champagne, B. J. Chem. Phys. 2007, 127, 084103. (63) Szalay, P. G.; Watson, T.; Perera, A.; Lotrich, V. F.; Bartlett, R. J. J. Phys. Chem. A 2012, 116, 6702−6710. (64) Caricato, M.; Trucks, G. W.; Frisch, M. J.; Wiberg, K. B. J. Chem. Theory Comput. 2011, 7, 456−466. (65) Silva-Junior, M. R.; Sauer, S. P. A.; Schreiber, M.; Thiel, W. Mol. Phys. 2010, 108, 453−465. (66) Sauer, S. P.; Pitzner-Frydendahl, H. F.; Buse, M.; Jensen, H. J. A.; Thiel, W. Mol. Phys. 2015, 113, 2026−2045. (67) Kánnár, D.; Szalay, P. G. J. Chem. Theory Comput. 2014, 10, 3757−3765. (68) Harbach, P. H. P.; Wormit, M.; Dreuw, A. J. Chem. Phys. 2014, 141, 064113. (69) Labhart, H. Helv. Chim. Acta 1957, 40, 1410−1421. (70) Casida, M. E. J. Mol. Struct.: THEOCHEM 2009, 914, 3−18. (71) Ullrich, C. Time-Dependent Density-Functional Theory: Concepts and Applications; Oxford Graduate Texts; Oxford University Press: New York, 2012. (72) Blase, X.; Attaccalite, C.; Olevano, V. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 115103. (73) Faber, C.; Attaccalite, C.; Olevano, V.; Runge, E.; Blase, X. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 115123. (74) Kaplan, F.; Harding, M. E.; Seiler, C.; Weigend, F.; Evers, F.; van Setten, M. J. J. Chem. Theory Comput. 2016, 12, 2528−2541. PMID: 27168352. (75) Weigend, F.; Köhn, A.; Hättig, C. J. Chem. Phys. 2002, 116, 3175−3183. (76) Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; Van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; de Jong, W. A. Comput. Phys. Commun. 2010, 181, 1477−1489. (77) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215−241. (78) Silva-Junior, M. R.; Schreiber, M.; Sauer, S. P. A.; Thiel, W. J. Chem. Phys. 2010, 133, 174318. (79) Isegawa, M.; Peverati, R.; Truhlar, D. G. J. Chem. Phys. 2012, 137, 244104. (80) Baerends, E. J.; Gritsenko, O. V.; van Meer, R. Phys. Chem. Chem. Phys. 2013, 15, 16408−16425. (81) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision D.01; Gaussian Inc.: Wallingford, CT, 2009. (82) TURBOMOLE, V6.6 2014, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH, since 2007; available from http://www. turbomole.com (accessed June 13, 2016). (83) Perpète, E. A.; Wathelet, V.; Preat, J.; Lambert, C.; Jacquemin, D. J. Chem. Theory Comput. 2006, 2, 434−440. (84) These six states are 3 B3u in naphthalene, 4 A″ and 5 A′ of imidazole, 2 E′ of s-triazine as well as 7 A′ and 8 A′ of formamide (using the numbering conventions of ref 67). (85) Jacquemin, D.; Brémond, E.; Planchat, A.; Ciofini, I.; Adamo, C. J. Chem. Theory Comput. 2011, 7, 1882−1892.

M


GW

Recommend Documents