Ab Initio Bethe-Salpeter Equation Approach to Neutral Excitations in Molecules with Numeric Atom-Centered Orbitals

TL;DR

This paper introduces an all-electron GW+BSE implementation using numeric atom-centered orbitals for molecules, achieving high accuracy in excitation energies and demonstrating efficient basis set convergence for small organic molecules.

Contribution

The paper presents a novel all-electron GW+BSE approach with NAO basis sets for molecules, including benchmarks and convergence analysis, enhancing computational efficiency and accuracy.

Findings

Reproduces reference excitation energies with ~1 meV accuracy.

Achieves excellent convergence with valence correlation consistent NAO basis sets.

Demonstrates similar convergence in TDDFT within NAO formalism.

Abstract

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 meV 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.

Tables

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

$G^0{W}^0$ self-energy	Two-pole		16-Parameter Padé
	Singlet	Triplet	Singlet	Triplet
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

Table 2. Table 2: The number of unoccupied states, N u ​ n ​ o ​ c ​ c subscript 𝑁 𝑢 𝑛 𝑜 𝑐 𝑐 N_{unocc} , averaged over all molecules in Thiel’s set, when imposing threshold values E c ​ u ​ t subscript 𝐸 𝑐 𝑢 𝑡 E_{cut} =20, 40, 60, and 80 eV, or no threshold (“-”) in the BSE@ G 0 ​ W 0 superscript 𝐺 0 superscript 𝑊 0 G^{0}W^{0} @PBE calculations. The tier2+aug2 basis set is used.

Threshold (eV)	20	40	60	80	-
$N_{unocc}$	64	95	117	142	292

Table 3. Table 3: Timings of BSE@ G 0 ​ W 0 superscript 𝐺 0 superscript 𝑊 0 G^{0}W^{0} @PBE and LR-TDDFT-LDA@PBE calculations performed for the neutral singlet excitations of benzene, naphathalene and anthracene. The “tier2+aug” basis set is used and a cutoff energy E c ​ u ​ t subscript 𝐸 𝑐 𝑢 𝑡 E_{cut} = 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).

	BSE@$G^0{W}^0$@PBE				LR-TDDFT-LDA@PBE
Molecules	Basis	Build Mat.	Solver	Total	Basis	Build Mat.	Solver	Total
Benzene	270	9	1	280	293	7	1	301
Naphathalene	1200	44	5	1249	1184	39	5	1228
Anthracene	3089	142	16	3247	2907	129	19	3055

Table 4. Table 4: Timings as in Table 3 , but without any energy cutoff applied for unoccupied states.

	BSE@$G^0{W}^0$@PBE				LR-TDDFT-LDA@PBE
Molecules	Basis	Build Mat.	Solver	Total	Basis	Build Mat.	Solver	Total
Benzene	276	31	7	313	292	22	7	321
Naphathalene	1205	314	72	1592	1200	173	71	1443
Anthracene	3175	1531	407	5113	3010	824	405	4239

Table 5. Table 5: Memory (in Megabytes) used to store the BSE or LR-TDDFT matrix in the BSE@ G 0 ​ W 0 superscript 𝐺 0 superscript 𝑊 0 G^{0}W^{0} @PBE or LR-TDDFT-LDA@PBE calculations of neutral singlet excitations of benzene, naphathalene and anthracene. The column labeled by “Full” denotes the memory required if all unoccupied states enter into the BSE or LR-TDDFT matrix. The column labeled by “ E c ​ u ​ t subscript 𝐸 𝑐 𝑢 𝑡 E_{cut} = 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 “ E c ​ u ​ t subscript 𝐸 𝑐 𝑢 𝑡 E_{cut} = 40 eV” column. Other computational details are the same as those in Table 3 and 4 .

	BSE@$G^0{W}^0$@PBE			LR-TDDFT-LDA@PBE
Molecules	Full	$E_{cut}=40$eV	Ratio	Full	$E_{cut}=40$eV	Ratio
Benzene	400	45	8.9	400	49	8.2
Naphathalene	2549	261	9.8	2549	282	9.0
Anthracene	9129	876	10.4	9129	979	9.3

Equations

$${\hat{h}}^{\text{(g)KS}}|{\psi }_l⟩={ϵ}_l|{\psi }_l⟩,$$

$${\hat{h}}^{\text{(g)KS}}={\hat{t}}_{\text{s}}+{\hat{v}}_{\text{ext}}+{\hat{v}}_{\text{H}}+{\hat{v}}_{\text{xc}}.$$

$${\psi }_l=\sum _i{c}_{li}{\phi }_i,$$

$${\phi }_i(\mathbf{r})=\frac{u_i(r)}{r}{Y}_{lm}(\Omega ).$$

$${ϵ}_l^{GW}={ϵ}_l+⟨{\psi }_l|{\Sigma }^{GW}({ϵ}_l^{GW})-{\hat{v}}_{\text{xc}}|{\psi }_l⟩,$$

$${\Sigma }^{GW}(\text{r},{\text{r}}^{\prime },\omega )=\frac{i}{2\pi }\int d{\omega }^{\prime }G(\text{r},{\text{r}}^{\prime },\omega +{\omega }^{\prime })W(\text{r},{\text{r}}^{\prime },{\omega }^{\prime })e^{i{\omega }^{\prime }\eta }.$$

$$G^0(\text{r},{\text{r}}^{\prime },\omega )=\sum _l\frac{{\psi }_l(\text{r}){\psi }_l^{†}({\text{r}}^{\prime })}{\omega -{ϵ}_l-\text{i}\eta \text{sgn}({ϵ}_{\text{F}}-{ϵ}_l)}.$$

$$W^0(\text{r},{\text{r}}^{\prime },{\omega }^{\prime })=\int d{\text{r}}^{\prime \prime }{\epsilon }^{-1}(\text{r},{\text{r}}^{\prime \prime },{\omega }^{\prime })v({\text{r}}^{\prime \prime },{\text{r}}^{\prime }),$$

$${\Sigma }_{ij}(z)\approx \sum _{n=1}^2\frac{a_n}{z+b_n},$$

$${\Sigma }_{ij}(z)\approx \frac{a_0+a_{1}z+...+a_{(N-2)/2}(z{)}^{(N-2)/2}}{1+b_{1}z+...+b_{N/2}(z{)}^{N/2}},$$

$$L(12;1^{\prime }{2}^{\prime })=L_0(12;1^{\prime }{2}^{\prime })+\int d(3456)L_0(14;1^{\prime }3)K(35;46)L(62;52^{\prime }),$$

$$L_0(12;1^{\prime }{2}^{\prime })=G^0(1,2^{\prime })G^0(2,1^{\prime }),$$

$$K(35;46)=\frac{\delta M(3,4)}{\delta G(6,5)},$$

$$M(3,4)=v_{\mathrm{H}}(3)\delta (3,4)+\Sigma (3,4).$$

$$K({\mathbf{r}}_3,{\mathbf{r}}_5,{\mathbf{r}}_4,{\mathbf{r}}_6)=-iv({\mathbf{r}}_3-{\mathbf{r}}_5)\delta ({\mathbf{r}}_3-{\mathbf{r}}_4)\delta ({\mathbf{r}}_5-{\mathbf{r}}_6)+iW({\mathbf{r}}_3,{\mathbf{r}}_4,\omega =0)\delta ({\mathbf{r}}_3-{\mathbf{r}}_5)\delta ({\mathbf{r}}_4-{\mathbf{r}}_6),$$

$$\left[\begin{array}{cc}A& B\\ -B^{†}& -A^{†}\end{array}\right]\left[\begin{array}{c}{X}_s\\ Y_s\end{array}\right]=E_s\left[\begin{array}{c}{X}_s\\ Y_s\end{array}\right].$$

$$\begin{array}{rl}{A}_{ia}^{jb}=& ({ϵ}_a^{GW}-{ϵ}_i^{GW}){\delta }_{ij}{\delta }_{ab}\\ & -{\alpha }^{S/T}(ia|V|jb)+(ij|W(\omega =0)|ab),\end{array}$$

$$B_{ia}^{jb}=-{\alpha }^{S/T}(ia|V|bj)+(ib|W(\omega =0)|aj).$$

$$(ia|V|jb)=\sum _{pqrs}{c}_{ip}^{†}{c}_{aq}{c}_{jr}{c}_{bs}^{†}(pq|rs),$$

$$(ij|W|ab)=\sum _{pqrs}{c}_{ip}^{†}{c}_{jq}{c}_{ar}{c}_{bs}^{†}(pq|W|rs).$$

$$(pq|rs)=\iint \frac{{\phi }_p(\mathbf{r}){\phi }_q(\mathbf{r}){\phi }_r({\mathbf{r}}^{\prime }){\phi }_s({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime },$$

$$A{X}_s=E_s{X}_s.$$

$$f_s=\frac{2}{3}{E}_s\sum _{\mu =x,y,z}(d_{s,\mu }{)}^2.$$

$$d_{s,\mu }=\sqrt{2}\sum _{ia}⟨{\psi }_i|\hat{\mu }|{\psi }_a⟩(X_{s,ia}+Y_{s,ia}).$$

$$\Omega F_s=E_s^2{F}_s.$$

$${\Omega }_{ia,jb}(\omega )={\delta }_{i,j}{\delta }_{a,b}({ϵ}_a-{ϵ}_i{)}^2+2\sqrt{({ϵ}_a-{ϵ}_i)}{K}_{ia,jb}(\omega )\sqrt{({ϵ}_b-{ϵ}_j)},$$

$$K_{ia,jb}(\omega )=\int \int {\psi }_i^{†}(\mathbf{r}){\psi }_a(\mathbf{r})[\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}+{\tilde{f}}_{\mathrm{xc}}[n_0](\mathbf{r},{\mathbf{r}}^{\prime },\omega )]{\psi }_j({\mathbf{r}}^{\prime }){\psi }_b^{†}({\mathbf{r}}^{\prime })d\text{r}d{\text{r}}^{\prime }.$$

$${\tilde{f}}_{\mathrm{xc}}[n_0](\text{r},{\text{r}}^{\prime },\omega )=\int d(t-t^{\prime })e^{i\omega (t-t^{\prime })}\frac{\delta v_{\mathrm{xc}}[n](\text{r},t)}{\delta n({\text{r}}^{\prime },t^{\prime })}{|}_{n_0}.$$

$${\tilde{f}}_{\mathrm{xc}}^A[n_0](\text{r},{\text{r}}^{\prime })=\frac{\delta v_{\mathrm{xc}}[n](\text{r})}{\delta n({\text{r}}^{\prime })}{|}_{n_0}.$$

$${\phi }_p(\mathbf{r}){\phi }_q(\mathbf{r})\approx \sum _{\mu }{C}_{pq}^{\mu }{P}_{\mu }(\mathbf{r}),$$

Full text

I Introduction

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)Stanton and Bartlett (1993); Bartlett and Musiał (2007); Sekino and Bartlett (1984) or complete active space second-order perturbation theory (CASPT2),Roos, Taylor, and Siegbahn (1980); Andersson, Malmqvist, and Roos (1992); Finley et al. (1998); Ghigo, Roos, and Malmqvist (2004) the quantum Monte Carlo method,Grossman et al. (2001); Tiago et al. (2008); Lee, DuBois, and Kanai (2014); Blunt et al. (2015) linear-response time-dependent density functional theory (LR-TDDFT)Runge and Gross (1984); Marques and Gross (2004); Casida (2009); Casida et al. (1998) or the Bethe-Salpeter equation (BSE) approach in the context of many-body perturbation theory (MBPT).Hedin (1965); Salpeter and Bethe (1951); Hybertsen and Louie (1986); Rohlfing and Louie (1998, 2000); Strinati (1988); Onida, Reining, and Rubio (2002); Blase, Duchemin, and Jacquemin (2018) 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, Stanton and Bartlett (1993); Bartlett and Musiał (2007); Sekino and Bartlett (1984); Andersson, Malmqvist, and Roos (1992); Finley et al. (1998); Ghigo, Roos, and Malmqvist (2004); Schreiber et al. (2008) 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.Jacquemin et al. (2010); Silva-Junior et al. (2008); Imamura and Nakai (2006); Baer and Neuhauser (2005); Refaely-Abramson et al. (2012); Okuno et al. (2012). However, LR-TDDFT calculations can encounter problems for charge-transfer (CT) excitationsTozer and Handy (1998); Dreuw and Head-Gordon (2004), especially when used with a simple XC functional such as the adiabatic local density approximation (LDA)Kohn and Sham (1965) and generalized gradient approximations (GGAs).Perdew, Burke, and Ernzerhof (1996) In LR-TDDFT, including long-range exact exchange in the XC functional can mitigate this problem.Kümmel (2017); Baer and Neuhauser (2005); Refaely-Abramson et al. (2012); Okuno et al. (2012)

The BSE approach is founded upon MBPT, based on Green’s function ($G$) theory and the idea of using the screened Coulomb interaction $W$.Hedin (1965) The BSE formalism was originally proposed in the field of nuclear physics in 1950s.Salpeter and Bethe (1951) Combined with the $GW$ approximation in MBPT,Hedin (1965) the BSE approach has been shown to successfully approximate the optical spectra of solidsRohlfing and Louie (1998); Albrecht et al. (1998); Benedict, Shirley, and Bohn (1998a); Shirley (1998); Spataru et al. (2004) and later work demonstrates similar applicability to excitations in atoms and molecules.Rohlfing and Louie (2000); Tiago and Chelikowsky (2006); Ma, Rohlfing, and Molteni (2009); Blase et al. (2016); Bruneval, Hamed, and Neaton (2015); Jacquemin, Duchemin, and Blase (2015, 2017); Krause and Klopper (2017); Azarias et al. (2017); Hirose, Noguchi, and Sugino (2017); Li et al. (2017a); Rangel et al. (2017) The $GW$ approachHedin (1965); Aulbur, Jönsson, and Wilkins (2000); Strinati (1988) 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.Aryasetiawan and Gunnarsson (1995); Faleev, van Schilfgaarde, and Kotani (2004); Kaplan, Weigend, and van Setten (2016); Baumeier, Rohlfing, and Andrienko (2014); Caruso et al. (2012); Marom et al. (2012); Bruneval (2012a); Pham et al. (2013); Caruso et al. (2013); Bruneval and Marques (2013); Faber et al. (2011); van Setten, Weigend, and Evers (2013); Faber et al. (2014); Atalla et al. (2013); Ren et al. (2012a) 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 for both solids and molecules, which is not captured by LR-TDDFT formalism without long range exchange component.Bruneval, Hamed, and Neaton (2015) CT excitations, which are problematic for LR-TDDFT especially with LDA and GGA functionals, can be efficiently and accurately predicted by the BSE approach.Blase and Attaccalite (2011); Rocca, Lu, and Galli (2010); Baumeier, Andrienko, and Rohlfing (2012); Hirose, Noguchi, and Sugino (2017) This has been demonstrated for CT excitations including both intermolecular and intramolecular types in systems such as simple dipeptides,Faber et al. (2013); Rocca, Lu, and Galli (2010) and more complex fullerene/polymer aggregates.Baumeier, Andrienko, and Rohlfing (2012); Faber et al. (2012)

Calculations of BSE excitation energies within MBPT usually adopt a three-step procedure: (i) Evaluate the Kohn-Sham (KS) Kohn and Sham (1965) or generalized Kohn-Sham (gKS) Seidl et al. (1996) DFT orbitals. (ii) Apply self-energy corrections at the $G^{0}W^{0}$ level or self-consistent $GW$ level ($G^{0}W^{0}$@DFT or $GW$@DFT; $G^{0}$ stands for the Green’s function of a non-interacting reference system and $W^{0}$ is the screened Coulomb interaction of that reference system).Hedin (1965); Aulbur, Jönsson, and Wilkins (2000); Strinati (1988); Rohlfing and Louie (2000) (iii) Solve the BSE (in practice, an approximate version thereof, see below) based on the $G^{0}W^{0}$ or $GW$ quasiparticle energies and screened and unscreened Coulomb integrals of (g)KS orbitals (BSE@$G^{0}W^{0}$@DFT or BSE@$GW$@DFT).Salpeter and Bethe (1951); Jacquemin, Duchemin, and Blase (2015); Bruneval, Hamed, and Neaton (2015); Krause and Klopper (2017) BSE implementations exist in different computational packages based on different basis functions, e.g., MolGWBruneval et al. (2016) and Turbomole,Furche et al. (2014); Krause and Klopper (2017) which are based on Gaussian-type orbital (GTO) basis sets, or BerkeleyGW,Deslippe et al. (2012) Yambo,Marini et al. (2009) Exciting,Gulans et al. (2014) ABINIT,Gonze et al. (2009, 2016) VASPKresse and Furthmüller (1996) and Quantum Espresso,Giannozzi et al. (2009) which are based on plane waves.

The present work introduces an accurate implementation of the BSE formalism utilizing compact and efficient numeric atom-centered orbital (NAO) basis setsBlum et al. (2009); Zhang et al. (2013) in the context of the all-electron electronic structure code FHI-aims.Blum et al. (2009); Havu et al. (2009); Ren et al. (2012a); Ihrig et al. (2015) 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.Ren et al. (2012a) 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.Ren et al. (2012b); Ihrig et al. (2015) Our current implementation also uses the ELSI infrastructureYu et al. (2018) and the ELPA eigenvalue solverMarek et al. (2014) for parallel eigenvalue solutions.

The paper is organized as follows. In Section II, we introduce the $GW$+BSE formalism in the context of MBPT. In Section III, we discuss the details of our implementation. In Section IV, we demonstrate the numerical correctness of our BSE implementation by comparing the excitation energies computed by FHI-aims and by the MolGW codeBruneval, Hamed, and Neaton (2015); Bruneval et al. (2016) with GTO correlation-consistent basis sets for Thiel’s molecular benchmark set.Schreiber et al. (2008); Silva-Junior et al. (2008) 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)Blum et al. (2009) or valence correlation consistent NAO basis sets (NAO-VCC-nZ)Zhang et al. (2013) 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.Kohn and Sham (1965); Perdew and Wang (1992)

II Methods

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.Kresse and Furthmüller (1996); Deslippe et al. (2012); Gonze et al. (2016); Gulans et al. (2014); Marini et al. (2009); Bruneval, Hamed, and Neaton (2015); Jacquemin, Duchemin, and Blase (2015); Krause and Klopper (2017)

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

[TABLE]

[TABLE]

Equation (1) states the electronic (g)KS single-particle equations for the effective single-particle orbitals $\psi_{l}$ and eigenvalues $\epsilon_{l}(l=1,2,...,N_{orbit})$. Equation (2) details the gKS Hamiltonian $\hat{h}^{\textrm{(g)KS}}$, including the effective single-particle kinetic energy (with relativistic corrections) $\hat{t}_{\textrm{s}}$, the external potential $\hat{v}_{\textrm{ext}}$, the electrostatic or Hartree potential of the electron density $\hat{v}_{\textrm{H}}$ and the exchange-correlation potential $\hat{v}_{\textrm{xc}}$. The underlying orbitals $\{\psi_{l}\}$ are here expanded in a basis of NAOs $\{\varphi_{i},i=1,2,...,N_{basis}\}$,

[TABLE]

where the NAOs are of the formBlum et al. (2009)

[TABLE]

$\mathbf{r}$ is a position vector with respect to the nucleus, $r$ is its modulus, and $\Omega$ the corresponding solid angle. In the FHI-aims code, the $u_{i}$ are numerically tabulated functions, defined as cubic splines in units of a logarithmic grid. $Y_{lm}$ are the real-valued spherical harmonics, and $l,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 non-periodic case, $\psi_{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 energiesRen et al. (2012a)

[TABLE]

where $\epsilon^{GW}_{l}$ is the quasiparticle energy. By convention, the arguments $\epsilon^{GW}_{l}$ used to evaluate the self-energy $\Sigma^{GW}$ on the right-hand side are updated self-consistently until they match the $\epsilon^{GW}_{l}$ values obtained on the left-hand side, even though the function $\Sigma^{GW}(\omega)$ 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 HedinHedin (1965):

[TABLE]

In the single-shot perturbative $GW$ (i.e., $G^{0}W^{0}$) approach, the Green’s function $G$ is approximated by the non-interacting Green’s function $G^{0}$, which is calculated from single-particle orbitals $\psi_{l}$ and orbital energies $\epsilon_{l}$:Ren et al. (2012a)

[TABLE]

$\epsilon_{\textrm{F}}$ is the Fermi energy and $\eta$ is a positive infinitesimal. The screened Coulomb potential $W^{0}$ is calculated from the dielectric function $\varepsilon$ asRen et al. (2012a)

[TABLE]

where the dielectric function $\varepsilon$ is obtained at the RPA level, using DFT results. The $G^{0}W^{0}$ self-energy can be calculated using an exact analytic treatment on the real axis, which is the case in the MolGW package.Bruneval et al. (2016); Bruneval (2012b) We refer to Section III.C of Ref. 90 for the details of this formalism. This treatment is limited to small systems. Instead, two-poleRojas, Godby, and Needs (1995) and PadéVidberg and Serene (1977) approximations are implemented in the FHI-aims code for the evaluation of the self-energy on the real axis.Ren et al. (2012a) Both of these approximations are based on an exact treatment of $G^{0}$, $W^{0}$, and the self-energy $\Sigma^{G^{0}W^{0}}$ on the imaginary frequency axis. $\Sigma^{G^{0}W^{0}}$ 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 ($G^{0}$, $W^{0}$, $\Sigma^{G^{0}W^{0}}$) on the imaginary frequency axis significantly reduces the number of frequency points needed, compared to a full frequency integration along the real frequency axis.van Setten et al. (2015) Specifically, the self-energy is approximated to have the following mathematical form in the complex plane in the two-pole approximation:Rojas, Godby, and Needs (1995)

[TABLE]

where the values of $a_{n}$ and $b_{n}$ depend on the indices $i$ and $j$. In the Padé approximation, the self-energy is expressed asVidberg and Serene (1977)

[TABLE]

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 non-unique 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 coworkersGolze et al. (2018) 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 $G^{0}W^{0}$ input data to the BSE are expected from the CD approach. On the other hand, the analytical continuation of $\Sigma$ according to Eqs. (9) or (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 $G^{0}W^{0}$ eigenvalues scales with the size of the system.Golze et al. (2018)

Here we perform one-shot perturbative $G^{0}W^{0}$ calculations based on a fixed DFT or HF reference. The quasiparticle energy in equation (5) is thus rewritten as $\epsilon^{G^{0}W^{0}}_{l}$. Some studies investigate the effect of iterating the $GW$ equations by updating the eigenvalues in equation (7) by those calculated from equation (5), whereas the wavefunctions $\psi_{l}$ are kept at the DFT level.Jacquemin, Duchemin, and Blase (2015); Li et al. (2017b); Jacquemin, Duchemin, and Blase (2016) This procedure, denoted as eigenvalue self-consistent $GW$, is reported to give better agreement with experimental results and wavefunction-based reference methods compared with single-shot $G^{0}W^{0}$ for some systems.Jacquemin, Duchemin, and Blase (2016); Li et al. (2017b); Jacquemin, Duchemin, and Blase (2015, 2017)

(iii) The BSE is a Dyson-like equation for the two-particle correlation function $L$:Rohlfing and Louie (2000); Strinati (1988); Bruneval, Hamed, and Neaton (2015)

[TABLE]

where the set of variables 1, 2, etc. are short for position, time, and spin $(\mathbf{r}_{1},t_{1},\sigma_{1})$, $(\mathbf{r}_{2},t_{2},\sigma_{2})$, etc. $L(12;1^{\prime}2^{\prime})$ is the electron-hole correlation function which describes the probability amplitude of an electron propagating from $1^{\prime}$ to 2 and a hole propagating from 1 to $2^{\prime}$.Strinati (1988) $L_{0}(12;1^{\prime}2^{\prime})$ represents the correlation function of the non-interacting system as defined below in equation (12). $K(35;46)$, usually called the electron-hole interaction kernel, is the screened interaction between the electron and the hole (including bare exchange). $L_{0}$ and $K$ can be expressed in the following equations:Rohlfing and Louie (2000)

[TABLE]

[TABLE]

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

[TABLE]

By applying equation (14) to equation (13), performing a time-energy Fourier transformation and ignoring the dynamical properties of $W$,Rohlfing and Louie (2000) the BSE kernel can be simplified to:

[TABLE]

where the variables 3, 4, etc. are reduced to $\mathbf{r}_{3},\mathbf{r}_{4}$, etc. $v$ is the bare Coulomb interaction. $W$ is the screened Coulomb interaction, with the frequency-dependence ignored.Rohlfing and Louie (2000) 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 $\omega=0$ value of $W$ enters into 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:Rohlfing and Louie (2000); Bruneval, Hamed, and Neaton (2015); Jacquemin, Duchemin, and Blase (2015)

[TABLE]

Here, $E_{s}$ are the optical excitation eigenvalues and ($X_{s},Y_{s}$) are the eigenvectors. The $X_{s}$ and $Y_{s}$ are expressed in electron-hole space of the unperturbed system with elements $X_{s,ia}$ and $Y_{s,ia}$, i.e., the actual BSE wavefunctions are linear combinations of the product of (g)KS orbitals. The excitation wavefunctions $X_{s},Y_{s}$ can be taken to be real-valued in finite (molecular) systems without an external field. Blocks $A$ and $-A^{\dagger}$ correspond to resonant transitions from occupied to unoccupied orbitals and the antiresonant transitions, respectively.Bechstedt (2015) Blocks $B$ and $-B^{\dagger}$ describe the coupling between blocks $A$ and $-A^{\dagger}$. In the BSE, the diagonal matrices $A(-A^{\dagger})$ and off-diagonal matrices $B(-B^{\dagger})$ are defined as:Rohlfing and Louie (2000); Bruneval, Hamed, and Neaton (2015); Jacquemin, Duchemin, and Blase (2015)

[TABLE]

[TABLE]

The indices $i$ and $j$ denote occupied states and $a$ and $b$ denote unoccupied states. $\epsilon_{a}^{GW}$ and $\epsilon_{i}^{GW}$ are the quasiparticle energies denoted as $\epsilon_{l}^{GW}$ in equation (5). The coefficient $\alpha^{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:Bruneval et al. (2016)

[TABLE]

[TABLE]

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

[TABLE]

and the same convention (Mulliken notation) for $\mathbf{r}$ and $\mathbf{r}^{\prime}$ is used in the notation of the screened Coulomb integrals $(pq|W|rs)$ as well. The neglect of the coupling blocks $B$($-B^{\dagger}$) in Eq. (16) is known as the Tamm-Dancoff approximation (TDA).Benedict, Shirley, and Bohn (1998b, c) In the TDA, which also we compare below, the relevant equation becomes simply

[TABLE]

The oscillator strength $f_{s}$ can be calculated from the eigenvalues and eigenvectors obtained by solving the BSE eigenvalue problem:Jacquemin et al. (2016)

[TABLE]

$d_{s,\mu}$ can be calculated asJacquemin et al. (2016)

[TABLE]

Since we are dealing with finite systems, the dipole operator $\hat{\mu}$ is simply taken to be the position operator, i.e., $\hat{\mu}\equiv(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,Casida (1996) which is formally equivalent to Eq. (16). Here, the LR-TDDFT formalism becomes

[TABLE]

In LR-TDDFT, $\Omega$ is called the Casida matrix, which has the same dimension as $A$ or $B$ in Eq. (16). $E_{s}^{2}$ are squares of the neutral many-body excitation energies and $F_{s}$ are the eigenvectors of this eigenvalue problem, which can also be related to the oscillator strengths via the dipole operator.Casida (1995) The Casida matrix can be written in a basis of products of (g)KS orbitals as

[TABLE]

where $\delta$ denotes the Kronecker delta. The kernel $K_{ia,jb}$ is defined as

[TABLE]

$\tilde{f}_{\mathrm{xc}}[n_{0}]$ is the exchange correlation kernel and $n_{0}$ is the (g)KS ground state electron density. The kernel is formally defined through the functional derivative of the time-dependent Kohn-Sham exchange and correlation potential $v_{\mathrm{xc}}[n](\textbf{r},t)$ with respect to the time-dependent density $n(\textbf{r}^{\prime},t^{\prime})$ such that

[TABLE]

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

[TABLE]

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

III Implementation

In our implementation, the two-electron Coulomb interaction in equations (17), (18) and (26), the static screened Coulomb interaction in equations (17) and (18), as well as the two-electron integrals of the exchange correlation kernel in equation (26), are calculated employing the RI approach.Whitten (1973); Dunlap, Connolly, and Sabin (1979); Feyereisen, Fitzgerald, and Komornicki (1993); Vahtras, Almlöf, and Feyereisen (1993); Weigend et al. (1998); Ren et al. (2012a) The RI represents pair products of atomic basis functions $\varphi_{p}(\mathbf{r})\cdot\varphi_{q}(\mathbf{r})$ in terms of auxiliary basis functions (ABFs)

[TABLE]

where $P_{\mu}(\mathbf{r})(\mu=1,2,...,N_{aux})$ are the ABFs and $C^{\mu}_{pq}$ are the expansion coefficients. The construction of the ABFs in FHI-aims is explained in Ref. 66 and in detail in Ref. 85. The evaluation of the integrals (21) then reduces to

[TABLE]

[TABLE]

The computation of the expansion coefficients $C^{\mu}_{pq}$ requires three-center integrals involving the ABFs and the pair products of the NAOs:

[TABLE]

[TABLE]

$(\nu|\mu)^{-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:

[TABLE]

using

[TABLE]

$(\nu|\mu)^{-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 memory. The screened Coulomb interaction $W$ can be represented in terms of the ABFs in a similar way to the Coulomb interaction $V$:

[TABLE]

The dielectric function $\varepsilon$ can be calculated as

[TABLE]

where $\chi^{0}$ is the non-interacting density response function. In real space and for a non-spinpolarized system, $\chi^{0}$ is defined as

[TABLE]

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

IV Results

IV.1 Numerical Validation

We quantify the precision of our BSE implementation by calculating the neutral excitation energies of the molecular benchmark set published by Schreiber et al.,Schreiber et al. (2008) 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 supporting information 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,Silva-Junior et al. (2008); Bruneval, Hamed, and Neaton (2015); Jacquemin, Duchemin, and Blase (2015) 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. We therefore do not compare to the BE results, but rather 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.Schreiber et al. (2008) (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.Schäfer, Horn, and Ahlrichs (1992) A previous study by Bruneval et al.Bruneval, Hamed, and Neaton (2015) indicates that BSE@$G^{0}W^{0}$@DFT-B3LYPBecke (1993) excitation energies for ethene and pyrrole, obtained with the TZVP basis set, overestimate the analogous results from the much larger aug-cc-pVQZ basis setKendall, Dunning, and Harrison (1992) 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 to MolGW using the TZVP basis set.Schäfer, Horn, and Ahlrichs (1992) We then discuss basis set convergence for low-lying singlet and triplet excitations using NAO basis sets.

In Figure 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=\{D_{i,n}\}$ in comparison to a reference set $R=\{R_{i,n}\}$ is defined as:

[TABLE]

$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 PBEPerdew, Burke, and Ernzerhof (1996) exchange-correlation functional is used for the initial DFT calculations and the $G^{0}W^{0}$ 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 MOLGW package at the level of 1 meV or below.

Next, Figure 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 $\mathrm{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 the Table S1-S5 of the Supplementary Material (SM). Since the calculation results of the FHI-aims and MolGW package agrees within 1 meV for the excitation energy and $10^{-4}$ for the oscillator strength, the values in Table S1-S5 are valid for both packages within the significant digits given. In short, Figures 1 and 2 validate our implementation.

IV.2 Effects of the frequency dependence of the self-energy in GW

We next investigate how the BSE energies are impacted by different treatments of the frequency dependence of the $G^{0}W^{0}$ self-energy $\Sigma^{G^{0}W^{0}}$. In the previous section we took the $G^{0}W^{0}$ quasiparticle energies calculated by MolGW as the input for the FHI-aims BSE calculations. $G^{0}W^{0}$ calculations in MolGW employ an exact analytic treatment of $\Sigma^{G^{0}W^{0}}$ on the real axis.Bruneval et al. (2016); Bruneval (2012b) A similarly precise result is expected from the CD technique by Golze et al.Golze et al. (2018) In the present section, we investigate the impact of two frequently employed inexact but potentially cost-savingGolze et al. (2018) approximations to the self-energy on the real axis, namely the two-pole approximation,Rojas, Godby, and Needs (1995) Eq. (9), and the Padé approximation,Vidberg and Serene (1977) Eq. (10), with 16 parameters.

In Figure 3, we show mean absolute errors MAE($i$) of the $G^{0}W^{0}$ quasiparticle energies of the HOMO-3 to LUMO+3 states, calculated either using the two-pole approximation (Fig. 3(a)) or the 16-parameter Padé approximation (Fig. 3(b)) and compared to the MolGW reference values. The $G^{0}W^{0}$ results are based on DFT calculations using the Perdew-Burke-Ernzerhof (PBE)Perdew, Burke, and Ernzerhof (1996) exchange-correlation functional and employ the Gaussian-type TZVP basis set.Schäfer, Horn, and Ahlrichs (1992) We see that the two-pole approximation gives MAE values of up to 0.3 eV in the investigated states. Although not a small value, this must be viewed in the context of plain DFT errors (no $GW$ approximation), which are typically of the order of several eV. The 16-parameter Padé approximation can reduce this error by a factor of two or better, 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 is closely in line with a broader analysis performed in Ref. 93. However, another practical advantage of the two-pole approximation is its relative numerical simplicity and therefore its relative robustness against numerical errors, compared to the Padé approximation. In practice, the Padé approximation can result in serious ambiguities for individual eigenvalues $\epsilon_{l}^{GW}$ if there are multiple solutions for the self-consistency condition between the left and the right sides of Eq. (5). Thus, the two-pole approximated self-energy can be preferable for simple reasons of stability, at the price of reduced accuracy compared to a formally exact $G^{0}W^{0}$ self-energy.

As shown in Figure 4, the different approximate $G^{0}W^{0}$ self-energy treatments affect the BSE excitation energies, which take the $G^{0}W^{0}$ quasiparticle energies as input. We compare BSE results based on $G^{0}W^{0}$ quasiparticle energies calculated using the self-energies of Eqs. (9) and (10) to BSE results based on the MolGW $G^{0}W^{0}$ eigenvalues with the exact analytic self-energy treatment.Bruneval et al. (2016); Bruneval (2012b) 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$) ($\approx$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 $G^{0}W^{0}$ quasiparticle energies (Figure 3), which constitute the input for the BSE calculations.

In addition to the MAE of the BSE@$G^{0}W^{0}$@PBE results between FHI-aims and MolGW, we can also define the state-dependent mean signed error MSE($i$),

[TABLE]

Just as in Eq. (40), $i$ and $n$ are indices for state and molecules, respectively, averaging over $N$ = 28 molecules in Thiel’s set. Table 1 summarizes the overall MSE and MAE (averaged over states $i$ in addition to averaging over molecules) of the BSE@$G^{0}W^{0}$@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@$G^{0}W^{0}$ excitation energies based on an exact $G^{0}W^{0}$ self-energy.

It is interesting to compare the errors incurred from the different $G^{0}W^{0}$ self-energy treatments to the errors associated with the BSE approach itself. Bruneval et al.Bruneval, Hamed, and Neaton (2015) show that the BSE@$G^{0}W^{0}$@PBE singlet excitation energies at the TZVP basis set level give a MSE of $-$0.8 eV and a MAE of 0.8 eV compared to the BE results of Schreiber *et al.*Schreiber et al. (2008) The MSE and MAE of the BSE@$G^{0}W^{0}$@PBE triplet excitation energies are around $-$1.2 eV and 1.2 eV, respectively.Bruneval, Hamed, and Neaton (2015) In comparison, the error incurred through the approximate $G^{0}W^{0}$ self-energies in Table 1 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 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@$G^{0}W^{0}$@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 $G^{0}W^{0}$ calculations and actual neutral excitation energies. This correction is much larger than the differences incurred above as a result of the approximate $G^{0}W^{0}$ self-energies. In our calculations, the BSE@$G^{0}W^{0}$@PBE results using both 16-parameter Padé approximation and two-pole approximations reduce the lowest $G^{0}W^{0}$ 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.Li et al. (2017a) show a $G^{0}W^{0}$@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@$G^{0}W^{0}$@HF results agree with the exact result to better than $\approx$0.01 Ha ($\approx$0.3 eV). For larger molecules, such as free-base porphyrin and tetraphenylporphyrin,Palummo et al. (2009) BSE@$G^{0}W^{0}$LDA excitations can match experimental absorption spectra to a similarly good degree (essentially exact within the remaining approximations made in Ref. Palummo et al. (2009)). Here, the BSE again corrects the simple $G^{0}W^{0}$ HOMO-LUMO energy difference by several eV, much more than the magnitude of changes due to analytical corrections to the $G^{0}W^{0}$ self-energies reported above.

IV.3 Basis Set Convergence

We now address the basis set convergence of the NAO basis sets for the BSE calculation, in comparison to cc-pVnZ basis setsDunning Jr. (1989) and aug-cc-pVnZ basis setsKendall, Dunning, and Harrison (1992) by Dunning and coworkers. Two types of NAO basis sets have been constructed in the context of FHI-aims. The first, denoted as “FHI-aims-2009”, was introduced in Ref. 82 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-meV/atomJensen et al. (2017) 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, shown in Table S6 in the SM. The second type of basis set is denoted as NAO-VCC-nZ with n = 2, 3, 4, 5.Zhang et al. (2013) The NAO-VCC-nZ basis functions are constructed in analogy to the cc-pVnZ correlation-consistent (cc) basis sets by Dunning,Dunning Jr. (1989) 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.Zhang et al. (2013) In the following, the above two types of NAO basis sets, as well as cc-pVnZ (n = 2, 3, 4, 5)Dunning Jr. (1989) and aug-cc-pVnZ (n = 2, 3, 4, 5) basis sets,Kendall, Dunning, and Harrison (1992) are compared to the aug-cc-pV5Z basis sets as the benchmark reference in the BSE calculations.

In Figure 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 $\Delta E_{i}$ 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 panel (b) and (c):

[TABLE]

The $\Delta E_{i}$ 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 (c), we see that the results obtained with the cc-pVnZ basis sets converge slowly towards 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,Kendall, Dunning, and Harrison (1992) 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:Bruneval, Hamed, and Neaton (2015); Jacquemin, Duchemin, and Blase (2015) E.g., in Fig. 3 of Ref. 43, the LR-TDDFT and BSE@$G^{0}W^{0}$@B3LYPBecke (1993) excitation energies of the Ethene and Pyrrole molecules using the aug-cc-pVDZ and aug-cc-pVQZ basis sets show differences of 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 (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 un-augmented 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.

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 function ($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 DFTJensen et al. (2017) 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.

Figure 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 $\bar{N}_{\textrm{basis}}$ of different basis sets, averaged over all molecules in Thiel’s set, is plotted as the $y$ axis in Fig. 6(a). $\bar{N}_{\textrm{basis}}$ is defined as:

[TABLE]

Here, $N_{\textrm{mol}}$ is the number of molecules in Thiel’s set and $N^{j}_{\textrm{basis}}$ is the basis size for 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@$G^{0}W^{0}$@PBE using the aug-cc-pV5Z and tier2+aug2 basis sets are also listed in Table S7 in the SM for all molecules in Thiel’s set.

IV.4 Convergence with respect to the BSE matrix size

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.Umari, Stenuit, and Baroni (2009) All-electron approaches to $G^{0}W^{0}$ 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.Jiang and Blaha (2016) 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. 112, 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.Roos, Taylor, and Siegbahn (1980); Pierloot (2003); Veryazov, Malmqvist, and Roos (2011); Stein and Reiher (2016) As 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 adsorbant.Besley (2004)

In Figure 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 $E_{cut}$ for unoccupied states, limiting the number of states $a$ and $b$ entering the matrix construction in Eqs. (17) and (18). Here, $E_{cut}$ 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 $E_{cut}$ are tabulated in Table 2. Figure 7 shows MAE($i$) values for the lowest ten singlet (subfigure 7a) or triplet (subfigure 7b) 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 $E_{cut}$=20 eV is about 20-30 meV. Setting $E_{cut}$=40 eV yields excitation energies with an error closer to 10 meV. Larger $E_{cut}$ values of 60 and 80 eV lead to further small improvements. In view of the remaining errors of these calculations (due to level of theory for underlying DFT, quasiparticle energies, 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 number of states needed in construction of the BSE matrix. Specifically, the time and memory complexity of constructing the BSE matrix in Eqs. (17) and (18) scale as  ${N}_{occ}^{2}\times{N}_{unocc}^{2}$, where $N_{occ}$ and $N_{unocc}$ denote the number of occupied and unoccupied (g)KS single-particle states, respectively. By setting $E_{cut}$=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 2). Additionally, the formal effort for solving the full BSE, Eq. (16), would scale as $O(N^{6})$, where $N$ is a measure of system size, due to the same considerations of how $N_{occ}$ and $N_{unocc}$ impact the matrix dimension. While imposing $E_{cut}$ 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 $E_{cut}$, without sacrificing much accuracy.

IV.5 Comparison to LR-TDDFT

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 Section IV.3. The LR-TDDFT in FHI-aims was implemented following Eqs. (25-29). The exchange-correlation kernel used is the LDA, employing the parametrization of the correlation energy by Perdew and Wang,Kohn and Sham (1965); Perdew and Wang (1992) provided by the Libxc library.Marques, Oliveira, and Burnus (2012); Lehtola et al. (2018) 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 $f_{xc}$ (Eq. (29)) 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.Perdew, Burke, and Ernzerhof (1996) We will use the notation “LR-TDDFT-LDA@PBE” in the following to denote this approach.

Fig. 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). Fig. 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 Figure 6. Excitation energy values derived from both the cc-pVnZ basis sets and the NAO-VCC-nZ basis sets converge slowly towards 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 towards 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 set are provided in Tables S8 and S9 in the SM for all molecules in Thiel’s set.

We finally compare the results of BSE@$G^{0}W^{0}$@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,Bruneval, Hamed, and Neaton (2015); Jacquemin, Duchemin, and Blase (2015) 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 PBEPerdew, Burke, and Ernzerhof (1996) exchange-correlation functional. To maintain consistency, the TDA is employed in the BSE calculations, as this is widely done for LR-TDDFT calculations as well. Figure 9 shows the MSE($i$) and the MAE($i$) between LR-TDDFT-LDA@PBE results and the BSE@$G^{0}W^{0}$@PBE results, averaged over all molecules in Thiel’s set. The BSE@$G^{0}W^{0}$@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@$G^{0}W^{0}$@PBE results for singlet and triplet excitations states is plotted in Fig. 9 (a) and (b), respectively. The MAE between LR-TDDFT-LDA@PBE and BSE@$G^{0}W^{0}$@PBE results for singlet and triplet excitations states is plotted in Fig. 9 (c) and (d), respectively. We see from Fig. 9 that the MAE between LR-TDDFT-LDA@PBE and BSE@$G^{0}W^{0}$@PBE for singlet excitation energies is less than 0.5 eV, where as the MAE for triplet excitations energies can be as large as 1.0 eV. Singlet excitation energy values obtained from LR-TDDFT-LDA@PBE tend to be lower than those obtained from BSE@$G^{0}W^{0}$@PBE, whereas LR-TDDFT-LDA@PBE triplet excitation energies appear to be larger by up to 1 eV than those obtained from BSE@$G^{0}W^{0}$@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.Bruneval, Hamed, and Neaton (2015); Jacquemin, Duchemin, and Blase (2015) Bruneval et al.Bruneval, Hamed, and Neaton (2015) compare the MSE of both BSE$G^{0}W^{0}$@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.Bruneval, Hamed, and Neaton (2015); Silva-Junior et al. (2008) However, there are several differences between their comparison and the comparison shown in this work. First, the dataset used by Bruneval and coworkers are the BEs of Thiel’s set,Schreiber et al. (2008) 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.Bruneval, Hamed, and Neaton (2015) Finally, we here use a basis set that is essentially converged for both BSE and LR-TDDFT calculations. Another set of comparable results are therefore those of Jacquemin and coworkers,Jacquemin, Duchemin, and Blase (2015) who compare the BSE@$G^{0}W^{0}$@PBE0 and LR-TDDFT@PBE0 in a benchmark paper using the aug-cc-pVTZ basis, which has similar convergence behavior as the tier2+aug2 basis set used here. Different MAE values between BSE@$G^{0}W^{0}$@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; 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 - 1.0 eV for different states of singlet and triplet excitation energies.

IV.6 Remarks on Time and Memory Cost of the BSE and LR-TDDFT Implementation

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 future work in our own implementation. In the present section, we consider a series of acene molecules of increasing size: benzene, naphathalene and anthracene. Both BSE and LR-TDDFT calculations were performed using the “tier2+aug” basis sets described in Sections IV.3 and IV.5, respectively, and employing the TDA.

In Table 3, we show the timings of BSE@$G^{0}W^{0}$@PBE and LR-TDDFT-LDA@PBE performed for the three acene molecules selected. We apply a cutoff energy $E_{cut}$ = 40 eV to limit the number of unoccupied states entering the BSE and Casida matrix construction as described in Section IV.4. 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 in the TDA 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 Equations 17 and 18, versus the LR-TDDFT matrix in Equation 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 $G^{0}W^{0}$ and DFT-PBE steps in the BSE@$G^{0}W^{0}$@PBE calculation and the DFT-PBE step in the LR-TDDFT-LDA@PBE calculation. Table 3 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 4, 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.

To demonstrate how the cutoff energy $E_{cut}$ = 40 eV reduces the computational timings, we show in Table 4 the timings for the same systems without using any cutoff energy for both BSE@$G^{0}W^{0}$@PBE and LR-TDDFT-LDA@PBE calculations. All other computational details remain the same as used to obtain Table 3. We see in Table 4 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 Table 3. In contrast, the timings for the product basis setup step stays 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 4. Specifically, the matrix building step for LR-TDDFT requires the calculations of the kernel $K_{ia,jb}$ (Equation 26 and 27). In our parallel implementation, the state indices $i,a,j,b$ are distributed among different processors, as are the RI two- and three-center integrals, and the calculation of $K_{ia,jb}$ requires significant inter-processor communication to get the correct state indices. For BSE, the analogous inter-processor 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 $E_{cut}$ = 40 eV reduces the memory requirements for the BSE and LR-TDDFT matrices. In Table 5, 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 $E_{cut}$ = 40 eV, indicating a reduction by a factor of $\sim$9-11. This is consistent with Section IV.4, where $E_{cut}$ = 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  ${N}_{occ}^{2}\times{N}_{unocc}^{2}$ (Eqs. (17) and (18), the overall reduction amounts to $(1/3)^{2}$=$1/9$, which is confirmed by the memory reduction shown in Table 5.

V Conclusions

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 moleculesSchreiber et al. (2008) demonstrate the numerical correctness of the implementation. Mean absolute errors of 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 $G^{0}W^{0}$ self-energy (the two-pole and 16-parameter Padé approximations), which impact the $G^{0}W^{0}$ quasiparticle energies entering the BSE. The MAE of the BSE@$G^{0}W^{0}$@PBE results with these analytical approximations is around 0.05 - 0.20 eV, compared with the exact $G^{0}W^{0}$ 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 the two-pole approximation offers an overall numerically more stable avenue. 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,Dunning Jr. (1989) FHI-aims-2009,Blum et al. (2009) NAO-VCC-nZ,Zhang et al. (2013) and aug-cc-pVnZKendall, Dunning, and Harrison (1992) 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,Kendall, Dunning, and Harrison (1992) called “tier2+aug2”. For both BSE@$G^{0}W^{0}$@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 $E_{cut}$=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@$G^{0}W^{0}$@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@$G^{0}W^{0}$@PBE and LR-TDDFT-LDA@PBE is quantified by a MAE in the range of 0.2 - 1.0 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@$G^{0}W^{0}$@DFT.

Supplementary Material

See 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; basis set convergence of excitation energies of the molecules in Thiel’s set using the BSE and LR-TDDFT.

Acknowledgements.

This work was financially supported by the NSF under Awards No. DMR-1729297 and No. DMR-1728921, as well as through the Research Triangle MRSEC (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 DE-AC02-06CH11357. We 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. XR acknowledges the financial support from Chinese National Science Foundation (Grant No. 11574283, 11874335) and the Max Planck Partner Group project.