The Bethe-Salpeter equation (BSE) based on GW quasiparticle levels is a successful approach for calculating the optical gaps and spectra of solids and also for predicting the neutral excitations of small molecules. We here present an all-electron implementation of the GW+BSE formalism for molecules, using numeric atom-centered orbital (NAO) basis sets. We present benchmarks for low-lying excitation energies for a set of small organic molecules, denoted in the literature as “Thiel’s set.” Literature reference data based on Gaussian-type orbitals are reproduced to about one millielectron-volt precision for the molecular benchmark set, when using the same GW quasiparticle energies and basis sets as the input to the BSE calculations. For valence correlation consistent NAO basis sets, as well as for standard NAO basis sets for ground state density-functional theory with extended augmentation functions, we demonstrate excellent convergence of the predicted low-lying excitations to the complete basis set limit. A simple and affordable augmented NAO basis set denoted “tier2+aug2” is recommended as a particularly efficient formulation for production calculations. We finally demonstrate that the same convergence properties also apply to linear-response time-dependent density functional theory within the NAO formalism.

Predicting the neutral (including optical) excitations of molecules and materials is of fundamental importance in photovoltaics, optoelectronics, and other technologically relevant areas. Several distinct types of computational formalisms are frequently employed in the community for this purpose, including wavefunction-based methods, e.g., equation-of-motion coupled cluster (EOM-CC),1–3 complete active space second-order perturbation theory (CASPT2),4–7 the quantum Monte Carlo method,8–11 linear-response time-dependent density functional theory (LR-TDDFT),12–15 or the Bethe-Salpeter equation (BSE) approach in the context of many-body perturbation theory (MBPT).16–23 EOM-CC and CASPT2 have been shown to produce highly accurate values for small and mid-sized molecules, when combined with sufficiently high-quality basis sets. They are therefore often used as a trusted reference,1–3,5–7,24 although their applicability to larger systems is somewhat limited by the associated computational cost. LR-TDDFT has been widely applied to predict optical excitations for molecules due to its computational efficiency and often reasonable accuracy, especially when combined with carefully designed exchange-correlation (XC) functionals.25–30 However, LR-TDDFT calculations can encounter problems for charge-transfer (CT) excitations,31,32 especially when used with a simple XC functional such as the adiabatic local density approximation (LDA)33 and generalized gradient approximations (GGAs).34 In LR-TDDFT, including long-range exact exchange in the XC functional can mitigate this problem.28–30,35

The BSE approach is founded upon MBPT, based on Green’s function (G) theory and the idea of using the screened Coulomb interaction W.16 The BSE formalism was originally proposed in the field of nuclear physics in the 1950s.17 Combined with the GW approximation in MBPT,16 the BSE approach has been shown to successfully approximate the optical spectra of solids19,22,36–38 and nanosystems39 and later work demonstrates similar applicability to excitations in atoms and molecules.20,40–51 The GW approach16,21,52,53 allows one to predict fundamental gaps—i.e., gaps between highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) quasiparticle states—as well as single-quasiparticle excitation spectra that are more accurate than those obtained by standard density functional theory (DFT) for a wide range of systems, including both solids and molecules.54–68 The description of optical excitations within the BSE approach then uses charged excitations, i.e., electron removal and addition excitations, from the GW approach as its input. The BSE method based on the GW method has several formal advantages over LR-TDDFT. The electron-hole interaction in the BSE approach has the correct asymptotic behavior in the so-called “Mulliken limit”66 for both solids and molecules [i.e., Ref. 66: “ECT0(D) = EA(acceptor) − IP(donor) − 1/D, where EA(acceptor) and IP(donor) are the acceptor electronic affinity and the donor ionization energy” for large acceptor-donor distance D, all in atomic units]. This limit is not captured by the LR-TDDFT formalism without a long range exchange component.43,66 CT excitations, which are problematic for LR-TDDFT especially with LDA and GGA functionals, can be efficiently and accurately predicted by the BSE approach.48,66,69–71 This has been demonstrated for CT excitations including both intermolecular and intramolecular types in systems such as simple dipeptides,70,72 and more complex fullerene/polymer aggregates.71,73

Calculations of BSE excitation energies within MBPT usually adopt a three-step procedure: (i) Evaluate the Kohn-Sham (KS)33 or generalized Kohn-Sham (gKS)74 DFT orbitals. (ii) Apply self-energy corrections at the G0W0 level or partially self-consistent23,51GW level (G0W0@DFT or GW@DFT; G0 stands for Green’s function of a noninteracting reference system and W0 is the screened Coulomb interaction of that reference system).16,20,21,52 (iii) Solve the BSE (in practice, an approximate version thereof, see below) based on the G0W0 or GW quasiparticle energies and screened and unscreened Coulomb integrals of (g)KS orbitals (BSE@G0W0@DFT or BSE@GW@DFT).17,43,44,46 BSE implementations exist in different computational packages based on different basis functions, e.g., MolGW,75 Fiesta,69,76 and Turbomole,46,77,78 which are based on Gaussian-type orbital (GTO) basis sets, BerkeleyGW,79 Yambo,80 ABINIT,81,82 VASP,83 and Quantum Espresso,84 which are based on plane waves, or Exciting85,86 and Elk,87 which are based on augmented plane waves.

The present work introduces an accurate implementation of the BSE formalism utilizing compact and efficient numeric atom-centered orbital (NAO) basis sets88,89 in the context of the all-electron electronic structure code FHI-aims.68,88,90,91 This paper focuses on spin-unpolarized systems (spin-polarized calculations are not yet possible in the present implementation). To obtain the two-electron Coulomb and screened Coulomb interaction matrix elements, we use an efficient and highly accurate variant of the resolution-of-identity (RI) technique.68 In FHI-aims, this RI technique is the numerical foundation for all methods beyond semilocal DFT, including Hartree-Fock (HF), hybrid density functionals, the random-phase approximation (RPA), second-order Møller-Plesset perturbation theory (MP2), and the GW method.68,91 Our current implementation also uses the ELSI infrastructure92 and the ELPA eigenvalue solver93 for parallel eigenvalue solutions.

This paper is organized as follows: In Sec. II, we introduce the GW+BSE formalism in the context of MBPT. In Sec. III, we discuss the details of our implementation. In Sec. IV, we demonstrate the numerical correctness of our BSE implementation by comparing the excitation energies computed by FHI-aims and by the MolGW code43,75 with GTO correlation-consistent basis sets for Thiel’s molecular benchmark set.24,26 In assessing our BSE implementation, we emphasize the dependence of the BSE results on the GW quasiparticle energies. We then study the convergence behavior of excitation energies to the complete basis set limit, combining standard NAO basis sets for ground state DFT (FHI-aims-2009)88 or valence correlation consistent NAO basis sets (NAO-VCC-nZ)89 with extended augmentation functions (NAO+aug). We demonstrate that the standard FHI-aims-2009 basis sets give essentially basis set converged numerical results for low-lying optical excitation energies when combined with a few extended augmentation functions (NAO+aug basis sets) that are also commonly included in Gaussian-type basis sets. Finally, similar convergence behavior is demonstrated for LR-TDDFT with adiabatic LDA as the kernel.33,94

Typical calculations of the neutral (optical) excitation energies of molecules using the BSE approach adopt the following three-step procedure, which is utilized by a wide range of electronic structure packages for calculations of neutral (optical) excitation properties in the framework of MBPT:43,44,46,79,80,82,83,85

(i) The initial step is performed by solving the self-consistent (g)KS equations with an approximate functional for the exchange-correlation energy Exc. Common choices for Exc are the LDA, GGAs, HF, and hybrid functionals. In KS theory (e.g., LDA and GGAs), we define v^xc as the functional derivative of Exc with respect to the electron density. In the gKS case (HF and hybrid functionals), v^xc is the functional derivative with respect to the set of orbitals ψl. In either case, the ψl are constructed as

ĥ(g)KSψl=εlψl,
(1)
ĥ(g)KS=t^s+v^ext+v^H+v^xc.
(2)

Equation (1) states the electronic (g)KS single-particle equations for the effective single-particle orbitals ψl and eigenvalues εl (l = 1, 2, …, Norbit). Equation (2) details the gKS Hamiltonian ĥ(g)KS, including the effective single-particle kinetic energy (with relativistic corrections) t^s, the external potential v^ext, the electrostatic or Hartree potential of the electron density v^H, and the exchange-correlation potential v^xc. The underlying orbitals {ψl} are here expanded on the basis of NAOs {φi, i = 1, 2, …, Nbasis} as

ψl=icliφi,
(3)

where the NAOs are of the form88 

φi(r)=ui(r)rYlm(Ω),
(4)

where r is a position vector with respect to the nucleus, r is its modulus, and Ω is the corresponding angular coordinate on the unit sphere. In the FHI-aims code, the ui are numerically tabulated functions, defined as cubic splines in units of a logarithmic grid. Ylm are the real-valued spherical harmonics, and l and m are implicitly included in the basis function index i. The eigenvalues and eigenfunctions produced by the initial step serve as a first guess for the quasiparticles and are used to evaluate the Coulomb interaction, the screened Coulomb interaction, and the GW self-energy in the subsequent GW and BSE@GW steps. Although, in the nonperiodic case, ψl can be chosen to be real-valued, we include complex conjugates in the derivations below.

(ii) A perturbative GW approach is then applied to obtain the quasiparticle energies as follows:68 

εlGW=εl+ψl|ΣGW(εlGW)v^xc|ψl,
(5)

where εlGW is the quasiparticle energy. By convention, the arguments εlGW used to evaluate the self-energy ΣGW on the right-hand side are updated self-consistently until they match the εlGW values obtained on the left-hand side, even though the function ΣGW(ω) itself is not further updated in the process. The GW self-energy is calculated from the Green’s function G and the screened Coulomb potential W following the GW approximation proposed by Hedin:16 

ΣGWr,r,ω=i2πdωGr,r,ω+ωWr,r,ωeiωη.
(6)

In the single-shot perturbative GW (i.e., G0W0) approach, the Green’s function G is approximated by the noninteracting Green’s function G0, which is calculated from single-particle orbitals ψl and orbital energies εl,68 

G0(r,r,ω)=lψl(r)ψl(r)ωεliηsgn(εFεl).
(7)

where εF is the Fermi energy and η is a positive infinitesimal. The screened Coulomb potential W0 is calculated from the dielectric function ε as68 

W0(r,r,ω)=drε1(r,r,ω)v(r,r),
(8)

where the dielectric function ε is obtained at the RPA level, using DFT results. The G0W0 self-energy can be calculated using an exact analytic treatment on the real axis, which is the case in the MolGW package.60,75 We refer to Sec. III.C of Ref. 60 for the details of this formalism. This treatment is limited to small systems. Instead, two-pole95 and Padé96 approximations are implemented in the FHI-aims code for the evaluation of the self-energy on the real axis.68 Both of these approximations are based on an exact treatment of G0, W0, and the self-energy ΣG0W0 on the imaginary frequency axis. ΣG0W0 is then extended to the real axis by performing an analytic fit of the data on the imaginary axis to a function with a form that has poles on the real axis. This process is usually referred to as “analytic continuation.” The smooth behavior of all quantities (G0, W0, ΣG0W0) on the imaginary frequency axis significantly reduces the number of frequency points needed, compared to a full frequency integration along the real frequency axis.97 Specifically, the self-energy is approximated to have the following mathematical form in the complex plane in the two-pole approximation:95 

Σij(z)n=12anz+bn,
(9)

where the values of an and bn depend on the indices i and j. In the Padé approximation, the self-energy is expressed as96 

Σij(z)a0+a1z++a(N2)/2(z)(N2)/21+b1z++bN/2(z)N/2,
(10)

where N denotes the total number of parameters in the Padé approximation. We note already here that the Padé approximation can be more accurate than the two-pole approximation to represent the true self-energy but that the Padé approximation is also, in practice, more prone to numerical problems, including nonunique solutions that can be difficult to control without manual inspection of all resulting eigenvalues. In addition to the two approaches mentioned above, another, more elaborate approach to evaluate the self-energy directly on the real axis by contour deformation (CD) was implemented in FHI-aims by Golze and co-workers98 while this paper was being completed. We do not assess this approach here because our emphasis here is on the BSE, but we note that essentially exact G0W0 input data to the BSE are expected from the CD approach. On the other hand, the analytical continuation of Σ according to Eq. (9) or Eq. (10) is advantageous over the CD approach in terms of computational cost, both in terms of the base cost (often called the prefactor) and in terms of the scaling exponent with system size if the number of needed G0W0 eigenvalues scales with the size of the system.98 

Here, we perform one-shot perturbative G0W0 calculations based on a fixed DFT or HF reference. The quasiparticle energy in Eq. (5) is thus rewritten as εlG0W0. Some studies investigate the effect of iterating the GW equations by updating the eigenvalues in Eq. (7) by those calculated from Eq. (5), whereas the wavefunctions ψl are kept at the DFT level.44,49,99 This procedure, denoted eigenvalue self-consistent GW, is reported to give better agreement with experimental results and wavefunction-based reference methods compared with single-shot G0W0 for some systems.44,45,49,99

(iii) The BSE is a Dyson-like equation for the two-particle correlation function L,20,21,43

L(12;12)=L0(12;12)+d(3456)L0(14;13)×K(35;46)L(62;52),
(11)

where the set of variables 1, 2, etc. is short for position, time, and spin (r1, t1, σ1), (r2, t2, σ2), etc. L(12; 1′2′) is the electron-hole correlation function which describes the probability amplitude of an electron propagating from 1′ to 2 and a hole propagating from 1 to 2′.21,L0(12; 1′2′) represents the correlation function of the noninteracting system defined below in Eq. (12). K(35; 46), usually called the electron-hole interaction kernel, is the screened interaction between the electron and the hole (including bare exchange). L0 and K can be expressed in the following equations:20 

L0(12;12)=G0(1,2)G0(2,1),
(12)
K(35;46)=δM(3,4)δG(6,5),
(13)

where G is the one-particle Green’s function and M is equal to the sum of the self-energy and the Hartree potential,

M(3,4)=vH(3)δ(3,4)+Σ(3,4).
(14)

By applying Eq. (14) to Eq. (13), performing a time-energy Fourier transformation and ignoring the dynamical properties of W,20 the BSE kernel can be simplified to

K(r3,r5,r4,r6)=iv(r3r5)δ(r3r4)δ(r5r6)+iW(r3,r4,ω=0)δ(r3r6)δ(r4r5),
(15)

where the variables 3, 4, etc. are reduced to r3, r4, etc., v is the bare Coulomb interaction, and W is the screened Coulomb interaction, with the frequency-dependence ignored.20 This approximation means that the actual BSE part (once the GW quasiparticle energies are fixed) is independent of a particular analytical continuation choice since only the ω = 0 value of W enters the approximated BSE.

In practical implementations, the BSE is usually rewritten in the following matrix form in a transition space spanned by the products of occupied and unoccupied orbitals:20,43,44

ABBAXsYs=EsXsYs,
(16)

where Es are the optical excitation eigenvalues and (Xs, Ys) are the eigenvectors. The Xs and Ys are expressed in electron-hole space of the unperturbed system with elements Xs,ia and Ys,ia, i.e., the actual BSE wavefunctions are linear combinations of the product of (g)KS orbitals. The excitation wavefunctions Xs, Ys can be taken to be real-valued in finite (molecular) systems without an external field. Blocks A and −A correspond to resonant transitions from occupied to unoccupied orbitals and the antiresonant transitions, respectively.100 Blocks B and −B describe the coupling between blocks A and −A. In the BSE, the diagonal matrices A(−A) and off-diagonal matrices B(−B) are defined as follows:20,43,44

Aiajb=(εaGWεiGW)δijδabαS/T(ia|V|jb)+(ij|W(ω=0)|ab),
(17)
Biajb=αS/T(ia|V|bj)+(ib|W(ω=0)|aj).
(18)

The indices i and j denote occupied states and a and b denote unoccupied states. εaGW and εiGW are the quasiparticle energies denoted εlGW in Eq. (5). The coefficient αS/T is equal to 2 for singlet excitations and 0 for triplet excitations. The index conventions for the bare and screened Coulomb interactions V and W are as follows:75 

(ia|V|jb)=pqrscipcaqcjrcbs(pq|rs),
(19)
(ij|W|ab)=pqrscipcjqcarcbs(pq|W|rs),
(20)

where (pq|rs) are the two-electron integrals in a basis set representation,

(pq|rs)=φp(r)φq(r)φr(r)φs(r)|rr|drdr,
(21)

and the same convention (Mulliken notation) for r and r′ is used in the notation of the screened Coulomb integrals (pq|W|rs) as well. The neglect of the coupling blocks B(−B) in Eq. (16) is known as the Tamm-Dancoff approximation (TDA).37,101 In the TDA, which also we compare below, the relevant equation becomes simply,

AXs=EsXs.
(22)

The oscillator strength fs can be calculated from the eigenvalues and eigenvectors obtained by solving the BSE eigenvalue problem as follows:102 

fs=23Esμ=x,y,z(ds,μ)2,
(23)

where ds,μ can be calculated as102 

ds,μ=2iaψi|μ^|ψa(Xs,ia+Ys,ia).
(24)

Since we are dealing with finite systems, the dipole operator μ^ is simply taken to be the position operator, i.e., μ^(x,y,z). For convenience, we reference the coordinates x, y, z to the center (average of atomic positions) of the molecule.

We will also compare our observations for the BSE to analogous results for LR-TDDFT, which is widely used in chemistry. We therefore briefly recapitulate the LR-TDDFT formalism, the mathematical structure of which is similar to the BSE, albeit with a two-point kernel instead of the four-point kernel of BSE. A deeper discussion of the mathematical similarities and differences of both levels of theory is given in Ref. 22. LR-TDDFT is often expressed as the Casida eigenvalue equation,103 which is formally equivalent to Eq. (16). Here, the LR-TDDFT formalism becomes

ΩFs=Es2Fs.
(25)

In LR-TDDFT, Ω is called the Casida matrix, which has the same dimension as A or B in Eq. (16). Es2 are squares of the neutral many-body excitation energies and Fs are the eigenvectors of this eigenvalue problem, which can also be related to the oscillator strengths via the dipole operator.104 The Casida matrix can be written in the basis of products of (g)KS orbitals as

Ωia,jb(ω)=δi,jδa,b(εaεi)2+2(εaεi)Kia,jb(ω)(εbεj),
(26)

where δ denotes the Kronecker delta. The kernel Kia,jb is defined as

Kia,jb(ω)=ψi(r)ψa(r)1|rr|+fxc[n0](r,r,ω)ψj(r)ψb(r)drdr,
(27)

where fxc[n0] is the exchange correlation kernel and n0 is the (g)KS ground state electron density. As in the case of the BSE approach, the LR-TDDFT method can also employ the TDA approximation as discussed in Ref. 105. Then, an analogous equation to Eq. (22) is solved but with the matrix A being

Aia,jb(ω)=δi,jδa,b(ϵaϵi)+Kia,jb(ω).
(28)

The kernel is formally defined through the functional derivative of the time-dependent Kohn-Sham exchange and correlation potential vxc[n](r, t) with respect to the time-dependent density n(r′, t′) such that

fxc[n0](r,r,ω)=d(tt)eiω(tt)δvxc[n](r,t)δn(r,t)|n0.
(29)

In practice, the so-called adiabatic approximation106 is employed in Eq. (29), as we do here, and the exchange-correlation kernel reads

fxcA[n0](r,r)=δvxc[n](r)δn(r)|n0.
(30)

This approximation makes the exchange-correlation kernel frequency-independent.

In our implementation, the two-electron Coulomb interaction in Eqs. (17), (18), and (26), the static screened Coulomb interaction in Eqs. (17) and (18), and the two-electron integrals of the exchange correlation kernel in Eq. (26) are calculated employing the RI approach.68,107–111 The RI represents pair products of atomic basis functions φp(r) ⋅ φq(r) in terms of auxiliary basis functions (ABFs),

φp(r)φq(r)μCpqμPμ(r),
(31)

where Pμ(r) (μ = 1, 2, …, Naux) are the ABFs and Cpqμ are the expansion coefficients. The construction of the ABFs in FHI-aims is explained in Ref. 68 and in detail in Ref. 91. The evaluation of the integrals (21) then reduces to

(pq|rs)μνCpqμ(μ|ν)Crsν,with
(32)
(μ|ν)=Pμ(r)Pν(r)|rr|drdr.
(33)

The computation of the expansion coefficients Cpqμ requires three-center integrals involving the ABFs and the pair products of the NAOs,

Cpqμ=ν(pq|ν)(ν|μ)1,where
(34)
(pq|ν)=φp(r)φq(r)Pν(r)|rr|drdr,
(35)

and (ν|μ)−1 denotes the inverse of the Coulomb matrix in ABF representation. Thus, the expensive computation of four-center integrals (pq|rs) is reduced to the computation of much cheaper three-center and two-center ones,

(pq|rs)μOpqμOrsμ,
(36)

using

Opqμ=νCpqν(ν|μ)1/2,
(37)

where (ν|μ)−1/2 denotes the square root of the inverse Coulomb matrix. This enables the efficient computation of the Coulomb matrix elements both in time and in memory. The screened Coulomb interaction W can be represented in terms of the ABFs in a similar way to the Coulomb interaction V as follows:

(pq|W|rs)=μνOpqμ(ε1)μνOrsν.
(38)

The dielectric function ε can be calculated as

εμν=δμναβ(μ|α)1/2χαβ0(β|ν)1/2,
(39)

where χ0 is the noninteracting density response function. In real space and for a non-spin-polarized system, χ0 is defined as

χ0(r,r,iω)=2iaψi(r)ψa(r)ψa(r)ψi(r)iω+εaεi+c.c.,
(40)

where “c.c.” denotes the complex conjugate. We refer to Ref. 68 for more details. The current BSE implementation in FHI-aims uses global RI.68 In the future, use of the localized RI formalism91 is expected to facilitate scalability to larger systems as well as support for extended (periodic) systems.

We quantify the precision of our BSE implementation by calculating the neutral excitation energies of the molecular benchmark set published by Schreiber et al.,24 known in the literature as “Thiel’s set.” This set (see Fig. 1 in Ref. 24) includes N = 28 small and medium-sized organic molecules, the largest of which is naphthalene with 18 atoms. The chemical elements represented in these organic molecules are H, C, N, and O. The atomic coordinates for the included molecules are taken from the supplementary material of Ref. 24. Schreiber et al. focused on obtaining “Best Estimate (BE)” values for singlet and triplet excitation energies of these molecules from ab initio calculations, including rather demanding multireference, coupled cluster, or complete active space approaches of their own or from the literature. While the BE values have been used as reference data by both Thiel’s and other groups for implementations, evaluation and development of a variety of methods,26,43,44 our present study is aimed at establishing the numerical precision of our approach at a fixed level of theory, i.e., BSE or LR-TDDFT. The primary objective of this paper is therefore not to compare to the BE results but rather to compare the BSE excitation energies calculated by our present NAO-based implementation to results obtained at the same level of theory, using the MolGW code as a benchmark. Regarding the basis set, Schreiber et al.24 (and therefore, also some of the results from other methods available for comparison in the literature) employed a relatively limited basis set level for correlated calculations, the TZVP basis set by Schäfer et al.112 A previous study by Bruneval et al.43 indicates that BSE@G0W0@DFT-B3LYP113 excitation energies for ethene and pyrrole, obtained with the TZVP basis set, overestimate the analogous results from the much larger aug-cc-pVQZ basis set114 by about 0.45–0.65 eV. Therefore, the goal of our following investigation is twofold. We first establish the numerical precision of our own implementation in comparison with MolGW using the TZVP basis set.112 We then discuss basis set convergence for low-lying singlet and triplet excitations using NAO basis sets.

In Fig. 1, we compare the numerical precision of the present BSE implementation and that in MolGW using the TZVP basis set. Specifically, we show state-resolved mean absolute error (MAE) values of the BSE-approximated energies of the lowest ten singlet and triplet excited states, respectively, of all molecules in Thiel’s set. The state-resolved MAE, MAE(i), of a given dataset D = {Di,n} in comparison with a reference set R = {Ri,n} is defined as

MAE(i)=1Nn=1N|Di,nRi,n|,
(41)

where i = 1, …, 10 is the index for the state and n = 1, …, N is the index for the molecules in Thiel’s set. For the MAE plotted in Fig. 1, the dataset D is the set of BSE excitation energies calculated using the FHI-aims implementation. The reference set R is the set of BSE excitation energies calculated by MolGW. The Perdew-Burke-Ernzerhof (PBE)34 exchange-correlation functional is used for the initial DFT calculations, and the G0W0 quasiparticle energies that enter the BSE are taken from MolGW in both sets. Panels (a) and (b) show the results for singlet excitations with and without TDA, respectively. Panels (c) and (d) show the results for triplet excitations with and without TDA, respectively. The BSE results from the present implementation agree with the results of the MolGW package at the level of 1 meV or below.

FIG. 1.

Mean absolute error MAE(i) [Eq. (41)] of the BSE@G0W0@PBE lowest 10 singlet/triplet excitation energies from our implementation, compared with reference values from MolGW. The GTO-type TZVP basis set112 is used. The G0W0 quasiparticle energies for the FHI-aims BSE calculations are here taken from MolGW, ensuring that the comparison is specific to the BSE part of the calculations. Panels (a) and (b) show the MAE of the singlet states with and without the Tamm-Dancoff approximation (TDA), respectively. Panels (c) and (d) show the MAE of the triplet states with and without TDA, respectively.

FIG. 1.

Mean absolute error MAE(i) [Eq. (41)] of the BSE@G0W0@PBE lowest 10 singlet/triplet excitation energies from our implementation, compared with reference values from MolGW. The GTO-type TZVP basis set112 is used. The G0W0 quasiparticle energies for the FHI-aims BSE calculations are here taken from MolGW, ensuring that the comparison is specific to the BSE part of the calculations. Panels (a) and (b) show the MAE of the singlet states with and without the Tamm-Dancoff approximation (TDA), respectively. Panels (c) and (d) show the MAE of the triplet states with and without TDA, respectively.

Close modal

Next, Fig. 2 compares BSE oscillator strengths, Eqs. (23) and (24), from both implementations for singlet states in the TDA. The MAEs for the oscillator strengths are at the level of 10−4 or below, i.e., numerically negligible. The actual value of the excitation energy and oscillator strength investigated in this section for all the molecules in Thiel’s set can be found in Tables S1–S5 of the supplementary material. Since the calculation results of the FHI-aims and MolGW package agree within 1 meV for the excitation energy and 10−4 for the oscillator strength, the values in Tables S1–S5 are valid for both packages within the significant digits given. In short, Figs. 1 and 2 validate our implementation.

FIG. 2.

Mean absolute error MAE(i) of the BSE@G0W0@PBE oscillator strengths from the present implementation compared to that of MolGW, using the TDA and the GTO-type TZVP basis set112 for the lowest 10 singlet excitation states, averaged over all molecules in Thiel’s set. As in Fig. 1, the G0W0 quasiparticle energies for the FHI-aims BSE calculations are here taken from MolGW.

FIG. 2.

Mean absolute error MAE(i) of the BSE@G0W0@PBE oscillator strengths from the present implementation compared to that of MolGW, using the TDA and the GTO-type TZVP basis set112 for the lowest 10 singlet excitation states, averaged over all molecules in Thiel’s set. As in Fig. 1, the G0W0 quasiparticle energies for the FHI-aims BSE calculations are here taken from MolGW.

Close modal

We next investigate how the BSE energies are impacted by different treatments of the frequency dependence of the G0W0 self-energy ΣG0W0. In Sec. IV A, we took the G0W0 quasiparticle energies calculated by MolGW as the input for the FHI-aims BSE calculations. G0W0 calculations in MolGW employ an exact analytic treatment of ΣG0W0 on the real axis.60,75 A similarly precise result is expected from the CD technique by Golze et al.98 In this section, we investigate the impact of two frequently employed inexact but potentially cost-saving98 approximations to the self-energy on the real axis, namely, the two-pole approximation,95 Eq. (9), and the Padé approximation,96 Eq. (10), with 16 parameters.

In Fig. 3, we show mean absolute errors MAE(i) of the G0W0 quasiparticle energies of the HOMO-3 to LUMO+3 states, calculated using either the two-pole approximation [Fig. 3(a)] or the 16-parameter Padé approximation [Fig. 3(b)] and compared to the MolGW reference values. The G0W0 results are based on DFT calculations using the PBE34 exchange-correlation functional and employ the Gaussian-type TZVP basis set.112 We see that the two-pole approximation gives MAE values up to 0.3 eV in the investigated states. Although not a small value, this must be viewed in the context of the overall change from plain DFT single-particle energies to correlated G0W0 quasiparticle energies, which are typically of the order of several eV. The 16-parameter Padé approximation can reduce this error by a factor of two or more, as shown in Fig. 3(b). Most of the quasiparticle energies agree with the MolGW results within 0.1 eV, except for the state LUMO+2, where the MAE is around 0.12 eV. The performance of both approximations can be compared with a broader analysis performed in Ref. 97 (the “GW100” benchmark). For a different set of molecules, Ref. 97 reports MAE values of 0.125 eV for HOMO levels and 0.071 eV for LUMO levels when comparing FHI-aims’ two-pole approximation with the 16-parameter Padé approximation. These values are somewhat worse than the precision of the two-pole approximation for our molecular set in Fig. 3(a). Additionally, Fig. 3(a) includes states away from the HOMO and LUMO (not considered in Ref. 97), with somewhat higher MAEs. Given the differences in assessed molecules and levels, the two-pole approximation precision reported in Fig. 3(a) is thus very much in line with that found in Ref. 97. For FHI-aims’ 16-parameter Padé approximation, Ref. 97 finds MAE values of 0.003 eV (HOMO) and 0.006 eV (LUMO) when comparing to an RI-free implementation in the Turbomole code and covering a subset of the smaller molecules in the GW100 benchmark. Although this finding is consistent with our general finding that the 16-parameter Padé approximation significantly improves over the average error associated with the two-pole approximation, the MAE values reported for the present test set in Fig. 3(b) still amount to a few tens of millielectronvolt for the HOMO and LUMO. The remaining discrepancy (though small) between our MAE values and those from the GW100 benchmark could partially be attributable to the different molecular test subsets considered in the assessment of the 16-parameter Padé approximation. Additionally, as also discussed in more depth in Sec. V A of Ref. 97, the Padé approximation can result in serious ambiguities for individual eigenvalues εlGW if there are multiple solutions for the self-consistency condition between the left and the right sides of Eq. (5). In our own assessment, we did not pursue the question of different solutions of Eq. (5) between the 16-parameter Padé approximation in FHI-aims and the numerically exact self-energy in MolGW. While our overall reported MAE values in Fig. 3(b) are evidently small, it is conceivable that multiple solutions of Eq. (5) for individual eigenvalues contribute to the observed MAE, as they would in a practical simulation, where existing implementations typically also do not screen specifically for multiple solutions. In fact, despite its higher MAE, a practical advantage of the two-pole approximation is its relative numerical simplicity and therefore its relative robustness against numerical ambiguities, compared to the Padé approximation. Thus, the two-pole approximated self-energy can be preferable for simple reasons of stability, at the price of reduced precision compared to a formally exact G0W0 self-energy. However, one should bear in mind that the two-pole approximation can in some cases also fail dramatically as was shown in Ref. 97.

FIG. 3.

Impact of two different self-energy extrapolation schemes in FHI-aims on G0W0 quasiparticle energies for the HOMO-3 to LUMO+3 states, compared to the formally exact self-energy treatment in MolGW (reference). Mean absolute errors MAE(i) for Thiel’s molecular benchmark set were computed using the GTO-type TZVP basis set.112 (a) MAE(i) using the two-pole approximation, Eq. (9). (b) MAE(i) using the Padé approximation, Eq. (10), for N = 16.

FIG. 3.

Impact of two different self-energy extrapolation schemes in FHI-aims on G0W0 quasiparticle energies for the HOMO-3 to LUMO+3 states, compared to the formally exact self-energy treatment in MolGW (reference). Mean absolute errors MAE(i) for Thiel’s molecular benchmark set were computed using the GTO-type TZVP basis set.112 (a) MAE(i) using the two-pole approximation, Eq. (9). (b) MAE(i) using the Padé approximation, Eq. (10), for N = 16.

Close modal

As shown in Fig. 4, the different approximate G0W0 self-energy treatments affect the BSE excitation energies, which take the G0W0 quasiparticle energies as input. We compare BSE results based on G0W0 quasiparticle energies calculated using the self-energies of Eqs. (9) and (10) to BSE results based on the MolGW G0W0 eigenvalues with the exact analytic self-energy treatment.60,75 Panels (a) and (b) show the MAE(i) values for the ten lowest singlet and triplet BSE excitation energies, respectively, using the two-pole approximation and averaged over all molecules in Thiel’s set. Panels (c) and (d) show the analogous results obtained using the 16-parameter Padé approximation. We see that the 16-parameter Padé approximation yields smaller MAE(i) (≈0.1 eV across all investigated states) than the two-pole approximation [MAE(i) of 0.1–0.4 eV]. The difference is a direct reflection of the different MAE(i) of the G0W0 quasiparticle energies (Fig. 3), which constitute the input for the BSE calculations.

FIG. 4.

Mean absolute errors MAE(i) of the lowest ten BSE@G0W0@PBE excitation energies, averaged over Thiel’s set of molecules, comparing results using G0W0 eigenvalues from analytically continued self-energies and results using the formally exact G0W0 self-energy in MolGW as a reference. TZVP basis sets112 and the TDA are used in all calculations. (a) MAE(i) of singlet states using G0W0 quasiparticle energies from the two-pole approximation, Eq. (9). (b) MAE(i) of triplet states using the two-pole approximation. (c) MAE(i) of singlet states using G0W0 quasiparticle energies from the 16-parameter Padé approximation, Eq. (10). (d) MAE(i) of triplet states using the 16-parameter Padé approximation.

FIG. 4.

Mean absolute errors MAE(i) of the lowest ten BSE@G0W0@PBE excitation energies, averaged over Thiel’s set of molecules, comparing results using G0W0 eigenvalues from analytically continued self-energies and results using the formally exact G0W0 self-energy in MolGW as a reference. TZVP basis sets112 and the TDA are used in all calculations. (a) MAE(i) of singlet states using G0W0 quasiparticle energies from the two-pole approximation, Eq. (9). (b) MAE(i) of triplet states using the two-pole approximation. (c) MAE(i) of singlet states using G0W0 quasiparticle energies from the 16-parameter Padé approximation, Eq. (10). (d) MAE(i) of triplet states using the 16-parameter Padé approximation.

Close modal

In addition to the MAE of the BSE@G0W0@PBE results between FHI-aims and MolGW, we can also define the state-dependent mean signed error MSE(i) as

MSE(i)=1Nn=1N(Di,nRi,n).
(42)

Just as in Eq. (41), i and n are indices for states and molecules, respectively, averaging over N = 28 molecules in Thiel’s set. Table I summarizes the overall MSE and MAE (averaged over states i in addition to averaging over molecules) of the BSE@G0W0@PBE results using the two-pole and 16-parameter Padé approximation vs MolGW as a reference. As noted above, the MAE of BSE results based on quasiparticle eigenvalues from the 16-parameter Padé approximation is less by a factor of two than the MAE of BSE results from two-pole approximated quasiparticle eigenvalues. The MSE indicates that BSE results from the two-pole quasiparticle eigenvalues overestimate the expected BSE@G0W0 excitation energies based on an exact G0W0 self-energy.

TABLE I.

Mean absolute error and mean signed error of the BSE@G0W0@PBE excitation energies using analytical continuations (two-pole or 16-parameter Padé approximation) for the G0W0 quasiparticle eigenvalues. The values are averaged over the lowest ten states and all molecules in Thiel’s set, using the exact G0W0 self-energy treatment in MolGW as a reference. Results are shown for both singlet and triplet excitation energies.

G0W0Two-pole16-Parameter Padé
self-energySingletTripletSingletTriplet
MSE 0.10 eV 0.15 eV 0.01 eV 0.04 eV 
MAE 0.16 eV 0.22 eV 0.06 eV 0.08 eV 
G0W0Two-pole16-Parameter Padé
self-energySingletTripletSingletTriplet
MSE 0.10 eV 0.15 eV 0.01 eV 0.04 eV 
MAE 0.16 eV 0.22 eV 0.06 eV 0.08 eV 

It is interesting to compare the errors incurred from the different G0W0 self-energy treatments to the errors associated with the BSE approach itself. Bruneval et al.43 show that the BSE@G0W0@PBE singlet excitation energies at the TZVP basis set level give an MSE of −0.8 eV and an MAE of 0.8 eV compared to the BE results of Schreiber et al.24 The MSE and MAE of the BSE@G0W0@PBE triplet excitation energies are around −1.2 eV and 1.2 eV, respectively.43 In comparison, the error incurred through the approximate G0W0 self-energies in Table I is rather small for the two-pole approximation and negligible for the 16-parameter Padé approximation. Additionally, the sign of the MSE from both self-energy approximations is the opposite of that of the MSE compared to the BE values, i.e., especially the two-pole approximation would actually reduce the MSE compared to the BE values, as a result of fortuitous error cancellation. However, this reduction should not be trusted systematically. For a true improvement over the reported small-molecule BSE excitation energies, it would obviously be preferable to pursue higher-level approaches than BSE@G0W0@PBE at the TZVP basis set level—in terms of the DFT starting point, in terms of the theoretical treatment of the neutral excitation, and in terms of the basis set.

Even for the small-molecule systems in Thiel’s set, the BSE correction still accurately captures the vast majority of the change between straight differences of HOMO and LUMO levels from G0W0 calculations and actual neutral excitation energies. This correction is much larger than the differences incurred above as a result of the approximate G0W0 self-energies. In our calculations, the BSE@G0W0@PBE results using both 16-parameter Padé approximation and two-pole approximations reduce the lowest G0W0 HOMO-LUMO gap by about 57%, averaged over all molecules in Thiel’s set. For the He atom, essentially basis set converged results by Li et al.49 show a G0W0@Hartree-Fock (HF) fundamental gap of 24.69 eV, compared to significantly lower lowest-lying triplet and singlet excitation energies of 19.82 eV and 21.22 eV, respectively. Even in this extreme case of an isolated two-electron atom, these essentially basis set converged BSE@G0W0@HF results agree with the exact result better than ≈0.01 Ha (≈0.3 eV). For larger molecules, such as free-base porphyrin and tetraphenylporphyrin,115 BSE@G0W0LDA excitations can match experimental absorption spectra to a similarly good degree (essentially exact within the remaining approximations made in Ref. 115). Here, the BSE again corrects the simple G0W0 HOMO-LUMO energy difference by several eV, much more than the magnitude of changes due to analytical corrections to the G0W0 self-energies reported above.

We now address the basis set convergence of the NAO basis sets for the BSE calculation, in comparison with cc-pVnZ basis sets116 and aug-cc-pVnZ basis sets114 by Dunning and co-workers. Two types of NAO basis sets have been constructed in the context of FHI-aims. The first, denoted “FHI-aims-2009,” was introduced in Ref. 88 and aimed at ground-state DFT calculations. The FHI-aims-2009 basis sets come in different tiers (i.e., levels) n = 1, 2, 3, 4 (in the case of H, a fourth tier does not exist and the molecular results for tier4 below employ tier3 for H). These basis sets allow one to achieve total-energy convergence corresponding to fast qualitative calculations to few-millielectronvolts per atom117 calculations in a single, hierarchical basis set scheme (i.e., for a given element, each basis set level contains the exact basis functions from all lower-accuracy basis set levels as a subset). For first- and second-row elements, FHI-aims’ “tight” settings employ all FHI-aims-2009 basis functions up to and including tier 2, as shown in Table S6 in the supplementary material. The second type of basis set is denoted NAO-VCC-nZ with n = 2, 3, 4, 5.89 The NAO-VCC-nZ basis functions are constructed in analogy to the cc-pVnZ correlation-consistent (CC) basis sets by Dunning,116 but employing the numerically tabulated shape of NAOs (nodeless hydrogen-like radial functions with a numerical confinement potential applied to the tails). The NAO-VCC-nZ basis sets are optimized to be suitable for converging electronic total-energy calculations based on valence-only correlation methods that include sums over unoccupied states, e.g., RPA or MP2.89 In the following, the above two types of NAO basis sets, as well as cc-pVnZ (n = 2, 3, 4, 5)116 and aug-cc-pVnZ (n = 2, 3, 4, 5) basis sets,114 are compared to the aug-cc-pV5Z basis sets as the benchmark reference in the BSE calculations.

In Fig. 5, we show the difference between the BSE excitation energy computed with different basis sets and that of the aug-cc-pV5Z basis set for the lowest five singlet excitations of the ethene molecule in Thiel’s set. The different investigated basis set types and levels are identified on the x axes of all three panels. The size of the different basis sets for the ethene molecule is plotted on the y axis in panel (a). The difference ΔEi between the BSE excitation energy computed using different basis sets and that computed using the aug-cc-pV5Z basis set is identified on the y axes in panels (b) and (c),

ΔEi=EibasisEiaugccpV5Z.
(43)

The ΔEi of the lowest five singlet states (i = 1–5) are plotted in Fig. 5(b) without the TDA and in Fig. 5(c) with the TDA. In both Figs. 5(b) and 5(c), we see that the results obtained with the cc-pVnZ basis sets converge slowly toward the reference value as the basis size increases. The remaining discrepancy is greater than 0.5 eV even for the very expensive cc-pV5Z basis set. Although the magnitude of this discrepancy is unsatisfactory, its occurrence is not surprising. The stated reason for developing the augmented cc basis sets was an improved description of electron affinities,114 a key constituent of the BSE@GW@DFT excitation energies discussed here. Accordingly, the results obtained with the aug-cc-pVnZ basis sets converge much faster. This is in line with the literature:43,44 e.g., in Fig. 3 of Ref. 43, the LR-TDDFT and BSE@G0W0@B3LYP113 excitation energies of the ethene and pyrrole molecules using the aug-cc-pVDZ and aug-cc-pVQZ basis sets show differences less than 0.1 eV for ethene and about 0.2 eV for pyrrole, for both LR-TDDFT and BSE. In Figs. 5(b) and 5(c), the results obtained with the NAO-VCC-nZ basis sets, which are constructed analogously to the cc-pVnZ basis sets, display a similarly unsatisfactory convergence behavior as that of the results from cc-pVnZ basis sets. The other type of NAO basis sets, i.e., the FHI-aims-2009 “tier” basis sets, behaves slightly differently compared to the NAO-VCC-nZ basis sets. The FHI-aims-2009 tier2 basis set improves the BSE excitation energies significantly compared to the FHI-aims-2009 tier1 results, in fact even slightly better than the two unaugmented cc basis set prescriptions. However, the FHI-aims-2009 tier3 and tier4 results are very similar to those obtained using the tier2 basis set, displaying no further significant improvement.

FIG. 5.

Difference between the BSE excitation energy computed with different basis sets and that of the aug-cc-pV5Z basis set for the lowest 5 singlet excitations of the ethene molecule in Thiel’s set. Panel (a) shows the basis size of all basis sets. Panel (b) and panel (c) give the difference of the BSE excitation energy without and with TDA, respectively. The excitation energy calculated using the aug-cc-pV5Z basis set by Dunning and co-workers114 is used as the reference value. The “tier” notation corresponds to the basis set tiers of the “FHI-aims-2009” basis sets,88 either unmodified or with additional augmentation functions from the “aug-cc” basis sets as described in the text.

FIG. 5.

Difference between the BSE excitation energy computed with different basis sets and that of the aug-cc-pV5Z basis set for the lowest 5 singlet excitations of the ethene molecule in Thiel’s set. Panel (a) shows the basis size of all basis sets. Panel (b) and panel (c) give the difference of the BSE excitation energy without and with TDA, respectively. The excitation energy calculated using the aug-cc-pV5Z basis set by Dunning and co-workers114 is used as the reference value. The “tier” notation corresponds to the basis set tiers of the “FHI-aims-2009” basis sets,88 either unmodified or with additional augmentation functions from the “aug-cc” basis sets as described in the text.

Close modal

The results discussed so far confirm the importance of the augmentation functions. We thus extend the FHI-aims-2009 tier2 basis set with different numbers of Gaussian-type augmentation functions obtained from the aug-cc-pV5Z basis set. The label “tier2+aug1” in Fig. 5 denotes the basis set generated by adding the first Gaussian augmentation function from the aug-cc-pV5Z basis set (angular momentum quantum number l = 0) to the FHI-aims-2009 tier2 basis set. Similarly, the label “tier2+aug2” denotes the basis set obtained by combining the FHI-aims-2009 tier2 basis set with the first two Gaussian augmentation functions (l = 0, 1) from the aug-cc-pV5Z basis set. We see that the addition of the augmentation functions significantly improves the accuracy of the FHI-aims-2009 tier2 results compared to the reference aug-cc-pV5Z values. Specifically, the tier2+aug2 basis set already yields essentially basis set converged values for the excitation energies of ethene shown in Fig. 5, with a remaining discrepancy of 0.1 eV or less compared to aug-cc-pV5Z. This conclusion is independent of whether or not the TDA is used in the BSE calculation, as shown in Fig. 5(c). As an important main result, the tier2+aug2 basis sets can here provide rather well converged values, comparable to the aug-cc-pV5Z reference values for low-lying neutral excitations. As will be shown below, this result can be generalized to the remainder of the molecules in Thiel’s set. Interestingly, the tier2+aug basis sets thus provide a recipe allowing one to use a basis set that is precise but affordable for ground-state DFT117 and, with a very limited modification, sufficient to achieve highly converged BSE results for low-lying excitations at the same time. For the lowest-energy excitations, which are often those of the greatest interest, we can thus use a very similar NAO basis set prescription as in ground-state DFT.

It is instructive to analyze the origins of the BSE@G0W0@PBE basis set convergence behavior in Fig. 5 in terms of the underlying electronic structure steps. We therefore next discuss the convergence of the HOMO and LUMO levels of ethene as well as of the HOMO-LUMO gap for DFT-PBE (Fig. S1 in the supplementary material) and for G0W0@PBE (Fig. S2 in the supplementary material). Not surprisingly, the DFT-PBE results are all significantly better converged than the BSE@G0W0@PBE results in Fig. 5. The LUMO converges slightly less well than the HOMO, not unexpectedly since it lies only a bit more than 1 eV below the vacuum level and is thus weakly bound. Much in line with the motivation of Kendall et al.114 to address electron affinities in their construction of the augmented cc basis sets, the addition of two augmentation functions converges the LUMO to the visual resolution of Fig. S1(b), i.e., a scale of roughly 0.01 eV. As is well known,68 the convergence of G0W0 quasiparticle energies is much slower with basis set size than that of the DFT single-particle energies. All the correlation consistent basis sets show noticeable improvements of the quasiparticle HOMO, LUMO, and gap values in Fig. S2, up to the largest basis sets investigated here. However, they do not yet reach visual convergence for the LUMO and, thus, for the gap, even for the largest basis sets. The same is true for the FHI-aims-2009 “tier” basis sets, although this group of basis sets was constructed for ground-state DFT88 and is thus, technically, not designed to converge explicitly correlated methods with the same degree of systematicity as the cc basis sets.89 The augmentation functions added to “tier2” do not contribute to the convergence of the quasiparticle HOMO, whereas they do contribute somewhat to the convergence of the quasiparticle LUMO. In the cases of cc basis sets, basis set extrapolation strategies could be employed as done in Ref. 97 for G0W0 quasiparticle levels. However, the BSE@G0W0@PBE results using augmented basis sets in Fig. 5 appear to converge well without extrapolation, which is initially surprising given the slow convergence of the quasiparticle energies. Regarding BSE@G0W0@PBE convergence, the most striking insight arises from the quasiparticle gap behavior in Fig. S2(c). The gap for all nonaugmented basis sets converges slowly, as seen in Fig. 5. In contrast, the slow convergence of quasiparticle HOMO and LUMO cancels nearly exactly for the augmented cc basis sets as well as for the tier2+aug2 and tier2+aug3 basis sets. Since G0W0 quasiparticle energy differences are key terms in BSE@G0W0 excitation energies, the quasiparticle gap behavior in Fig. S2(c) indicates that, at least for ethene, both error cancellation between HOMO and LUMO G0W0 correlation energies and the better description of the LUMO by the augmentation functions contribute to low-lying BSE@G0W0 excitation energies. As a result, the latter is largely numerically converged already for the relatively small tier2+aug2 basis set as well as for the augmented cc basis sets.

Fig. 6 shows the convergence of the MAE(i) of the five lowest-lying excitation energies as a function of the size of various basis sets in BSE calculations, employing the TDA and averaged over all molecules in Thiel’s set. The reference is, again, the aug-cc-pV5Z basis set. The different investigated basis set types and levels are identified on the x axes of panels (a) and (b). The basis size N¯basis of different basis sets, averaged over all molecules in Thiel’s set, is plotted as the y axis in Fig. 6(a). N¯basis is defined as

N¯basis=1Nmolj=1NmolNbasisj,
(44)

where Nmol is the number of molecules in Thiel’s set and Nbasisj is the basis size for the molecule j in the benchmark set. We see that the convergence of different basis sets is similar to the earlier observations made for ethene in Fig. 5. Specifically, the tier2+aug2 basis set produces rather well converged results for all molecules investigated here. The lowest five singlet excitation energies calculated by BSE@G0W0@PBE using the aug-cc-pV5Z and tier2+aug2 basis sets are also listed in Table S7 in the supplementary material for all molecules in Thiel’s set.

FIG. 6.

(a) Size of all basis sets shown, averaged over all molecules in Thiel’s set [see Eq. (44)]. (b) MAE(i) values of the lowest five BSE@G0W0@PBE excitation energies computed with different basis sets, using results from the aug-cc-pV5Z basis set as the reference. The TDA is employed.

FIG. 6.

(a) Size of all basis sets shown, averaged over all molecules in Thiel’s set [see Eq. (44)]. (b) MAE(i) values of the lowest five BSE@G0W0@PBE excitation energies computed with different basis sets, using results from the aug-cc-pV5Z basis set as the reference. The TDA is employed.

Close modal

In the BSE calculations presented in this work so far, we include the orbital pairs of all occupied and unoccupied orbitals in the construction of the BSE matrix. The dimension of the matrix problem Eq. (16) thus grows quadratically with the basis set size and also (for a fixed basis set level) with molecular size. Solving this full BSE matrix problem produces an excitation spectrum of the studied molecule that includes very high excitation energies. In many practical applications, however, we are interested in only the low-lying part of the excitation spectrum. In such a situation, one might suspect that the high-lying unoccupied quasiparticle states do not contribute much to the low-lying optical excitations. This can, in fact, not be entirely true, since single-quasiparticle like GW observables such as the ionization potential can depend significantly on high-lying parts of the spectrum being included in unoccupied-state sums in the relevant perturbation expressions.118 All-electron approaches to G0W0 band gaps suffer from similar convergence issues with basis set size, specifically the basis set resolution in those regions of space that are closer to a nucleus.119 However, neutral excitations are not the same objects, and the question thus remains how many quasiparticle states, especially the high-lying unoccupied orbitals, should be included in the construction of the BSE matrix, in order to obtain sufficiently precise results for low-lying optical excitation energies. Previous investigations on leaving away states in the calculation of excitation energies exist. For instance, this was done in Ref. 115, already referenced above. As another common example, various studies show how to efficiently select the desirable orbital space in the study of complete active space approaches.4,120–122 As a final example, in the calculation of the electronic spectra of molecules in solution or surfaces, Besley developed an approach within LR-TDDFT and single excitation configuration interaction that limits electronic excitations to include only those between orbitals localized on the solute or adsorbent.123 

In Fig. 7, we show the errors incurred in the BSE low-lying singlet and triplet excitation energies obtained when applying different values of a cutoff energy Ecut for unoccupied states, limiting the number of states a and b entering the matrix construction in Eqs. (17) and (18). Here, Ecut denotes a cutoff energy above which the high-lying unoccupied quasiparticle states are omitted from the BSE matrix (however, such cutoffs are not applied in the construction of the quantities entering the BSE matrix). The average numbers of unoccupied states included for different choices of Ecut are tabulated in Table II. Figure 7 shows MAE(i) values for the lowest ten singlet (a) or triplet (b) excitation energies for the different cutoff energy values, using the tier2+aug2 basis set for all calculations and the excitation energies from the full calculation (no cutoff imposed) as a reference. In these calculations, all occupied states are kept in the construction of the BSE matrix, i.e., no cutoff threshold is applied to the occupied quasiparticle states. Figure 7 shows that the error of the results calculated with Ecut = 20 eV is about 20–30 meV. Setting Ecut = 40 eV yields excitation energies with an error closer to 10 meV. Larger Ecut values of 60 and 80 eV lead to further small improvements. In view of the remaining errors of these calculations (due to the level of theory for underlying DFT, quasiparticle energies, and neutral excitation formalism itself), these results suggest that one may use 40 eV as a threshold beyond which the impact of high-lying unoccupied quasiparticle states becomes negligible in BSE calculations of low-lying excitations. This can reduce the computational effort significantly because of the reduction of the number of states needed in construction of the BSE matrix. Specifically, the time complexity and memory complexity of constructing the BSE matrix in Eqs. (17) and (18) scale as Nocc2×Nunocc2, where Nocc and Nunocc denote the number of occupied and unoccupied (g)KS single-particle states, respectively. By setting Ecut = 40 eV and for the tier2+aug2 basis sets, the number of unoccupied states is reduced to about 1/3 of the analogous number if no threshold energy values are used (see Table II). Additionally, the formal effort for solving the full BSE, Eq. (16), would scale as O(N6), where N is a measure of system size, due to the same considerations of how Nocc and Nunocc impact the matrix dimension. While imposing Ecut will not reduce the formal scaling, the actual computational effort will nevertheless be reduced substantially in the limit of large systems where the BSE solution must eventually dominate. In short, both the time and the memory requirements of the BSE calculation of low-lying excitations are expected to be reduced significantly by imposing Ecut, without sacrificing much accuracy.

FIG. 7.

(a) MAE(i) of the BSE@G0W0@PBE excitation energies of the lowest ten singlet states using different energy cutoff values Ecut (20, 40, 60, 80 eV) for unoccupied states, using the case where all unoccupied states are included as a reference. The basis set is the tier2+aug2 basis set, and the values are averaged over all molecules in Thiel’s set. (b) Analogous MAE(i) for triplet states.

FIG. 7.

(a) MAE(i) of the BSE@G0W0@PBE excitation energies of the lowest ten singlet states using different energy cutoff values Ecut (20, 40, 60, 80 eV) for unoccupied states, using the case where all unoccupied states are included as a reference. The basis set is the tier2+aug2 basis set, and the values are averaged over all molecules in Thiel’s set. (b) Analogous MAE(i) for triplet states.

Close modal
TABLE II.

The number of unoccupied states, Nunocc, averaged over all molecules in Thiel’s set, when imposing threshold values Ecut = 20, 40, 60, and 80 eV or no threshold (“…”) in the BSE@G0W0@PBE calculations. The tier2+aug2 basis set is used.

Threshold (eV)20406080
Nunocc 64 95 117 142 292 
Threshold (eV)20406080
Nunocc 64 95 117 142 292 

In this section, we first investigate the basis set convergence of LR-TDDFT excitation energies for the molecules in Thiel’s set, using the strategy already employed for the BSE in Sec. IV C. The LR-TDDFT in FHI-aims was implemented following Eqs. (25)–(30). The exchange-correlation kernel used is the LDA, employing the parametrization of the correlation energy by Perdew and Wang,33,94 provided by the Libxc library.124,125 Note that the LR-TDDFT formalism leaves one with the freedom to choose different prescriptions of the XC functional for (i) the self-consistent solutions of single-particle orbitals and energies and (ii) the XC kernel fxc [Eq. (30)] used in the actual LR-TDDFT construction. In this section, the exchange correlation functional for the self-consistent solutions of single-particle orbitals is PBE.34 We will use the notation “LR-TDDFT-LDA@PBE” in the following to denote this approach. To maintain consistency, the TDA is employed in the LR-TDDFT calculations shown in the main text of this paper, as this is widely done for BSE calculations as well. Additionally, LR-TDDFT calculations without making the TDA are included in the SM.

Figure 8 provides state-dependent MAE(i) values for the lowest five LR-TDDFT excitation energies computed with different basis sets, averaged over Thiel’s set. Excitation energies computed with the aug-cc-pV5Z basis sets are used as reference values. The different investigated basis set types and levels are identified on the x axes of both panels (a) and (b). Figure 8(a) for singlet states and Fig. 8(b) for triplet states show essentially identical behavior. The overall convergence pattern associated with all the basis sets investigated here is also very similar to the behavior seen for the BSE in Fig. 6. Excitation energy values derived from both the cc-pVnZ basis sets and the NAO-VCC-nZ basis sets converge slowly toward the reference value as the basis size increases. The largest error is about 0.5 eV for the very expensive cc-pV5Z basis set for singlet excitation energies and 0.3 eV for triplet excitation energies. The aug-cc-pVnZ basis sets, again, converge much faster. As for the BSE, the FHI-aims-2009 basis sets improve significantly toward the reference results from tier1 to tier2, but not further using tier3 and tier4. Finally, by adding augmentation functions to the FHI-aims-2009 tier2 basis sets, one can obtain well converged LR-TDDFT results by including only one or two augmentation basis functions. The tier2+aug2 basis set, already discussed for BSE calculations above, leads to compellingly good basis set convergence as a production recipe. The lowest five singlet and triplet excitation energies calculated by LR-TDDFT-LDA@PBE using aug-cc-pV5Z and tier2+aug2 basis sets are provided in Tables S8 and S9 in the supplementary material for all molecules in Thiel’s set.

FIG. 8.

(a) MAE(i) of LR-TDDFT-LDA@PBE excitation energies for the lowest five singlet excitation states, averaged for all molecules in Thiel’s set, computed with different basis sets and using the aug-cc-pV5Z basis set as a reference. (b) Analogous MAE(i) for triplet states. The TDA is employed in both subfigures.

FIG. 8.

(a) MAE(i) of LR-TDDFT-LDA@PBE excitation energies for the lowest five singlet excitation states, averaged for all molecules in Thiel’s set, computed with different basis sets and using the aug-cc-pV5Z basis set as a reference. (b) Analogous MAE(i) for triplet states. The TDA is employed in both subfigures.

Close modal

As pointed out by a reviewer of this work, it is initially somewhat surprising that the convergence of LR-TDDFT and of BSE@G0W0 is qualitatively so similar since LR-TDDFT contains no quasiparticle energies and no screened Coulomb interaction. Indeed, in our discussion of ethene above (Fig. 5 as well as Figs. S1 and S2 in the supplementary material), we traced the slow convergence of BSE@G0W0 excitation energies to the quasiparticle energies as well as the need for augmentation functions to better capture the quasiparticle LUMO level. However, the underlying concepts of BSE and LR-TDDFT are of course similar—see, for instance, Ref. 22 for a detailed derivation. Indeed, the derivation of LR-TDDFT is based on the linear response function χ of the system [see Eq. (5.7) and following in Ref. 22 for details], and the linear response χ0 of the underlying generalized Kohn-Sham system is common to both BSE@G0W0 and LR-TDDFT. It was shown by Bruneval126 that treating only χ0 in a larger basis set can enhance the convergence of G0W0 quasiparticle eigenvalues, indicating that this is the quantity that determines the common convergence behavior of both approaches.

We next compare the results of BSE@G0W0@PBE with those of LR-TDDFT-LDA@PBE, implemented in FHI-aims and using the tier2+aug2 basis set validated above. While the primary focus of the present work is numerical and basis set convergence, we provide this comparison here because LR-TDDFT is widely used in quantum chemistry as a computationally efficient method for optical excitation calculations. We note that similar comparisons can be found in the literature,43,44 albeit not using the same basis sets. In our comparison, the underlying DFT calculations for both BSE and LR-TDDFT are carried out using the PBE34 exchange-correlation functional. Figure 9 shows the MSE(i) and the MAE(i) between LR-TDDFT-LDA@PBE results and BSE@G0W0@PBE results, averaged over all molecules in Thiel’s set. The BSE@G0W0@PBE results are taken as the reference. The lowest ten singlet and triplet excitation energies are compared. The MSE between LR-TDDFT-LDA@PBE and BSE@G0W0@PBE results for singlet and triplet excitation states is plotted in Figs. 9(a) and 9(b), respectively. The MAE between LR-TDDFT-LDA@PBE and BSE@G0W0@PBE results for singlet and triplet excitation states is plotted in Figs. 9(c) and 9(d), respectively. We see from Fig. 9 that the MAE between LR-TDDFT-LDA@PBE and BSE@G0W0@PBE for singlet excitation energies is less than 0.5 eV, whereas the MAE for triplet excitations energies can be as large as 0.7 eV. Singlet excitation energy values obtained from LR-TDDFT-LDA@PBE tend to be lower than those obtained from BSE@G0W0@PBE, whereas LR-TDDFT-LDA@PBE triplet excitation energies appear to be larger by up to 0.7 eV than those obtained from BSE@G0W0@PBE. Previous studies of BSE and LR-TDDFT show that the excitation energies computed by BSE and LR-TDDFT depend highly on the DFT starting point.43,44 Bruneval et al.43 compare the MSE of both BSE@G0W0@B3LYP and LR-TDDFT@B3LYP excitation energies with the BEs of Thiel’s set, showing that the MSE of LR-TDDFT is about 0.4 eV lower than the MSE of BSE.26,43 However, there are several differences between their comparison and the comparison shown in this work. First, the dataset used by Bruneval and co-workers is the BEs of Thiel’s set,24 which contains 103 singlet and 63 triplet excitation energies. In our work, we have a larger dataset to include the lowest ten singlet and triplet states of each molecule in Thiel’s set. Second, the BSE and LR-TDDFT calculations analyzed in this section rely on the TDA, which is not employed in the comparison performed by Bruneval et al.43 Finally, we here use a basis set that is essentially converged for both BSE and LR-TDDFT calculations. Another set of comparable results is therefore those of Jacquemin and co-workers,44 who compare BSE@G0W0@PBE0 and LR-TDDFT@PBE0 in a benchmark paper using the aug-cc-pVTZ basis set, which has similar convergence behavior as the tier2+aug2 basis set used here. Different MAE values between BSE@G0W0@PBE0 and LR-TDDFT@PBE0 excitation energies are reported for different categories of molecules in Thiel’s set: 0.27 eV for unsaturated aliphatic hydrocarbons; 0.51 eV for aromatic compounds; 0.37 eV for aldehydes, ketones, and amides; and 0.47 eV for nucleobases. The reported values are comparable to the values that we find in Fig. 9, in the range of 0.2–0.7 eV for different states of singlet and triplet excitation energies.

FIG. 9.

Mean signed errors MSE(i) and mean absolute errors MAE(i) between LR-TDDFT-LDA@PBE and BSE@G0W0@PBE (reference) results, averaged for all molecules in Thiel’s set, using the tier2+aug2 basis set in the TDA for the lowest ten singlet and triplet states. Panels (a) and (b) show the MSE of the singlet and triplet states, respectively. Panels (c) and (d) show the MAE of the singlet and triplet states, respectively.

FIG. 9.

Mean signed errors MSE(i) and mean absolute errors MAE(i) between LR-TDDFT-LDA@PBE and BSE@G0W0@PBE (reference) results, averaged for all molecules in Thiel’s set, using the tier2+aug2 basis set in the TDA for the lowest ten singlet and triplet states. Panels (a) and (b) show the MSE of the singlet and triplet states, respectively. Panels (c) and (d) show the MAE of the singlet and triplet states, respectively.

Close modal

Finally, we briefly return to the overall accuracy of both BSE@G0W0@PBE and LR-TDDFT-LDA@PBE compared to BE24 values. In the literature and for TZVP basis sets, the reported error for BSE@G0W0@PBE can be around 0.8 eV (singlets) or 1.2 eV (triplets),43 as already alluded to in Sec. IV B. In Fig. 10, we show results using our own recommended basis set, tier2+aug2, comparing to BE values for the lowest singlets and (where available) triplets for all molecules in Thiel’s set. Interestingly, the BSE@G0W0@PBE energies appear to deviate even somewhat more from the BE values as a result of the, arguably improved, basis set than those reported for the TZVP basis set.43 In comparison, the LR-TDDFT-LDA@PBE level of theory performs better for this set of molecules [Fig. 10(b)], shown using the TDA, albeit still off on average by several tenths of electronvolts. This observation does not depend on the TDA, as shown in Fig. S3 in the SM, which compares LR-TDDFT-LDA@PBE data with and without the TDA to the same set of BE data as used in Fig. 10. As noted earlier, the primary objective of this paper is not a discussion of the accuracy compared to the BE value but rather the assessment of basis set convergence and other methodological aspects of our approach. Nevertheless, the results shown in Fig. 10 should be borne in mind by users of the approaches. Significant improvements to both BSE and LR-TDDFT for molecules are of course possible and documented in the literature, e.g., by changing the underlying density functional approximation44 and/or (for BSE) the underlying level of the GW approach applied.23,51

FIG. 10.

Deviations of the BSE@G0W0@PBE and LR-TDDFT-LDA@PBE results for the lowest singlet and triplet excitation energies of each molecule in Thiel’s set, compared with the BE values.24 Panels (a) and (b) show the comparison for BSE and LR-TDDFT, respectively. The molecule indices follow the order in the original paper by Schreiber et al.24 The BE values for the lowest triplet state of molecules 14–18 and 25–28 are not available; therefore, the comparison of the lowest triplet state is omitted for these molecules. The TDA is employed in both subfigures.

FIG. 10.

Deviations of the BSE@G0W0@PBE and LR-TDDFT-LDA@PBE results for the lowest singlet and triplet excitation energies of each molecule in Thiel’s set, compared with the BE values.24 Panels (a) and (b) show the comparison for BSE and LR-TDDFT, respectively. The molecule indices follow the order in the original paper by Schreiber et al.24 The BE values for the lowest triplet state of molecules 14–18 and 25–28 are not available; therefore, the comparison of the lowest triplet state is omitted for these molecules. The TDA is employed in both subfigures.

Close modal

While the present work does not include a performance optimization of either the BSE or the LR-TDDFT implementations reported above, we provide some individual timings indicative of the computational cost of our (not fully optimized) BSE and LR-TDDFT implementations. The results should only be understood as qualitative indicators of the current implementation, since no dedicated optimization was carried out, but nevertheless give some idea of the relative cost of different steps at present and indicate avenues for the future work in our own implementation. In this section, we consider a series of acene molecules of increasing size: benzene, naphthalene, and anthracene. Both BSE and LR-TDDFT calculations were performed using the “tier2+aug2” basis sets described in Secs. IV C and IV E, respectively, and employing the TDA.

In Table III, we show the timings of BSE@G0W0@PBE and LR-TDDFT-LDA@PBE performed for the three acene molecules selected. We apply a cutoff energy Ecut = 40 eV to limit the number of unoccupied states entering the BSE and Casida matrix construction as described in Sec. IV D. The calculations are performed using a single node with 44 cores (Intel Xeon, Broadwell microarchitecture, 2.4 GHz) on the Dogwood cluster at University of North Carolina, Chapel Hill. The total time is decomposed into three parts, i.e., RI basis setup denoted by “Basis” (precomputed three-center and two-center integrals in Sec. III), BSE/LR-TDDFT matrix construction or building denoted by “Build Mat.,” and solving the BSE/LR-TDDFT matrix as eigenvalue problems in Eqs. (22) and (25) denoted by “Solver.” The “Build Mat.” timing accounts for the computational effort in building the BSE matrix as outlined in Eqs. (17) and (18) vs the LR-TDDFT matrix in Eq. (26). We note that the timings comprise the BSE/LR-TDDFT step only, whereas the timings for the preceding steps are not counted, such as the G0W0 and DFT-PBE steps in the BSE@G0W0@PBE calculation and the DFT-PBE step in the LR-TDDFT-LDA@PBE calculation. Table III shows that the majority of the total timing is attributed to the product basis setup step for both BSE and LR-TDDFT calculations. Both the total timing and the timings for each step are comparable for BSE and LR-TDDFT calculations, which is expected due to the similar formalisms. However, as will be seen below in Table IV, the difference of timings for the matrix building step between BSE and TDDFT becomes more significant in our present implementations if the matrix size is increased.

TABLE III.

Timings of BSE@G0W0@PBE and LR-TDDFT-LDA@PBE calculations performed for the neutral singlet excitations of benzene, naphthalene, and anthracene. The “tier2+aug2” basis set is used, and a cutoff energy Ecut = 40 eV for unoccupied states is applied. Calculations are performed using one node with 44 cores on the Dogwood cluster at UNC-Chapel Hill (Intel Broadwell microarchitecture, 2.4 GHz). All values are given in seconds (s). The TDA is employed for BSE but not for LR-TDDFT; however, the matrix dimensions are identical in Casida’s Equation vs. in BSE(TDA).

BSE@G0W0@PBELR-TDDFT-LDA@PBE
MoleculesBasisBuild Mat.SolverTotalBasisBuild Mat.SolverTotal
Benzene 270 280 293 301 
Naphthalene 1200 44 1249 1184 39 1228 
Anthracene 3089 142 16 3247 2907 129 19 3055 
BSE@G0W0@PBELR-TDDFT-LDA@PBE
MoleculesBasisBuild Mat.SolverTotalBasisBuild Mat.SolverTotal
Benzene 270 280 293 301 
Naphthalene 1200 44 1249 1184 39 1228 
Anthracene 3089 142 16 3247 2907 129 19 3055 
TABLE IV.

Timings as in Table III, but without any energy cutoff applied for unoccupied states.

BSE@G0W0@PBELR-TDDFT-LDA@PBE
MoleculesBasisBuild Mat.SolverTotalBasisBuild Mat.SolverTotal
Benzene 276 31 313 292 22 321 
Naphthalene 1205 314 72 1592 1200 173 71 1443 
Anthracene 3175 1531 407 5113 3010 824 405 4239 
BSE@G0W0@PBELR-TDDFT-LDA@PBE
MoleculesBasisBuild Mat.SolverTotalBasisBuild Mat.SolverTotal
Benzene 276 31 313 292 22 321 
Naphthalene 1205 314 72 1592 1200 173 71 1443 
Anthracene 3175 1531 407 5113 3010 824 405 4239 

To demonstrate how the cutoff energy Ecut = 40 eV reduces the computational timings, we show in Table IV the timings for the same systems without using any cutoff energy for both BSE@G0W0@PBE and LR-TDDFT-LDA@PBE calculations. All other computational details remain the same as those used to obtain Table III. We see in Table IV that, by using all the unoccupied states in the BSE and LR-TDDFT calculations, the timings for building and solving matrix steps increase significantly compared to those in Table III. In contrast, the timings for the product basis setup step stay almost unchanged, which is expected since the number of basis functions and auxiliary basis functions is not affected by applying the cutoff energy to limit the number of unoccupied states. As mentioned above, the matrix building step in the present implementation reveals a difference between the BSE and the LR-TDDFT timings in Table IV. Specifically, the matrix building step for LR-TDDFT requires the calculations of the kernel Kia,jb [Eqs. (26) and (27)]. In our parallel implementation, the state indices i, a, j, and b are distributed among different processors, as are the RI two- and three-center integrals, and the calculation of Kia,jb requires significant interprocessor communication to get the correct state indices. For BSE, the analogous interprocessor communication has to be conducted twice in the calculations of (ia|V | jb) and (ij|W|ab) because the state index order is different in (ia|V | jb) and (ij|W|ab), and thus, they are calculated in separate steps. As a result, the timing of the matrix building step in the BSE is significantly larger than that in the LR-TDDFT calculations.

Finally, we verify how the cutoff energy Ecut = 40 eV reduces the memory requirements for the BSE and LR-TDDFT matrices. In Table V, we compare the memory used to store the BSE or LR-TDDFT matrix using all unoccupied states and the memory requirements for the matrices when applying Ecut = 40 eV, indicating a reduction by a factor of ∼9–11. This is consistent with that in Sec. IV D, where Ecut = 40 eV was shown to reduce the number of unoccupied states to about 1/3 of the full number of states if no cutoff energy values are used. Since the BSE and LR-TDDFT matrix size scales as Nocc2×Nunocc2 [Eqs. (17) and (18)], the overall reduction amounts to (1/3)2 = 1/9, which is confirmed by the memory reduction shown in Table V.

TABLE V.

Memory (in Megabytes) used to store the BSE or LR-TDDFT matrix in the BSE@G0W0@PBE or LR-TDDFT-LDA@PBE calculations of neutral singlet excitations of benzene, naphthalene, and anthracene. The column labeled “Full” denotes the memory required if all unoccupied states enter the BSE or LR-TDDFT matrix. The column labeled “Ecut = 40 eV” denotes the memory required when applying the cutoff energy of 40 eV for unoccupied states. The column labeled “Ratio” denotes the ratio of the value in the “Full” column over the value in the “Ecut = 40 eV” column. Other computational details are the same as those in Tables III and IV.

BSE@G0W0@PBELR-TDDFT-LDA@PBE
MoleculesFullEcut = 40 eVRatioFullEcut = 40 eVRatio
Benzene 400 45 8.9 400 49 8.2 
Naphthalene 2549 261 9.8 2549 282 9.0 
Anthracene 9129 876 10.4 9129 979 9.3 
BSE@G0W0@PBELR-TDDFT-LDA@PBE
MoleculesFullEcut = 40 eVRatioFullEcut = 40 eVRatio
Benzene 400 45 8.9 400 49 8.2 
Naphthalene 2549 261 9.8 2549 282 9.0 
Anthracene 9129 876 10.4 9129 979 9.3 

We describe an implementation of the Bethe-Salpeter equation approach to neutral excitations in small molecules based on numeric atom-centered basis sets in an all-electron electronic structure framework (the FHI-aims code). Benchmarks performed using Thiel’s set of small molecules24 demonstrate the numerical correctness of the implementation. Mean absolute errors less than 1 meV are obtained compared to the reference values computed using the MolGW code when exactly the same basis set and underlying technical approximations are used. We next investigate the impact of analytical approximations to the G0W0 self-energy (the two-pole and 16-parameter Padé approximations) on the G0W0 quasiparticle energies entering the BSE. The MAE of the BSE@G0W0@PBE results with these analytical approximations is around 0.05–0.20 eV, compared with the exact G0W0 self-energy used in the MolGW reference code. The 16-parameter Padé approximation is more precise than the two-pole approximation where it can be used, but (due to the possible occurrence of multiple solutions in the quasiparticle equation, discussed elsewhere97) the two-pole approximation offers an overall numerically more stable avenue (with the caveats laid out in Ref. 97). Ultimately, the differences due to either approximation are smaller than typical basis set errors and errors due to the level of theory itself as assessed in other benchmark publications. The basis set convergence behavior of the predicted low-lying excitations is investigated for the cc-pVnZ,116 FHI-aims-2009,88 NAO-VCC-nZ,89 and aug-cc-pVnZ114 literature basis sets, as well as for a simple modification of the FHI-aims-2009 tier2 basis set by adding two augmentation functions from the aug-cc basis sets,114 called “tier2+aug2.” For both BSE@G0W0@PBE and LR-TDDFT-LDA@PBE, adequate precision requires the use of augmentation functions, as expected from the literature. The “tier2+aug2” basis set provides high precision for both BSE and LR-TDDFT calculations while remaining applicable in production calculations. Furthermore, the convergence is investigated with respect to the number of unoccupied states included in the BSE or LR-TDDFT matrix construction. A threshold of Ecut = 40 eV is suggested, above which the unoccupied states are discarded in the construction of either the BSE or the LR-TDDFT matrix. This threshold significantly reduces the time and memory consumption while maintaining high precision of the results for low-lying excitations. Finally, BSE@G0W0@PBE and LR-TDDFT-LDA@PBE results are compared using the tier2+aug2 basis set for Thiel’s set of molecules. The difference between BSE@G0W0@PBE and LR-TDDFT-LDA@PBE is quantified by an MAE in the range of 0.2–0.7 eV for different singlet and triplet states calculated for molecules in Thiel’s set. In agreement with the literature, deviations from higher-level “best estimate” values are of a similar magnitude; one likely mitigation strategy is the selection of a better starting-point density functional for BSE@G0W0@DFT.

See the supplementary material for excitation energies of the molecules in Thiel’s set using the BSE with and without the TDA; corresponding oscillator strengths of singlets within the TDA; definition of the “tier2” basis sets and numerical settings; DFT and G0W0 single-particle HOMO, LUMO, and gap data for the ethene molecule; and basis set convergence of excitation energies of the molecules in Thiel’s set using the BSE and LR-TDDFT; comparison of LR-TDDFT results with and without the TDA to BE data. We furthermore provide numerical versions of the data shown in the supplementary material, Tables S1–S5 and S7–S9, in a comma-separated value (csv) format as file “si_tab.csv.”

This work was financially supported by the NSF under Award Nos. DMR-1729297 and DMR-1728921, as well as through the Research Triangle MRSEC (Grant No. DMR-11-21107). An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility (ALCF), which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357. The authors would like to thank the University of North Carolina at Chapel Hill and the Research Computing group for providing computational resources and support that have contributed to these research results. X.R. acknowledges the financial support from the Chinese National Science Foundation (Grant Nos. 11574283 and 11874335) and the Max Planck Partner Group project. The authors are grateful to the reviewers of this work for their careful reading and numerous helpful suggestions that the authors were fortunate to be able to include in this work.

1.
J. F.
Stanton
and
R. J.
Bartlett
, “
The equation of motion coupled-cluster method. A systematic biorthogonal approach to molecular excitation energies, transition probabilities, and excited state properties
,”
J. Chem. Phys.
98
,
7029
7039
(
1993
).
2.
R. J.
Bartlett
and
M.
Musiał
, “
Coupled-cluster theory in quantum chemistry
,”
Rev. Mod. Phys.
79
,
291
352
(
2007
).
3.
H.
Sekino
and
R. J.
Bartlett
, “
A linear response, coupled-cluster theory for excitation energy
,”
Int. J. Quantum Chem.
26
,
255
265
(
1984
).
4.
B. O.
Roos
,
P. R.
Taylor
, and
P. E.
Siegbahn
, “
A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach
,”
Chem. Phys.
48
,
157
173
(
1980
).
5.
K.
Andersson
,
P.-Å.
Malmqvist
, and
B. O.
Roos
, “
Second-order perturbation theory with a complete active space self-consistent field reference function
,”
J. Chem. Phys.
96
,
1218
1226
(
1992
).
6.
J.
Finley
,
P.-A.
Malmqvist
,
B. O.
Roos
, and
L.
Serrano-Andrés
, “
The multi-state CASPT2 method
,”
Chem. Phys. Lett.
288
,
299
306
(
1998
).
7.
G.
Ghigo
,
B. O.
Roos
, and
P.
Malmqvist
, “
A modified definition of the zeroth-order Hamiltonian in multiconfigurational perturbation theory (CASPT2)
,”
Chem. Phys. Lett.
396
,
142
149
(
2004
).
8.
J. C.
Grossman
,
M.
Rohlfing
,
L.
Mitas
,
S. G.
Louie
, and
M. L.
Cohen
, “
High accuracy many-body calculational approaches for excitations in molecules
,”
Phys. Rev. Lett.
86
,
472
475
(
2001
).
9.
M. L.
Tiago
,
P. R. C.
Kent
,
R. Q.
Hood
, and
F. A.
Reboredo
, “
Neutral and charged excitations in carbon fullerenes from first-principles many-body theories
,”
J. Chem. Phys.
129
,
084311
(
2008
).
10.
D.
Lee
,
J. L.
DuBois
, and
Y.
Kanai
, “
Importance of excitonic effect in charge separation at quantum-dot/organic interface: First-principles many-body calculations
,”
Nano Lett.
14
,
6884
6888
(
2014
).
11.
N. S.
Blunt
,
S. D.
Smart
,
G. H.
Booth
, and
A.
Alavi
, “
An excited-state approach within full configuration interaction quantum Monte Carlo
,”
J. Chem. Phys.
143
,
134117
(
2015
).
12.
E.
Runge
and
E. K. U.
Gross
, “
Density-functional theory for time-dependent systems
,”
Phys. Rev. Lett.
52
,
997
1000
(
1984
).
13.
M.
Marques
and
E.
Gross
, “
Time-dependent density functional theory
,”
Annu. Rev. Phys. Chem.
55
,
427
455
(
2004
).
14.
M. E.
Casida
, “
Time-dependent density-functional theory for molecules and molecular solids
,”
J. Mol. Struct.
914
,
3
18
(
2009
).
15.
M. E.
Casida
,
C.
Jamorski
,
K. C.
Casida
, and
D. R.
Salahub
, “
Molecular excitation energies to high-lying bound states from time-dependent density-functional response theory: Characterization and correction of the time-dependent local density approximation ionization threshold
,”
J. Chem. Phys.
108
,
4439
4449
(
1998
).
16.
L.
Hedin
, “
New method for calculating the one-particle Green’s function with application to the electron-gas problem
,”
Phys. Rev.
139
,
A796
A823
(
1965
).
17.
E. E.
Salpeter
and
H. A.
Bethe
, “
A relativistic equation for bound-state problems
,”
Phys. Rev.
84
,
1232
1242
(
1951
).
18.
M. S.
Hybertsen
and
S. G.
Louie
, “
Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies
,”
Phys. Rev. B
34
,
5390
5413
(
1986
).
19.
M.
Rohlfing
and
S. G.
Louie
, “
Electron-hole excitations in semiconductors and insulators
,”
Phys. Rev. Lett.
81
,
2312
2315
(
1998
).
20.
M.
Rohlfing
and
S. G.
Louie
, “
Electron-hole excitations and optical spectra from first principles
,”
Phys. Rev. B
62
,
4927
4944
(
2000
).
21.
G.
Strinati
, “
Application of the Green’s functions method to the study of the optical properties of semiconductors
,”
Riv. Nuovo Cimento
11
,
1
86
(
1988
).
22.
G.
Onida
,
L.
Reining
, and
A.
Rubio
, “
Electronic excitations: Density-functional versus many-body Green’s-function approaches
,”
Rev. Mod. Phys.
74
,
601
659
(
2002
).
23.
X.
Blase
,
I.
Duchemin
, and
D.
Jacquemin
, “
The Bethe–Salpeter equation in chemistry: Relations with TD-DFT, applications and challenges
,”
Chem. Soc. Rev.
47
,
1022
1043
(
2018
).
24.
M.
Schreiber
,
M. R.
Silva-Junior
,
S. P. A.
Sauer
, and
W.
Thiel
, “
Benchmarks for electronically excited states: CASPT2, CC2, CCSD, and CC3
,”
J. Chem. Phys.
128
,
134110
(
2008
).
25.
D.
Jacquemin
,
E. A.
Perpéte
,
I.
Ciofini
, and
C.
Adamo
, “
Assessment of functionals for TD-DFT calculations of singlet-triplet transitions
,”
J. Chem. Theory Comput.
6
,
1532
1537
(
2010
).
26.
M. R.
Silva-Junior
,
M.
Schreiber
,
S. P. A.
Sauer
, and
W.
Thiel
, “
Benchmarks for electronically excited states: Time-dependent density functional theory and density functional theory based multireference configuration interaction
,”
J. Chem. Phys.
129
,
104103
(
2008
).
27.
Y.
Imamura
and
H.
Nakai
, “
Time-dependent density functional theory (TDDFT) calculations for core-excited states: Assessment of an exchange functional combining the Becke88 and van Leeuwen–Baerends-type functionals
,”
Chem. Phys. Lett.
419
,
297
303
(
2006
).
28.
R.
Baer
and
D.
Neuhauser
, “
Density functional theory with correct long-range asymptotic behavior
,”
Phys. Rev. Lett.
94
,
043002
(
2005
).
29.
S.
Refaely-Abramson
,
S.
Sharifzadeh
,
N.
Govind
,
J.
Autschbach
,
J. B.
Neaton
,
R.
Baer
, and
L.
Kronik
, “
Quasiparticle spectra from a nonempirical optimally tuned range-separated hybrid density functional
,”
Phys. Rev. Lett.
109
,
226405
(
2012
).
30.
K.
Okuno
,
Y.
Shigeta
,
R.
Kishi
,
H.
Miyasaka
, and
M.
Nakano
, “
Tuned CAM-B3LYP functional in the time-dependent density functional theory scheme for excitation energies and properties of diarylethene derivatives
,”
J. Photochem. Photobiol. A: Chem.
235
,
29
34
(
2012
).
31.
D. J.
Tozer
and
N. C.
Handy
, “
Improving virtual Kohn-Sham orbitals and eigenvalues: Application to excitation energies and static polarizabilities
,”
J. Chem. Phys.
109
,
10180
10189
(
1998
).
32.
A.
Dreuw
and
M.
Head-Gordon
, “
Failure of time-dependent density functional theory for long-range charge-transfer excited states: The zincbacteriochlorin-bacteriochlorin and bacteriochlorophyll-spheroidene complexes
,”
J. Am. Chem. Soc.
126
,
4007
4016
(
2004
).
33.
W.
Kohn
and
L. J.
Sham
, “
Self-consistent equations including exchange and correlation effects
,”
Phys. Rev.
140
,
A1133
A1138
(
1965
).
34.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
35.
S.
Kümmel
, “
Charge-transfer excitations: A challenge for time-dependent density functional theory that has been met
,”
Adv. Energy Mater.
7
,
1700440
(
2017
).
36.
S.
Albrecht
,
L.
Reining
,
R.
Del Sole
, and
G.
Onida
, “
Ab initio calculation of excitonic effects in the optical spectra of semiconductors
,”
Phys. Rev. Lett.
80
,
4510
4513
(
1998
).
37.
L. X.
Benedict
,
E. L.
Shirley
, and
R. B.
Bohn
, “
Optical absorption of insulators and the electron-hole interaction: An ab initio calculation
,”
Phys. Rev. Lett.
80
,
4514
4517
(
1998
).
38.
E. L.
Shirley
, “
Ab initio inclusion of electron-hole attraction: Application to x-ray absorption and resonant inelastic x-ray scattering
,”
Phys. Rev. Lett.
80
,
794
797
(
1998
).
39.
C. D.
Spataru
,
S.
Ismail-Beigi
,
L. X.
Benedict
, and
S. G.
Louie
, “
Excitonic effects and optical spectra of single-walled carbon nanotubes
,”
Phys. Rev. Lett.
92
,
077402
(
2004
).
40.
M. L.
Tiago
and
J. R.
Chelikowsky
, “
Optical excitations in organic molecules, clusters, and defects studied by first-principles Green’s function methods
,”
Phys. Rev. B
73
,
205334
(
2006
).
41.
Y.
Ma
,
M.
Rohlfing
, and
C.
Molteni
, “
Excited states of biological chromophores studied using many-body perturbation theory: Effects of resonant-antiresonant coupling and dynamical screening
,”
Phys. Rev. B
80
,
241405
(
2009
).
42.
X.
Blase
,
P.
Boulanger
,
F.
Bruneval
,
M.
Fernandez-Serra
, and
I.
Duchemin
, “
GW and Bethe-Salpeter study of small water clusters
,”
J. Chem. Phys.
144
,
034109
(
2016
).
43.
F.
Bruneval
,
S. M.
Hamed
, and
J. B.
Neaton
, “
A systematic benchmark of the ab initio Bethe-Salpeter equation approach for low-lying optical excitations of small organic molecules
,”
J. Chem. Phys.
142
,
244101
(
2015
).
44.
D.
Jacquemin
,
I.
Duchemin
, and
X.
Blase
, “
Benchmarking the Bethe-Salpeter formalism on a standard organic molecular set
,”
J. Chem. Theory Comput.
11
,
3290
3304
(
2015
).
45.
D.
Jacquemin
,
I.
Duchemin
, and
X.
Blase
, “
Is the Bethe-Salpeter formalism accurate for excitation energies? Comparisons with TD-DFT, CASPT2, and EOM-CCSD
,”
J. Phys. Chem. Lett.
8
,
1524
1529
(
2017
).
46.
K.
Krause
and
W.
Klopper
, “
Implementation of the Bethe-Salpeter equation in the TURBOMOLE program
,”
J. Comput. Chem.
38
,
383
388
(
2017
).
47.
C.
Azarias
,
I.
Duchemin
,
X.
Blase
, and
D.
Jacquemin
, “
Bethe-Salpeter study of cationic dyes: Comparisons with ADC(2) and TD-DFT
,”
J. Chem. Phys.
146
,
034301
(
2017
).
48.
D.
Hirose
,
Y.
Noguchi
, and
O.
Sugino
, “
Quantitative characterization of exciton from GW+Bethe-Salpeter calculation
,”
J. Chem. Phys.
146
,
044303
(
2017
).
49.
J.
Li
,
M.
Holzmann
,
I.
Duchemin
,
X.
Blase
, and
V.
Olevano
, “
Helium atom excitations by the GW and Bethe-Salpeter many-body formalism
,”
Phys. Rev. Lett.
118
,
163001
(
2017
).
50.
T.
Rangel
,
S. M.
Hamed
,
F.
Bruneval
, and
J. B.
Neaton
, “
An assessment of low-lying excitation energies and triplet instabilities of organic molecules with an ab initio Bethe-Salpeter equation approach and the Tamm-Dancoff approximation
,”
J. Chem. Phys.
146
,
194108
(
2017
).
51.
X.
Gui
,
C.
Holzer
, and
W.
Klopper
, “
Accuracy assessment of GW starting points for calculating molecular excitation energies using the Bethe–Salpeter formalism
,”
J. Chem. Theory Comput.
14
,
2127
2136
(
2018
).
52.
W. G.
Aulbur
,
L.
Jönsson
, and
J. W.
Wilkins
, “
Quasiparticle calculations in solids
,”
Solid State Phys.
54
,
1
218
(
2000
).
53.
D.
Golze
,
M.
Dvorak
, and
P.
Rinke
, “
The GW compendium: A practical guide to theoretical photoemission spectroscopy
,”
Front. Chem.
7
,
377
(
2019
).
54.
F.
Aryasetiawan
and
O.
Gunnarsson
, “
Electronic structure of NiO in the GW approximation
,”
Phys. Rev. Lett.
74
,
3221
3224
(
1995
).
55.
S. V.
Faleev
,
M.
van Schilfgaarde
, and
T.
Kotani
, “
All-electron self-consistent GW approximation: Application to Si, MnO, and NiO
,”
Phys. Rev. Lett.
93
,
126406
(
2004
).
56.
F.
Kaplan
,
F.
Weigend
, and
M. J.
van Setten
, “
Quasi-particle self-consistent GW for molecules
,”
J. Chem. Theory Comput.
12
,
2528
2541
(
2016
).
57.
B.
Baumeier
,
M.
Rohlfing
, and
D.
Andrienko
,
J. Chem. Theory Comput.
10
,
3104
(
2014
).
58.
F.
Caruso
,
P.
Rinke
,
X.
Ren
,
M.
Scheffler
, and
A.
Rubio
,
Phys. Rev. B
86
,
081102
(
2012
).
59.
N.
Marom
,
F.
Caruso
,
X.
Ren
,
O. T.
Hofmann
,
T.
Körzdörfer
,
J. R.
Chelikowsky
,
A.
Rubio
,
M.
Scheffler
, and
P.
Rinke
, “
Benchmark of GW methods for azabenzenes
,”
Phys. Rev. B
86
,
245127
(
2012
).
60.
F.
Bruneval
, “
Ionization energy of atoms obtained from GW self-energy or from random phase approximation total energies
,”
J. Chem. Phys.
136
,
194107
(
2012
).
61.
T. A.
Pham
,
H.-V.
Nguyen
,
D.
Rocca
, and
G.
Galli
,
Phys. Rev. B
87
,
155148
(
2013
).
62.
F.
Caruso
,
P.
Rinke
,
X.
Ren
,
A.
Rubio
, and
M.
Scheffler
,
Phys. Rev. B
88
,
075105
(
2013
).
63.
F.
Bruneval
and
M. A. L.
Marques
,
J. Chem. Theory Comput.
9
,
324
(
2013
).
64.
C.
Faber
,
C.
Attaccalite
,
V.
Olevano
,
E.
Runge
, and
X.
Blase
,
Phys. Rev. B
83
,
115123
(
2011
).
65.
M. J.
van Setten
,
F.
Weigend
, and
F.
Evers
,
J. Chem. Theory Comput.
9
,
232
(
2013
).
66.
C.
Faber
,
P.
Boulanger
,
C.
Attaccalite
,
I.
Duchemin
, and
X.
Blase
, “
Excited states properties of organic molecules: From density functional theory to the GW and Bethe-Salpeter Green’s function formalisms
,”
Philos. Trans. R. Soc., A
372
,
20130271
(
2014
).
67.
V.
Atalla
,
M.
Yoon
,
F.
Caruso
,
P.
Rinke
, and
M.
Scheffler
,
Phys. Rev. B
88
,
165122
(
2013
).
68.
X.
Ren
,
P.
Rinke
,
V.
Blum
,
J.
Wieferink
,
A.
Tkatchenko
,
A.
Sanfilippo
,
K.
Reuter
, and
M.
Scheffler
, “
Resolution-of-identity approach to Hartree–Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions
,”
New J. Phys.
14
,
053020
(
2012
).
69.
X.
Blase
and
C.
Attaccalite
, “
Charge-transfer excitations in molecular donor-acceptor complexes within the many-body Bethe-Salpeter approach
,”
Appl. Phys. Lett.
99
,
171909
(
2011
).
70.
D.
Rocca
,
D.
Lu
, and
G.
Galli
, “
Ab initio calculations of optical absorption spectra: Solution of the Bethe–Salpeter equation within density matrix perturbation theory
,”
J. Chem. Phys.
133
,
164109
(
2010
).
71.
B.
Baumeier
,
D.
Andrienko
, and
M.
Rohlfing
, “
Frenkel and charge-transfer excitations in donor–acceptor complexes from many-body Green’s functions theory
,”
J. Chem. Theory Comput.
8
,
2790
2795
(
2012
).
72.
C.
Faber
,
P.
Boulanger
,
I.
Duchemin
,
C.
Attaccalite
, and
X.
Blase
, “
Many-body Green’s function GW and Bethe-Salpeter study of the optical excitations in a paradigmatic model dipeptide
,”
J. Chem. Phys.
139
,
194308
(
2013
).
73.
C.
Faber
,
I.
Duchemin
,
T.
Deutsch
, and
X.
Blase
,
Phys. Rev. B
86
,
155315
(
2012
).
74.
A.
Seidl
,
A.
Görling
,
P.
Vogl
,
J. A.
Majewski
, and
M.
Levy
, “
Generalized Kohn-Sham schemes and the band-gap problem
,”
Phys. Rev. B
53
,
3764
3774
(
1996
).
75.
F.
Bruneval
,
T.
Rangel
,
S. M.
Hamed
,
M.
Shao
,
C.
Yang
, and
J. B.
Neaton
, “
MOLGW 1: Many-body perturbation theory software for atoms, molecules, and clusters
,”
Comput. Phys. Commun.
208
,
149
161
(
2016
).
76.
X.
Blase
,
I.
Duchemin
,
P.
Boulanger
,
C.
Faber
,
C.
Attacalite
,
V.
Olevano
,
J.
Li
, and
G.
D’Avino
, The Fiesta Code (last retrieved October 23, 2019).
77.
F.
Furche
,
R.
Ahlrichs
,
C.
Hättig
,
W.
Klopper
,
M.
Sierka
, and
F.
Weigend
, “
Turbomole
,”
Wiley Interdiscip. Rev. Comput. Mol. Sci.
4
,
91
100
(
2014
).
78.
C.
Holzer
and
W.
Klopper
, “
Ionized, electron-attached, and excited states of molecular systems with spin–orbit coupling: Two-component GW and Bethe–Salpeter implementations
,”
J. Chem. Phys.
150
,
204116
(
2019
).
79.
J.
Deslippe
,
G.
Samsonidze
,
D. A.
Strubbe
,
M.
Jain
,
M. L.
Cohen
, and
S. G.
Louie
, “
BerkeleyGW: A massively parallel computer package for the calculation of the quasiparticle and optical properties of materials and nanostructures
,”
Comput. Phys. Commun.
183
,
1269
1289
(
2012
).
80.
A.
Marini
,
C.
Hogan
,
M.
Grüning
, and
D.
Varsano
, “
yambo: An ab initio tool for excited state calculations
,”
Comput. Phys. Commun.
180
,
1392
1403
(
2009
).
81.
X.
Gonze
,
B.
Amadon
,
P.-M.
Anglade
,
J.-M.
Beuken
,
F.
Bottin
,
P.
Boulanger
,
F.
Bruneval
,
D.
Caliste
,
R.
Caracas
,
M.
Côté
,
T.
Deutsch
,
L.
Genovese
,
P.
Ghosez
,
M.
Giantomassi
,
S.
Goedecker
,
D.
Hamann
,
P.
Hermet
,
F.
Jollet
,
G.
Jomard
,
S.
Leroux
,
M.
Mancini
,
S.
Mazevet
,
M.
Oliveira
,
G.
Onida
,
Y.
Pouillon
,
T.
Rangel
,
G.-M.
Rignanese
,
D.
Sangalli
,
R.
Shaltaf
,
M.
Torrent
,
M.
Verstraete
,
G.
Zerah
, and
J.
Zwanziger
, “
ABINIT: First-principles approach to material and nanosystem properties
,”
Comput. Phys. Commun.
180
,
2582
2615
(
2009
).
82.
X.
Gonze
,
F.
Jollet
,
F. A.
Araujo
,
D.
Adams
,
B.
Amadon
,
T.
Applencourt
,
C.
Audouze
,
J.-M.
Beuken
,
J.
Bieder
,
A.
Bokhanchuk
,
E.
Bousquet
,
F.
Bruneval
,
D.
Caliste
,
M.
Côté
,
F.
Dahm
,
F. D.
Pieve
,
M.
Delaveau
,
M. D.
Gennaro
,
B.
Dorado
,
C.
Espejo
,
G.
Geneste
,
L.
Genovese
,
A.
Gerossier
,
M.
Giantomassi
,
Y.
Gillet
,
D.
Hamann
,
L.
He
,
G.
Jomard
,
J. L.
Janssen
,
S. L.
Roux
,
A.
Levitt
,
A.
Lherbier
,
F.
Liu
,
I.
Lukačević
,
A.
Martin
,
C.
Martins
,
M.
Oliveira
,
S.
Poncé
,
Y.
Pouillon
,
T.
Rangel
,
G.-M.
Rignanese
,
A.
Romero
,
B.
Rousseau
,
O.
Rubel
,
A.
Shukri
,
M.
Stankovski
,
M.
Torrent
,
M. V.
Setten
,
B. V.
Troeye
,
M.
Verstraete
,
D.
Waroquiers
,
J.
Wiktor
,
B.
Xu
,
A.
Zhou
, and
J.
Zwanziger
, “
Recent developments in the ABINIT software package
,”
Comput. Phys. Commun.
205
,
106
131
(
2016
).
83.
G.
Kresse
and
J.
Furthmüller
, “
Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set
,”
Comput. Mater. Sci.
6
,
15
50
(
1996
).
84.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
G. L.
Chiarotti
,
M.
Cococcioni
,
I.
Dabo
,
A. D.
Corso
,
S.
de Gironcoli
,
S.
Fabris
,
G.
Fratesi
,
R.
Gebauer
,
U.
Gerstmann
,
C.
Gougoussis
,
A.
Kokalj
,
M.
Lazzeri
,
L.
Martin-Samos
,
N.
Marzari
,
F.
Mauri
,
R.
Mazzarello
,
S.
Paolini
,
A.
Pasquarello
,
L.
Paulatto
,
C.
Sbraccia
,
S.
Scandolo
,
G.
Sclauzero
,
A. P.
Seitsonen
,
A.
Smogunov
,
P.
Umari
, and
R. M.
Wentzcovitch
, “
QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials
,”
J. Phys.: Condens. Matter.
21
,
395502
(
2009
).
85.
A.
Gulans
,
S.
Kontur
,
C.
Meisenbichler
,
D.
Nabok
,
P.
Pavone
,
S.
Rigamonti
,
S.
Sagmeister
,
U.
Werner
, and
C.
Draxl
, “
Exciting: A full-potential all-electron package implementing density-functional theory and many-body perturbation theory
,”
J. Phys.: Condens. Matter.
26
,
363202
(
2014
).
86.
C.
Vorwerk
,
B.
Aurich
,
C.
Cocchi
, and
C.
Draxl
, “
Bethe–Salpeter equation for absorption and scattering spectroscopy: Implementation in the exciting code
,”
Electron. Struct.
1
,
037001
(
2019
).
87.
K.
Dewhurst
,
S.
Sharma
,
L.
Nordström
,
F.
Cricchio
,
O.
Granäs
,
H.
Gross
,
C.
Ambrosch-Draxl
,
C.
Persson
,
F.
Bultmark
,
C.
Brouder
,
R.
Armiento
,
A.
Chizmeshya
,
P.
Anderson
,
I.
Nekrasov
,
F.
Wagner
,
F.
Kalarasse
,
J.
Spitaler
,
S.
Pittalis
,
N.
Lathiotakis
,
T.
Burnus
,
S.
Sagmeister
,
C.
Meisenbichler
,
S.
Lebègue
,
Y.
Zhang
,
F.
Körmann
,
A.
Baranov
,
A.
Kozhevnikov
,
S.
Suehara
,
F.
Essenberger
,
A.
Sanna
,
T.
McQueen
,
T.
Baldsiefen
,
M.
Blaber
,
A.
Filanovich
,
T.
Björkman
,
M.
Stankovski
,
J.
Goraus
,
M.
Meinert
,
D.
Rohr
,
V.
Nazarov
,
K.
Krieger
,
A.
Davydov
,
F.
Eich
,
A.
Romero Castro
,
K.
Kitahara
,
J.
Glasbrenner
,
K.
Bussmann
,
I.
Mazin
,
M.
Verstraete
,
D.
Ernsting
,
S.
Dugdale
,
P.
Elliott
,
M.
Dulak
,
J. A.
Flores Livas
,
S.
Cottenier
,
Y.
Shinohara
,
M.
Fechner
,
Y.
Kvashnin
,
T.
Müller
,
A.
Gerasimov
,
M. D.
Le
,
J. L.
Bartolomé
,
R.
Wirnata
,
J.
Kumar
, and
A.
Shyichuk
, The ELK code (last retrieved October 23, 2019).
88.
V.
Blum
,
R.
Gehrke
,
F.
Hanke
,
P.
Havu
,
V.
Havu
,
X.
Ren
,
K.
Reuter
, and
M.
Scheffler
, “
Ab initio molecular simulations with numeric atom-centered orbitals
,”
Comput. Phys. Commun.
180
,
2175
2196
(
2009
).
89.
I. Y.
Zhang
,
X.
Ren
,
P.
Rinke
,
V.
Blum
, and
M.
Scheffler
, “
Numeric atom-centered-orbital basis sets with valence-correlation consistency from H to Ar
,”
New J. Phys.
15
,
123033
(
2013
).
90.
V.
Havu
,
V.
Blum
,
P.
Havu
, and
M.
Scheffler
, “
Efficient integration for all-electron electronic structure calculation using numeric basis functions
,”
J. Comput. Phys.
228
,
8367
8379
(
2009
).
91.
A. C.
Ihrig
,
J.
Wieferink
,
I. Y.
Zhang
,
M.
Ropo
,
X.
Ren
,
P.
Rinke
,
M.
Scheffler
, and
V.
Blum
, “
Accurate localized resolution of identity approach for linear-scaling hybrid density functionals and for many-body perturbation theory
,”
New J. Phys.
17
,
093020
(
2015
).
92.
V. W.
Yu
,
F.
Corsetti
,
A.
García
,
W. P.
Huhn
,
M.
Jacquelin
,
W.
Jia
,
B.
Lange
,
L.
Lin
,
J.
Lu
,
W.
Mi
,
A.
Seifitokaldani
,
A.
Vázquez-Mayagoitia
,
C.
Yang
,
H.
Yang
, and
V.
Blum
, “
ELSI: A unified software interface for Kohn–Sham electronic structure solvers
,”
Comput. Phys. Commun.
222
,
267
285
(
2018
).
93.
A.
Marek
,
V.
Blum
,
R.
Johanni
,
V.
Havu
,
B.
Lang
,
T.
Auckenthaler
,
A.
Heinecke
,
H.-J.
Bungartz
, and
H.
Lederer
, “
The ELPA library: Scalable parallel eigenvalue solutions for electronic structure theory and computational science
,”
J. Phys.: Condens. Matter.
26
,
213201
(
2014
).
94.
J. P.
Perdew
and
Y.
Wang
, “
Accurate and simple analytic representation of the electron-gas correlation energy
,”
Phys. Rev. B
45
,
13244
13249
(
1992
).
95.
H. N.
Rojas
,
R. W.
Godby
, and
R. J.
Needs
, “
Space-time method for ab initio calculations of self-energies and dielectric response functions of solids
,”
Phys. Rev. Lett.
74
,
1827
1830
(
1995
).
96.
H. J.
Vidberg
and
J. W.
Serene
, “
Solving the Eliashberg equations by means of N-point Padé approximants
,”
J. Low Temp. Phys.
29
,
179
192
(
1977
).
97.
M. J.
van Setten
,
F.
Caruso
,
S.
Sharifzadeh
,
X.
Ren
,
M.
Scheffler
,
F.
Liu
,
J.
Lischner
,
L.
Lin
,
J. R.
Deslippe
,
S. G.
Louie
,
C.
Yang
,
F.
Weigend
,
J. B.
Neaton
,
F.
Evers
, and
P.
Rinke
, “
GW100: Benchmarking G0W0 for molecular systems
,”
J. Chem. Theory Comput.
11
,
5665
5687
(
2015
).
98.
D.
Golze
,
J.
Wilhelm
,
M. J.
van Setten
, and
P.
Rinke
, “
Core-level binding energies from GW: An efficient full-frequency approach within a localized basis
,”
J. Chem. Theory Comput.
14
,
4856
4869
(
2018
).
99.
D.
Jacquemin
,
I.
Duchemin
, and
X.
Blase
, “
Assessment of the convergence of partially self-consistent BSE/GW calculations
,”
Mol. Phys.
114
,
957
967
(
2016
).
100.
F.
Bechstedt
,
Many-Body Approach to Electronic Excitations: Concepts and Application
(
Springer Berlin Heidelberg
,
2015
).
101.
L. X.
Benedict
,
E. L.
Shirley
, and
R. B.
Bohn
, “
Theory of optical absorption in diamond, Si, Ge, and GaAs
,”
Phys. Rev. B
57
,
R9385
R9387
(
1998
).
102.
D.
Jacquemin
,
I.
Duchemin
,
A.
Blondel
, and
X.
Blase
, “
Assessment of the accuracy of the Bethe-Salpeter (BSE/GW) oscillator strengths
,”
J. Chem. Theory Comput.
12
,
3969
3981
(
2016
).
103.
M. E.
Casida
,
Time-Dependent Density Functional Response Theory of Molecular Systems: Theory, Computational Methods, and Functionals
(
Elsevier
,
Amsterdam
,
1996
).
104.
M. E.
Casida
, “
Time-dependent density functional response theory for molecules
,” in
Recent Advances in Density Functional Methods
, edited by
D.
Chong
(
World Scientific
,
Singapore
,
1995
), pp.
155
192
, https://www.worldscientific.com/doi/pdf/10.1142/9789812830586_0005.
105.
S.
Hirata
and
M.
Head-Gordon
, “
Time-dependent density functional theory within the Tamm-Dancoff approximation
,”
Chem. Phys. Lett.
314
,
291
299
(
1999
).
106.
F.
Furche
and
K.
Burke
, “
Time-dependent density functional theory in quantum chemistry
,” in
Annual Reports in Computational Chemistry
, edited by
D. C.
Spellmeyer
(
Elsevier
,
2005
), Vol. 1, Chap. 2, pp.
19
30
.
107.
J. L.
Whitten
, “
Coulomb potential energy integrals and approximations
,”
J. Chem. Phys.
58
,
4496
(
1973
).
108.
B. I.
Dunlap
,
J. W. D.
Connolly
, and
J. R.
Sabin
, “
On some approximations of Xα method
,”
J. Chem. Phys.
71
,
3396
(
1979
).
109.
M.
Feyereisen
,
G.
Fitzgerald
, and
A.
Komornicki
, “
Use of approximate integrals in ab initio theory, an application in MP2 energy calculations
,”
Chem. Phys. Lett.
208
,
359
(
1993
).
110.
O.
Vahtras
,
J.
Almlöf
, and
M. W.
Feyereisen
, “
Integral approximations for LCAO-SCF calculations
,”
Chem. Phys. Lett.
213
,
514
(
1993
).
111.
F.
Weigend
,
M.
Häser
,
H.
Patzelt
, and
R.
Ahlrichs
, “
RI-MP2: Optimized auxiliary basis sets and demonstration of efficiency
,”
Chem. Phys. Lett.
294
,
143
(
1998
).
112.
A.
Schäfer
,
H.
Horn
, and
R.
Ahlrichs
, “
Fully optimized contracted Gaussian basis sets for atoms Li to Kr
,”
J. Chem. Phys.
97
,
2571
2577
(
1992
).
113.
A. D.
Becke
, “
Density-functional thermochemistry. III. The role of exact exchange
,”
J. Chem. Phys.
98
,
5648
(
1993
).
114.
R. A.
Kendall
,
T. H.
Dunning
, and
R. J.
Harrison
, “
Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions
,”
J. Chem. Phys.
96
,
6796
6806
(
1992
).
115.
M.
Palummo
,
C.
Hogan
,
F.
Sottile
,
P.
Bagala
, and
A.
Rubio
, “
Ab initio electronic and optical spectra of free-base porphyrins: The role of electronic correlation
,”
J. Chem. Phys.
131
,
084102
(
2009
).
116.
T. H.
Dunning
, Jr.
, “
Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen
,”
J. Chem. Phys.
90
,
1007
1023
(
1989
).
117.
S. R.
Jensen
,
S.
Saha
,
J. A.
Flores-Livas
,
W.
Huhn
,
V.
Blum
,
S.
Goedecker
, and
L.
Frediani
, “
The elephant in the room of density functional theory calculations
,”
J. Phys. Chem. Lett.
8
,
1449
1457
(
2017
).
118.
P.
Umari
,
G.
Stenuit
, and
S.
Baroni
, “
Optimal representation of the polarization propagator for large-scale GW calculations
,”
Phys. Rev. B
79
,
201104
(
2009
).
119.
H.
Jiang
and
P.
Blaha
, “
GW with linearized augmented plane waves extended by high-energy local orbitals
,”
Phys. Rev. B
93
,
115203
(
2016
).
120.
K.
Pierloot
, “
The CASPT2 method in inorganic electronic spectroscopy: From ionic transition metal to covalent actinide complexes
,”
Mol. Phys.
101
,
2083
2094
(
2003
).
121.
V.
Veryazov
,
P.
Malmqvist
, and
B. O.
Roos
, “
How to select active space for multiconfigurational quantum chemistry?
,”
Int. J. Quantum Chem.
111
,
3329
3338
(
2011
).
122.
C. J.
Stein
and
M.
Reiher
, “
Automated selection of active orbital spaces
,”
J. Chem. Theory Comput.
12
,
1760
1771
(
2016
).
123.
N. A.
Besley
, “
Calculation of the electronic spectra of molecules in solution and on surfaces
,”
Chem. Phys. Lett.
390
,
124
129
(
2004
).
124.
M. A.
Marques
,
M. J.
Oliveira
, and
T.
Burnus
, “
Libxc: A library of exchange and correlation functionals for density functional theory
,”
Comput. Phys. Commun.
183
,
2272
2281
(
2012
).
125.
S.
Lehtola
,
C.
Steigemann
,
M. J.
Oliveira
, and
M. A.
Marques
, “
Recent developments in libxc—A comprehensive library of functionals for density functional theory
,”
SoftwareX
7
,
1
5
(
2018
).
126.
F.
Bruneval
, “
Optimized virtual orbital subspace for faster GW calculations in localized basis
,”
J. Chem. Phys.
145
,
234110
(
2016
).

Supplementary Material