Reproducibility in $G_0W_0$ Calculations for Solids

TL;DR

This paper investigates the reproducibility of $G_0W_0$ calculations for solids, analyzing convergence, implementation differences, and demonstrating consistent quasiparticle energies across multiple computational codes.

Contribution

It provides the first comprehensive comparison of $G_0W_0$ implementations across different codes, addressing key approximations and convergence issues.

Findings

Consistent quasiparticle energies obtained across BerkeleyGW, Abinit, and Yambo.

Analysis of frequency-integration schemes and Coulomb divergence treatment.

Identification of factors affecting $G_0W_0$ calculation variability.

Abstract

Ab initio many-body perturbation theory within the $GW$ approximation is a Green's function formalism widely used in the calculation of quasiparticle excitation energies of solids. In what has become an increasingly standard approach, Kohn-Sham eigenenergies, generated from a DFT calculation with a strategically-chosen exchange correlation functional ``starting point'', are used to construct $G$ and $W$, and then perturbatively corrected by the resultant $GW$ self-energy. In practice, there are several ways to construct the $GW$ self-energy, and these can lead to variations in predicted quasiparticle energies. For example, for ZnO and TiO$_2$, reported $GW$ fundamental gaps can vary by more than 1 eV. In this work, we address the convergence and key approximations in contemporary $G_0W_0$ calculations, including frequency-integration schemes and the treatment of the Coulomb divergence in the exact-exchange term. We study several systems,and compare three different $GW$ codes: BerkeleyGW, Abinit and Yambo. We demonstrate, for the first time, that the same quasiparticle energies for systems in the condensed phase can be obtained with different codes, and we provide a comprehensive assessment of implementations of the $GW$ approximation.

Figures

Tables

Table 1. Table 1: VBM, CBM and fundamental energy-gap of silicon calculated within G ​ W 𝐺 𝑊 GW with several codes using different frequency-integration schemes. Band energies are shown with respect to the DFT VBM.

	QP energies of silicon (eV)
	PPM-GN				PPM-HL			FF-CD		FF-RA
	ABI	BGW	YMB		ABI	BGW		ABI	BGW	YMB
VBM	-0.64	-0.64	-0.64		-0.95	-0.95		-0.74	-0.79	-0.72
CBM	0.52	0.53	0.52		0.29	0.28		0.48	0.49	0.49
Gap	1.16	1.17	1.16		1.24	1.24		1.22	1.28	1.21

Table 2. Table 2: Absolute G ​ W 𝐺 𝑊 GW self-energy for gold at high-symmetry k points, obtained from a scalar-relativistic PBE DFT calculation. Calculations were performed with three different codes and with the PPM-GN.

	$GW$-PPM self-energy for gold (eV)
	Abinit	BerkeleyGW	Yambo
${\Gamma }_{12}$	-23.33	-23.35	-23.29
$X_5$	-24.25	-24.20	-24.20
$X_{4^{\prime }}$	-12.98	-13.08	-12.97

Table 3. Table 3: QP energies for rutile within a spherical-cutoff technique. This comparison is performed with small convergence parameters: a 6 × 6 × 10 6 6 10 6\times 6\times 10 k-point grid and ϵ cut = 20 subscript italic-ϵ cut 20 \epsilon_{\text{cut}}=20 Ry. The actual QP energies of rutile are shown in Table.  4 . We use different codes and frequency-integration schemes (see text) For PPM-HL, unfulfilled PPM modes ( ω unf. superscript 𝜔 unf. \omega^{\scriptsize\textrm{unf.}} ) are either ∗ * set to 1 1 1  hartree or † † \dagger neglected.

	Rutile: (unconverged) QP energies obtained with a spherical-cutoff method (eV)
	PPM-GN			PPM-HL∗		PPM-HL†		FF-CD		FF-RA
	BGW	Abinit	Yambo	BGW	Abinit	BGW	Abinit	BGW	Abinit	Yambo
VBM	1.66	1.66	1.66	1.53	1.58	1.27	1.32	1.59	1.59	1.59
CBM	5.47	5.47	5.47	5.62	5.58	5.98	5.94	5.45	5.45	5.43
Gap	3.81	3.81	3.81	4.09	4.00	4.71	4.62	3.86	3.86	3.84

Table 4. Table 4: We show the fundamental energy-gap of rutile calculated with G 0 ​ W 0 subscript 𝐺 0 subscript 𝑊 0 G_{0}W_{0} using different set of approximations within different codes, such as the frequency-integration scheme, basis set and pseudopotentials/all-electron. ∗ FF method in the complex plane Zhang  et al. ( 2015 ) .

Rutile TiO2 QP bandgap (eV)
Code	Potential	Freq.	$E_g$	Ref.
Yambo	NC-PP	PPM-GN	3.2	This work
Abinit	NC-PP	PPM-GN	3.2	This work
BGW	NC-PP	PPM-GN	3.2	This work
BGW	NC-PP	FF-CD	3.3	This work
BGW	NC-PP	PPM-HL	3.1	Malashevich et al.,2014
Tombo	AE	PPM-HL	4.0	Zhang et al., 2015, Zhang et al.,2016
Yambo	NC-PP	PPM-GN	3.6	Chiodo et al.,2010
SaX	NC-PP	PPM-GN	3.4	Patrick and Giustino,2012
	AE		4.8	Oshikiri et al.,2003
Yambo	NC-PP	FF-CD	3.3	Kang and Hybertsen,2010
Tombo	AE	FF∗	3.3	Zhang et al., 2015, Zhang et al.,2016

Table 5. Table 5: G ​ W 𝐺 𝑊 GW quasiparticle energies of ZnO within a spherical-cutoff technique. The three G ​ W 𝐺 𝑊 GW codes, Abinit , Yambo  and BGW , agree for the calculated QP energies. This comparison is performed with under-converged parameters: a 5 × 5 × 4 5 5 4 5\times 5\times 4 k-point grid and ϵ cut = 30 subscript italic-ϵ cut 30 \epsilon_{\text{cut}}=30  Ry.

	ZnO QP energies (unconverged) (eV)
	PPM-GN			FF-CD		FF-RA
	BGW	Abinit	Yambo	BGW	Abinit	Yambo
VBM	4.26	4.29	4.26	4.27	4.26	4.26
CBM	8.43	8.43	8.43	8.40	8.42	8.41
Gap	4.17	4.14	4.18	4.14	4.15	4.15

Table 6. Table 6: Fundamental bandgap of ZnO within G 0 ​ W 0 subscript 𝐺 0 subscript 𝑊 0 G_{0}W_{0} @LDA. The converged gap is extrapolated to the CBS, as detailed in the text. The reported bandgaps using different codes and techniques are shown for comparison. * semicore electrons were excluded to calculate the ground-state density required to fit the PPM-HL parameters, see Ref.  Samsonidze  et al. , 2014 . † † \dagger FF integration in the complex plane. + + Frequency integration method based on the random-phase approximation  Anisimov ( 2000 ); Usuda  et al. ( 2002 ) .

ZnO QP bandgap (eV)
Code	Potential	Freq.	$E_g$	Ref.
BGW	NC-PP	PPM-HL	3.4	Shih et al.,2010
Abinit	NC-PP	PPM-HL	3.6	Stankovski et al., 2011, Miglio et al.,2012
Tombo	AE	PPM-HL	4.5	Zhang et al.,2016
Abinit	NC-PP	PPM-HL	2.8	Larson et al.,2013
BGW	NC-PP	PPM-HL∗	3.0	Samsonidze et al.,2014
Abinit	NC-PP	PPM-GN	2.3	Stankovski et al., 2011, Miglio et al.,2012
Abinit	NC-PP	PPM-GN	2.6	Berger et al.,2012
	AE	FF+	2.4	Usuda et al.,2002
Abinit	NC-PP	FF-CD	2.4	Stankovski et al., 2011, Miglio et al.,2012
Vasp	PAW	FF-RA	2.5	Lim et al.,2012
	AE	FF-CD	2.8	Friedrich et al.,2011a; *friedrich_erratum:_2011
Vasp	NC-PAW	FF-RA	2.8	Klimes̆ et al.,2014
Tombo	AE	FF†	2.8	Zhang et al.,2016
BGW	NC-PP	FF-CD	2.8	This work

Table 7. Table 7: Testing norm-conserving pseudopotentials for rutile. For each PP type we show the radii per angular momentum (s, p or d), the plane-wave energy-cutoff ( E cut. subscript 𝐸 cut. E_{\textrm{\scriptsize cut.}} ) (see text), and the corresponding DFT and G ​ W 𝐺 𝑊 GW gap of rutile. We use G 0 ​ W 0 subscript 𝐺 0 subscript 𝑊 0 G_{0}W_{0} PPM-GN with a DFT-PBE starting point, at unconverged G ​ W 𝐺 𝑊 GW parameters (see text).

Pseudopotentials for TiO2
PP	PP radii		$E_{\text{cut.}}$	DFT gap	$GW$ gap
type	(Bohr)		(Ry)	(eV)	(eV)
FHI111 FHI PPs: Ti is defined in the OPIUM-v3.8 user guideopi and O is from the FHI98 library Fuchs and Scheffler (1999); abi .	Ti:	s 1.48, p 1.62, d 1.70	60	1.78	3.12
HGH222HGH PPs from Refs. Hartwigsen et al., 1998; abi, .	Ti:	s 0.34, p 0.24, d 0.24	280	1.88	3.23
	O:	s 0.22, p 0.21
PD333PD PPs from pseudo-dojo-v.2 pse ; van Setten et al. (2018).	Ti:	s 1.35, p 1.30, d 1.65	60	1.88	3.23
	O:	s 1.25, p 1.35

Table 8. Table 8: Sensitiveness of the G ​ W 𝐺 𝑊 GW gap of ZnO with respect to the choice of PPs. Same as TiO 2 in Table  7 . We use G 0 ​ W 0 subscript 𝐺 0 subscript 𝑊 0 G_{0}W_{0} FF-CD with a DFT-LDA starting point at unconverged G ​ W 𝐺 𝑊 GW parameters (see text). Note that the G ​ W 𝐺 𝑊 GW gaps of ZnO shown in this table agree with the converged gap ( = 2.8 absent 2.8 =2.8  eV) due to spurious cancellation of errors.

Pseudopotentials for ZnO
PP	PP radii		$E_{\text{cut.}}$	DFT gap	$GW$ gap
type	(Bohr)		(Ry)	(eV)	(eV)
FHI444FHI PPs from Ref. Stankovski et al., 2011.	Zn:	s 0.80, p 0.80, d 0.80	300	0.67	2.76
	O:	s 1.20, p 1.20
RRKJ555RRKJRappe et al. (1990) PPs from Ref. Samsonidze et al., 2014.	Zn:	s 1.00, p 1.00, d 0.85	300	0.73	2.87
	O:	s 1.10, p 1.10
HGH666HGH PPs from Refs. Hartwigsen et al., 1998; abi, .	Zn:	s 0.40, p 0.53, d 0.25	300	0.73	2.90
	O:	s 0.22, p 0.21
PD777PD PPs from pseudo-dojo-v.2 pse ; van Setten et al. (2018).	Zn:	s 1.35, p 1.65, d 1.85	60	0.78	2.82
	O:	s 1.25, p 1.35
PD888PD PPs generated with the ONCVP code Hamann (2013).	Zn:	s 0.80, p 0.80, d 0.60	500	0.74	2.84
	O:	s 0.80, p 0.80

Equations

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

$$W(\mathbf{r},{\mathbf{r}}^{\prime };\omega )=\int d{\mathbf{r}}^{\prime \prime }{ϵ}^{-1}(\mathbf{r},{\mathbf{r}}^{\prime \prime };\omega )v({\mathbf{r}}^{\prime \prime },{\mathbf{r}}^{\prime }).$$

$$E_i^{\text{QP}}=E_i^{\text{KS}}+⟨{\psi }_i^{\text{KS}}|\Sigma (E_i^{\text{QP}})-V_{xc}|{\psi }_i^{\text{KS}}⟩$$

$$E_i^{\text{QP}}=E_i^{\text{KS}}+Z_i⟨{\psi }_i^{\text{KS}}|\Sigma (E_i^{\text{KS}})-V_{xc}|{\psi }_i^{\text{KS}}⟩,$$

$$Z_i={\left[1-⟨{\psi }_i^{\text{KS}}|{\frac{\partial \Sigma (\omega )}{\partial \omega }|}_{\omega =E_i^{\text{KS}}}|{\psi }_i^{KS}⟩\right]}^{-1}.$$

$$\text{Im}{ϵ}_{\mathbf{G},{\mathbf{G}}^{\prime }}^{-1}(\mathbf{q},\omega )$$

$$\text{Re}{ϵ}_{\mathbf{G},{\mathbf{G}}^{\prime }}^{-1}(\mathbf{q},\omega )$$

$$⟨i\mathbf{k}|{\Sigma }_x|j\mathbf{k}⟩=-\sum _{\mathbf{q},\mathbf{G}}v(\mathbf{q}+\mathbf{G}){\mathcal{F}}_{ij\mathbf{k}}(\mathbf{q}+\mathbf{G})$$

$${\mathcal{F}}_{ij\mathbf{k}}(\mathbf{q}+\mathbf{G})=\sum _{v\in occ.}{M}_{iv\mathbf{k}}(\mathbf{q}+\mathbf{G})M_{jv\mathbf{k}}^{\ast }(\mathbf{q}+\mathbf{G}).$$

$$⟨i\mathbf{k}|{\Sigma }_x|j\mathbf{k}⟩$$

$$f(N_{\text{eps.}},$$

$$\left(\frac{a_1}{N_{\text{eps.}}}+b_1\right)\left(\frac{a_2}{N_{\text{PW}}}+b_2\right)\left(\frac{a_3}{N_{\text{sig.}}}+b_3\right),$$

Full text

I Introduction

Quantitative prediction of charged single-particle excitations in otherwise interacting many-particle systems such as solids is a key component of the design and discovery of materials and the fundamental understanding of matter at the atomistic level. A rigorous formalism for computing such particle-like excitations is many-body perturbation theory, in which electron addition/removal energies are solutions to an effective non-Hermitian single-particle Hamiltonian with a non-local energy-dependent potential, or self-energy operator $\Sigma$. In the so-called $GW$ method, Hedin (1965) the self-energy $\Sigma$ is approximated, to lowest order in the screened Coulomb interaction $W$, as $iGW$, where $G$ is the one-electron Green’s function. In a standard approach, $G$ and $W$ are constructed from a (either regular or generalized Seidl et al. (1996)) Kohn-Sham (KS) eigensystem, computed via density functional theory (DFT), and the KS eigenvalues are corrected perturbatively with a one-shot $G_{0}W_{0}$ self-energy, where the subscript indicates that $G$ and $W$ are not updated self-consistently. By accounting for the screening of the crystal environment, $GW$ is naturally applicable to solids and has proven quite effective in predicting quasiparticle energies of a wide range of crystals.Hybertsen and Louie (1985, 1986); Aryasetiawan and Gunnarsson (1998); Aulbur et al. (2000); Onida et al. (2002); Louie (2006) However, because of the complexity, computational cost, and the number of convergence parameters involved, numerical approximations are required in $GW$ calculations, and varying algorithms in different codes can sometimes yield distinct results.

Crystalline silicon is probably the most-studied test-bed solid for $GW$. Having high crystal symmetry and containing only sp-bonded orbitals, silicon is a relatively-simple system, for which $GW$ within standard approximations yields accurate quasiparticle energies and sizable self-energy corrections.Hybertsen and Louie (1985, 1986) Transition metals (TMs) and transition metal oxides (TMOs), with localized $d$ or $f$ electrons, present a bigger numerical challenge for $GW$. When dealing with TMs, care should be taken in the technical details and approximations used within $GW$. For instance, the convergence criteria,Shih et al. (2010) and the choice of frequency-integration scheme Shaltaf et al. (2008); Kang and Hybertsen (2010); Stankovski et al. (2011); Miglio et al. (2012); Larson et al. (2013) and pseudopotentials Klimes̆ et al. (2014) can yield substantially different results. Several $GW$ works for rutile TiO2 have reported gaps ranging from 3.1 to 4.8 eV,Oshikiri et al. (2003); Kang and Hybertsen (2010); Chiodo et al. (2010); Patrick and Giustino (2012); Malashevich et al. (2014); Zhang et al. (2015) while for ZnO gaps published so far range from 2.6 to 4.5 eV.Shih et al. (2010); Stankovski et al. (2011); Friedrich et al. (2011a); *friedrich_erratum:_2011; Berger et al. (2012); Samsonidze et al. (2014); Klimes̆ et al. (2014); Zhang et al. (2016) Thanks to advances in computational resources and algorithms, recent works have explored convergence beyond past limits,Shih et al. (2010); Stankovski et al. (2011); Friedrich et al. (2011a); *friedrich_erratum:_2011; Berger et al. (2012); Klimes̆ et al. (2014); Zhang et al. (2016) and accurate pseudopotentials specific for $GW$ have been proposed,Klimes̆ et al. (2014); pse a general agreement on the $GW$ quasiparticle energies with different codes has yet to be perceived as being achieved for the difficult cases. The growing popularity of $GW$, the multiple dedicated codes used for $GW$, and the existing challenges and discrepancies encountered when performing $GW$ on increasingly chemically complex systems, such as TMs and TMOs, make it imperative to have reproducibility of predictions from different $GW$ codes.

In this work, we report the results of a detailed comparison of three different plane-wave-based $GW$ codes, and we find that predictions from these codes can agree very well, under given similarly physically sound approximations. For purposes of assessment, we study the representative solids Si, Au, TiO2, and ZnO with the open-source $GW$ codes Abinit (ABI) Gonze et al. (2016), BerkeleyGW (BGW),Deslippe et al. (2012) and Yambo (YMB).Marini et al. (2009) Our benchmark calculations provide a framework for users and developers to document the precision of new applications and methodological improvements, and provides standards for the reproducibility of $GW$ calculations.

II The $GW$ method in practice

The $GW$ method is an interacting Green’s function formalism which accounts for the response of the system to addition or removal of a single electron in an interacting $N$-electrons system, via a non-Hermitian, non-local, and frequency-dependent self-energy operator

[TABLE]

where $\eta$ is a positive infinitesimal and the bare Coulomb potential $v$ and the inverse of the dielectric matrix $\epsilon^{-1}$ are used to construct the screened Coulomb potential

[TABLE]

In the so-called one-shot $GW$, also known as $G_{0}W_{0}$, the quasiparticle energies $E^{\textrm{\scriptsize{QP}}}$ are solved perturbatively from a mean-field Kohn-Sham (KS) starting point; that is, $G_{0}$ and $W_{0}$ are constructed from the KS mean-field. In this approach, which implicitly assumes the KS wavefunctions $\psi^{\textrm{\scriptsize{KS}}}$ are close to the QP wavefunctions $\psi^{\textrm{\scriptsize{QP}}}$, the QP energy of the $i$th state is given byHybertsen and Louie (1985, 1986)

[TABLE]

where $V_{xc}$ is the KS exchange-correlation potential, and $\Sigma$ is evaluated at the QP energy $E^{\textrm{\scriptsize{QP}}}_{i}$. A common approximation is to linearize $\Sigma$ in the QP energy with a first-order Taylor expansion around $E_{i}^{\textrm{\scriptsize{KS}}}$, such that

[TABLE]

with the renormalization factor

[TABLE]

As discussed later, the standard linearlization scheme should be used with care as it can lead to relatively large deviations (up to $0.2$ eV in ZnO) in predicted QP energies.

A source of deviation among $GW$ results with different codes is the numerical integration scheme used to evaluate the frequency depencence of $\Sigma$ in Eq. (1).Miglio et al. (2012); Stankovski et al. (2011); Kang and Hybertsen (2010); Shaltaf et al. (2008) A common practice to reduce computational cost is to approximate the dielectric function with a single-pole via a generalized plasmon-pole model (PPM). For each set of momentum components ($\mathbf{q},\mathbf{G},\mathbf{G^{\prime}}$), the inverse dielectric function $\epsilon^{-1}$ in this approximation takes the form

[TABLE]

where the matrices $A_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q})$ and $\tilde{\omega}_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q})$ are to be determined. Hybertsen and Louie (1986) In the Hybertsen-Louie (HL) approach, the PPM parameters are determined from sum rules and by evaluating the dielectric function at $\omega=0$. Hybertsen and Louie (1986) In the Godby-Needs (GN) scheme, the parameters are set by calculating $\epsilon^{-1}$ at two frequencies: $\omega=0$ and an imaginary frequency close to the plasma frequency.Godby and Needs (1989) Both Abinit and Yambo use the PPM-GN scheme as default; BerkeleyGW uses a PPM-HL version modified to deal with non-centrosymmetric systems.Zhang et al. (1989); Deslippe et al. (2012) When calculating $\epsilon(\mathbf{q},\mathbf{q}^{\prime};\omega=0)$ to find the PPM-HL parameters from Eq. (7), it may happen that the dielectric function cannot be satisfactorily approximated by a single-pole model for certain $(\mathbf{q},\mathbf{G},\mathbf{G}^{\prime})$ leading to imaginary frequencies $\omega_{\mathbf{G},\mathbf{G}^{\prime}}(\mathbf{q})$. Such modes, referred to here as unfulfilled PPM modes $\omega^{\scriptsize\textrm{unf.}}$, are neglected in the original version of the PPM-HL.Hybertsen and Louie (1986) Other treatments of the unfulfilled modes are also possible. For example, these frequencies can be given an arbitrary value of $\omega^{\scriptsize\textrm{unf.}}=1$ Ha, which was the default behavior in Abinit and Yambo.

Beyond PPMs, it is increasingly standard for $GW$ codes to use full-frequency (FF) methods, in which the frequency convolution in Eq. (1) is evaluated numerically. A straightforward integration method on the real axis (FF-RA) is available in codes such as Yambo and BerkeleyGW. However, such an integration of $\Sigma$ in Eq. (1) presents numerical challenges since $G$ and $W$ possess poles close to the real axis. To avoid this difficulty, in the full-frequency contour-deformation (FF-CD) method, the integration contour in Eq. (1) is deformed into the complex plane, into a region where the integrand is smooth; the alternative integration path must be supplemented with the residues from the poles of $G$, as explained in detail in Refs. Lundqvist, 1968; Aryasetiawan and Gunnarsson, 1998; Giantomassi et al., 2011. The FF-CD method is available in Abinit and has been recently implemented into BerkeleyGW.Del Ben et al. (2019); Del Ben et al. (2019) For other FF methods we refer the reader to Refs. Shishkin and Kresse, 2006; Liu et al., 2015; Zhang et al., 2015.

The self-energy is usually split into a frequency-independent exchange part $\Sigma_{x}$ and a correlation part $\Sigma_{c}$, so that $\Sigma(\mathbf{r},\mathbf{r^{\prime}};\omega)=\Sigma_{x}(\mathbf{r},\mathbf{r^{\prime}})+\Sigma_{c}(\mathbf{r},\mathbf{r^{\prime}};\omega)$,Bruneval (2005) where the matrix element of $\Sigma_{x}$ between two Bloch states reads:

[TABLE]

and

[TABLE]

Here, $M_{iv\mathbf{k}}=\langle i\mathbf{k}|e^{i(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}}|v\mathbf{k-q}\rangle$ are matrix elements for states $i$ and $v$ at k-point $\mathbf{k}$. The expression for $\Sigma_{c}$ is given in Ref. Stankovski et al., 2011.

The exchange term, also present in the evaluation of Fock exchange for hybrid functionals in DFT, features a divergence in the Coulomb potential $v(\mathbf{q}+\mathbf{G})=4\pi e^{2}/|\mathbf{q}+\mathbf{G}|^{2}$ as $\mathbf{q}\rightarrow 0$ for $\mathbf{G}=0$. Several schemes have been proposed to treat the divergence of the Coulomb term. Hybertsen and Louie (1987); Baroni and Resta (1986); Gygi and Baldereschi (1986); Pick et al. (1970); Massidda et al. (1993); Carrier et al. (2007); Spencer and Alavi (2008); Rozzi et al. (2006); Ismail-Beigi (2006); Marini et al. (2009); Deslippe et al. (2012) For instance, in the spherical-cutoff technique, the Coulomb interaction is attenuated beyond $R_{c}$ and $v(0)$ is replaced with $2\pi e^{2}R_{c}^{2}$, where the sphere of radius $R_{c}$ has volume equal to that of the unit cell times the number of k-points. Spencer and Alavi (2008) In Abinit by default the Coulomb singularity is approached by an auxiliary-function integration method detailed in Ref. Carrier et al., 2007. Other codes avoid the Coulomb singularity by replacing the value of $q\rightarrow 0$ in Eq. (8) by an integral around $q\simeq 0$.Pulci et al. (1998); Marini et al. (2009); Deslippe et al. (2012) This method is applicable to any $\mathbf{q}$ point in the BZ by assuming,

[TABLE]

where the integral is performed over the BZ region $R_{\mathbf{q+G}}$, which is associated with a volume $\Omega(R_{\mathbf{q+G}})$, and centered around each $\mathbf{q}+\mathbf{G}$ point. This method gives the effect of a larger sampling of points around $\mathbf{q}$ assuming that $\mathcal{F}(\mathbf{q+G})$ is constant over that region.

In the “random integration method” (RIM) implemented in Yambo Marini et al. (2009) and “Monte Carlo averaging” (MC average) technique used in BerkeleyGW Deslippe et al. (2012) the integral is evaluated using a stochastic scheme. In both codes a stochastic scheme is also used to evaluate every term of the form $\int d^{n}\mathbf{q}f(\mathbf{q})v(\mathbf{q})$ in $\Sigma_{c}$, as the scheme can straightforwardly account for integration of arbitrary potentials in regions $R_{\mathbf{q+G}}$ with arbitrary boundaries. Moreover, with the MC averaging scheme, the analytical behavior of $W(\mathbf{q}\rightarrow 0)$ is also appropriately adjusted depending on whether the system behaves like a metal, semiconductor, or displays a graphene-like linearly vanishing density of states; it is also adjusted based on the dimensionality of the system, as discussed in Ref. Deslippe et al., 2012. These stochastic integration methods have shown success in accurately computing the Coulomb singularity and in improving the convergence of $\Sigma$ with respect to k-point sampling.Marini et al. (2009); Deslippe et al. (2012) To facilitate a complete comparison, we also implemented the MC averaging method into Abinit for the present work, as will be discussed below.

Aside from the physical model employed for the dielectric matrix and the treatment of the Coulomb divergence, we emphasize that several parameters must be converged in order to achieve meaningful $GW$ results. Both the calculation of $\epsilon$ and $\Sigma_{c}$ involve unrestricted sums over bands that are truncated up to $N_{\scriptsize\textrm{eps.}}$ and $N_{\scriptsize\textrm{sig.}}$, respectively. Additionally, the codes discussed here use plane-wave basis sets; the number of plane-wave basis functions, $N_{\scriptsize\textrm{PW}}$, used to evaluate $\epsilon$ and $\Sigma$, is expanded up to an energy-cutoff $\epsilon_{\scriptsize\textrm{cut}}$. These three parameters $N_{\scriptsize\textrm{eps.}}$, $N_{\scriptsize\textrm{sig.}}$, and $N_{\scriptsize\textrm{PW}}$ are interdependent, and their convergence needs to be addressed simultaneously.Shih et al. (2010); Stankovski et al. (2011) Here, we extrapolate the $GW$ QP gaps (energy eigenvalue differences) to the complete basis set (CBS) limit with a function of the formDel Ben et al. (2019)

[TABLE]

where $a_{1}$, $a_{2}$, $a_{3}$, $b_{1}$, $b_{2}$, and $b_{3}$ are constants to be determined. Other important convergence parameters and considerations include the k-point sampling of the Brillouin zone, pseudopotential choice, basis used to describe the wavefunctions, and in the case of full-frequency calculations, the frequency sampling on the real and imaginary axis.

III Technical details

In what follows, we compare $GW$ calculations for several materials using three codes implementing the same approaches. For all materials considered, we fix the lattice parameters to the experimental values. These are, for Si in the diamond structure, fcc Au, rutile TiO2, and wurtzite ZnO, respectively, 5.43 Å, 4.08 Å, ($a=4.60$, $c=2.9$) Å, and ($a=3.25$, $c=5.20$) Å. We use norm-conserving Fritz-Haber Institute pseudopotentials with 6, 4, 12 and 20 valence electrons for O, Si, Ti and Zn, respectively. For Au, we use Optimized Norm-Conversving Vanderbilt Pseudoptentials (ONCVP) Hamann (2013) with 19 valence electrons. We use a Perdew-Burke-Ernzerhof (PBE) Perdew et al. (1996) starting point for $GW$, except for ZnO in which the Local Density Approximation (LDA) is used for the sake of comparison to previous works. Our DFT calculations use a k-point mesh and a plane-wave energy-cutoff which ensure that the total energies are converged within 50 meV per unit cell. The k-point mesh is consistent with that for $GW$ calculations, see below; we use a plane-wave energy cutoff to represent wavefunctions of 40, 88, 300 and 300 Ry for silicon, gold, TiO2 and ZnO, respectively. The $GW$ parameters are carefully set to converge quasiparticle energies to 0.1 eV; for silicon, we use a $\Gamma$-centered Monhorst-Pack grid of $12\times 12\times 12$ k-points, $\epsilon_{\scriptsize\textrm{cut}}=20$ Ry and 300 unoccupied states; for gold, we use a mesh of $16\times 16\times 16$ k-points, $\epsilon_{\scriptsize\textrm{cut}}=32$ Ry, and 400 unoccupied states; for rutile TiO2, we use a shifted k-grid of $6\times 6\times 10$ k-points and the number of unoccupied states and $\epsilon_{\scriptsize\textrm{cut}}$ value were extrapolated to the CBS, as detailed in the supplemental materials (SI); and for wurtzite ZnO, we use a shifted k-grid of $8\times 8\times 5$ k-points, and the unoccupied states and $\epsilon_{\scriptsize\textrm{cut}}$ are also extrapolated to the CBS. We summarize in Table I of the Supplementary Information (SI) all convergence parameters used for tables and figures in this manuscript.

IV Results and discussion

IV.1 Silicon

We calculate the $GW$ quasiparticle corrections to the bandstructure of bulk silicon, a typical system for $GW$ calculations. We use a common pseudopotential for all $GW$ calculations, as defined in Section III. The effect of the pseudopotential approximation for silicon is discussed in Ref. Gómez-Abal et al., 2008.

We first study the accuracy of common approximations to treat the Coulomb divergence, which influences the rate of convergence with respect to k-points. In Fig. 1, we show the convergence of the matrix elements of $\Sigma_{x}$ for the valence band maximum (VBM) and conduction band minimum (CBM) at $\Gamma$. We consider different techniques to treat the Coulomb singularity, in particular the MC average in BerkeleyGW for only $\mathbf{q=G}=0$ (black lines, default up to version 1.1 of BerkeleyGW) and for all $\mathbf{G}$ vectors and q-points in the BZ (blue lines, default starting from version 1.2); the RIM for $q=0$ only (brown lines) and all BZ (orange line) in Yambo   and the auxiliary-function treatment Carrier et al. (2007) in Abinit (pink lines). As expected, both the convergence rate with respect to k-points and the converged number of k-points can differ with the choice of method to treat the Coulomb singularity. In this case the RIM and MC average approaches converge fastest, with a grid of $8\times 8\times 8$ k-points being sufficient to converge the $\Sigma_{x}$ matrix elements for the VBM and CBM within 0.05 eV.

In Table 1, we show converged $G_{0}W_{0}$@PBE QP energies for bulk silicon using two different frequency integration schemes and different $GW$ codes. In fact, we find the same QP energies within 0.05 eV for all codes considered here. With respect to the frequency-integration schemes, we find that the PPM in the GN or HL fashions provide a gap for Si within 0.1 eV with respect to the full frequency (FF-CD) reference. Importantly, for a given frequency-integration scheme, the QP energies obtained with the different codes considered here agree within a tolerance better than 0.05 eV, demonstrating that the same $GW$ corrections can be found with different codes.

We highlight that the VBM, CBM, and gap energies calculated with BerkeleyGW and Abinit with FF-CD agree with the energies obtained with Yambo and FF-RA. This result serves as a numerical verification of the equivalence between the implemented FF-CD and FF-RA integration schemes, which was demonstrated exactly only for the electron gas. Lundqvist (1968)

IV.2 Gold

We now revisit the $G_{0}W_{0}$ corrections to the scalar-relativistic band structure of bulk gold, a relatively difficult case for $GW$ due to convergence issues, the non-negligible influence of semicore orbitals on the band structure, and relativistic effects. Rangel et al. (2012); Bernardi et al. (2015) In what follows, we neglect spin-orbit interactions. We first converge the number of bands and $\epsilon_{\scriptsize\textrm{cut}}$, as detailed in the SI; 400 unoccupied states and $\epsilon_{\scriptsize\textrm{cut}}=32$ Ry ensures a convergence of 0.15 eV in the QP gaps between occupied and unoccupied bands across the Brillouin zone in a relatively large window of energies up to $\sim 15$ eV above the Fermi level. Secondly, we uniformly increase the k-point mesh up to $16\times 16\times 16$. We observe differences in $k$-point convergence rate that can be traced to the specific numerical methods used. BerkeleyGW uses a zero-temperature formalism, and a long wavelength limit of the head ($\mathbf{G=G^{\prime}}=0$ component) of the inverse dielectric matrix is $\epsilon^{-1}_{00}(q\rightarrow 0)\sim q^{2}$ specific to metals. This in turn modifies the MC averaging scheme, since the head of the screened Coulomb potential $W_{00}(\mathbf{q})$ is now a finite and smooth function for $q\rightarrow 0$ Deslippe et al. (2012). On the other hand, Abinit and Yambo use finite-temperature occupation factors, requiring a smearing parameter. Here we use Gaussian smearing with a broadening of 0.010 Ry.

In Fig. 2, we show the matrix elements of $\Sigma$ calculated with sets of k-points of increasing size; here we set $\epsilon_{\scriptsize\textrm{cut}}=32$ Ry and $N=400$. As mentioned, the rate of convergence depends on the treatment of band occupations and the Coulomb singularity. While Abinit and Yambo use partial occupations consistent with the underlying DFT code, BerkeleyGW uses a zero-temperature scheme where the bands are either fully-occupied or fully-empty. Moreover, BerkeleyGW uses a particular metal-screening scheme to treat $\epsilon(q\rightarrow 0)$ as described in Ref. Deslippe et al., 2012. With these different approaches, as expected, the self-energy can converge at different rates with respect to the k-point sampling (see Fig. 2). Importantly, when using a relatively-dense mesh of $16\times 16\times 16$ k-points, the codes considered here agree within 0.1 eV in the predicted self-energy of the VBM/CBM at $\Gamma$, demonstrating that for metals the codes predict the same QP energies when convergence is reached.

In Table 2 we show that the matrix elements of $\Sigma$ for bands around the Fermi level calculated with the different codes. The scalar-relativistic DFT band structure and the Brillouin zone are shown in the SI. The $GW$ corrections agree within 0.05 eV, corroborating that at convergence different codes give the same QP energies.

IV.3 Rutile TiO2

Rutile has been the subject of several $GW$ studies, and the reported $G_{0}W_{0}$ gaps range from 3.1 to 4.8 eV Oshikiri et al. (2003); Chiodo et al. (2010); Kang and Hybertsen (2010); Patrick and Giustino (2012); Malashevich et al. (2014); Zhang et al. (2015, 2016). Part of the reported disagreement comes from the treatment of the frequency dependence of $\Sigma$. As detailed in Ref. Kang and Hybertsen, 2010, the fundamental gap calculated with certain PPMs can deviate considerably (by up to 1.1 eV) from a full-frequency reference. The sensitivity of the TiO2 gap to the manner in which the frequency dependence of $\Sigma$ is treated makes rutile an interesting case to investigate the effect and accuracy of PPM and FF methods. As mentioned previously, we use FHI-type pseudopotentials including semicore states consistently in all calculations performed with different codes. Although the choice of pseudopotentials for $GW$ is not studied in this work, we found that our results for rutile are somewhat modified (by less than 0.1 eV) relative to those obtained with other PPs, such as GaussianHartwigsen et al. (1998) and pseudo-dojo-v0.2pse PPs (see Appendix A for more details).

We first examine the $G_{0}W_{0}$@PBE QP energies of rutile TiO2 obtained from different codes, frequency-integration schemes, and in the case of PPMs, choices for $\omega^{\scriptsize\textrm{unf.}}$, as shown in Table 3. The PPM-GN predicts the VBM, CBM, and gap of rutile within $0.1$ eV of the FF reference. The accurate performance of the PPM-GN has been observed consistently for other systems, including other transition metal oxides Shaltaf et al. (2008); Stankovski et al. (2011); Miglio et al. (2012).

We now examine the PPM-HL and in particular the effect of the different choices for $\omega^{\scriptsize\textrm{unf.}}$. Interestingly, when the terms with unfulfilled PPM modes are set to $1$-Ha, the PPM-HL yields results within $0.1$ eV of the PPM-GN and FF approaches, and when neglecting components with $\omega^{\scriptsize\textrm{unf.}}$ the results tend to deviate by up to $0.8$ eV from the FF reference. This clearly indicates that the performance of PPMs for rutile is highly sensitive to the treatment of unfulfilled PPM modes. For rutile, $\omega^{\scriptsize\textrm{unf.}}$ make up an alarming proportion of the dielectric function ($\sim$54% of the matrix elements), which suggests the need for a full-frequency treatment of $\epsilon$, in agreement with Ref. Kang and Hybertsen, 2010. The fraction of unfulfilled PPM modes is therefore an important indicator of whether a full frequency approach is required.

We now compare the $G_{0}W_{0}$ self-energy calculated with different codes in Table 3. When using the PPM-HL, the self-energy can deviate by up to $0.1$ eV for the different codes used here, due to different variants of the PPM-HL being implemented; while Abinit implements the original version of PPM-HL in Ref. Hybertsen and Louie, 1986, BerkeleyGW uses a modified version of the PPM to deal with non-centrosymmetric systems as detailed in Ref. Deslippe et al., 2012. Assessing these small variations in the PPM is beyond the scope of this work. When using the PPM-GN or FF methods, the agreement is better than 20 meV, similar to the silicon case. Importantly, we find that the quasiparticle energies predicted by the different codes agree within $0.1$ eV when using the same treatment of the frequency-dependence.

To converge the $GW$ gap of rutile we extrapolate the interdependent $GW$ parameters ($\epsilon_{\scriptsize\textrm{cut}}$, $N_{\text{sig.}}$ and $N_{\text{eps.}}$) to the CBS limit, as described above and in the SI. The converged bandgap is 3.3 eV for the different codes used here; this result also agrees with previous full-frequency calculations of Refs. Kang and Hybertsen, 2010; Zhang et al., 2015, 2016, as reported in Table 4.

IV.4 Wurtzite Zinc Oxide

Historically, ZnO has been a challenging and controversial system for $GW$. For ZnO, the $GW$ result is strongly affected by the slow convergence of the $\Sigma_{c}$ term Shih et al. (2010). Convergence issues are further aggravated when using PPMs Stankovski et al. (2011), although these PPM-related issues may be partially remedied as illustrated in Ref. Samsonidze et al., 2014. Here we only show results with FF methods and the PPM-GN (validated against FF referencesStankovski et al. (2011)). For more on the PPM approximation for ZnO, we refer the reader to Refs. Stankovski et al., 2011; Miglio et al., 2012; Larson et al., 2013; Samsonidze et al., 2014. Other discrepancies in the $GW$ gap of ZnO arise from the use of incomplete basis-sets and different pseudopotentials, such as projector-augmented waves Klimes̆ et al. (2014). Due to these issues, the reported $G_{0}W_{0}$@LDA gaps with different approximations and codes range from 2.3 – 4.5 eV (see Table 6).

We start by showing that the different codes used here agree on the gap of ZnO, for a given pseudopotential. Again, although pseudopotential issues are not discussed here, we find that our results are insensitive (within 0.1 eV) to the choice of PPs tested in this work, as discussed in Appendix A. In Table 5, we show underconverged QP energies for ZnO calculated using a spherical-cutoff scheme within $G_{0}W_{0}$@LDA. All ZnO results in Table V are computed at the same number of bands, dielectric matrix cutoffs, and k-point grid for the purposes of comparison. However, these parameters are underconverged.) We use the GN method, FF-CD method with Abinit and BerkeleyGW, and the FF-RA method with Yambo. We set $\epsilon_{\scriptsize\textrm{cut}}=30$ Ry, a Coulomb cutoff radius of 19.7177 Bohr, a plasma frequency of 38.82 eV (for PPM-GN), a $\Gamma$-centered homogeneous grid of $5\times 5\times 4$ k-points and 34 bands, and show that the unconverged $GW$ gap of ZnO calculated with the different codes is consistent within 0.1 eV.

Linearizing the self-energy to the QP energy, especially when using coarse $k$-grids, can be inaccurate. An illustration of the difference between the linearized and graphically-solved QP energies is given in Fig. 3a. For the VBM, the linearized and graphical solutions can differ by $\sim 0.2$ eV; for an unconverged set of $GW$ parameters ($5\times 5\times 4$ $k$ grid, $\epsilon_{\scriptsize\textrm{cut}}=40$ Ry and $N_{\scriptsize\textrm{eps.}}=N_{\scriptsize\textrm{sig.}}=2000$), we find $E_{i}^{\textrm{\scriptsize{QP}}}(\Sigma[E_{i}^{\textrm{\scriptsize{QP}}}])=3.7$ eV and $E_{i}^{\textrm{\scriptsize{QP}}}\left(Z\Sigma[E^{\textrm{\scriptsize{KS}}}]\right)=3.9$ eV, where $E_{i}^{\textrm{\scriptsize{KS}}}=5.3$ eV. In Fig. 3b we show the QP bandgap as a function of the number of bands used to evaluate $\Sigma$. We use shifted grids of $5\times 5\times 4$ and $8\times 8\times 5$ k-point grids, $N_{\scriptsize\textrm{eps.}}=2000$ and $\epsilon_{\scriptsize\textrm{cut}}=40$ Ry. Within the coarser grid the actual (blue dots) and linearized (cyan dots) solutions can disagree by more than $0.1$ eV due to features in $\Sigma(\omega)$, as shown in Fig. 3 (a). These features are smoothed out when using a finer grid, reducing the discrepancy associated with linearization.

Having demonstrated good agreement between different codes for ZnO QP energies, we then proceed to converge the gap of ZnO only with BerkeleyGW, excluding the other codes due to our limits on computational resources. To accelerate the convergence with respect to k-points, we use a shifted grid, a common practice well-documented in the pastSottile et al. (2003). Using the finest grid of k-points (that is the $8\times 8\times 5$ grid), we proceed to converge the $N_{\scriptsize\textrm{eps.}}$, $N_{\scriptsize\textrm{sig.}}$ and $\epsilon_{\scriptsize\textrm{cut}}$ by extrapolating to the CBS limit (see SI). As shown in Fig. 4, the bandgap converges linearly with respect to $N_{\scriptsize\textrm{PW}}^{-1}$ and a relatively high $\epsilon_{\scriptsize\textrm{cut}}>80$ Ry is needed to assure convergence within 0.05 eV. At convergence, we find the $G_{0}W_{0}$@LDA gap of ZnO is $2.8$ eV, in agreement with recent calculations, as shown in Table 6.

Finally, we compare our $G_{0}W_{0}$ bandgaps with the corresponding electronic gaps measured in photoemission experiments. Here we use full-frequency $G_{0}W_{0}$ approaches (FF-CD or FF-RA). Note that when comparing to experiment the lattice-renormalization effect should also be included Botti and Marques (2013); Antonius et al. (2014), e.g., the measured/calculated zero-point renormalization (ZPR) of silicon is 62–64 meV, 150 meV for TiO2 and 156–164 meV for ZnO Cardona and Thewalt (2005); Monserrat (2016). Our calculated indirect gap of 1.21–1.28 eV for silicon (without renormalization) is therefore in good agreement with the experimental gap of 1.17 eV Kittel (1986). Our result is also in agreement with the seminal work of Ref. Hybertsen and Louie, 1986. As mentioned above, since we neglect spin-orbit effects in this work, we do not compare the $GW$ bandstructure of gold to experiment. Our calculated gap of 3.3 eV of rutile TiO2 is also in good agreement with the experimental gap of $3.3\pm 0.5$ eV Tezuka et al. (1994); See and Bartynski (1994). On the other hand, our $GW$ gap of ZnO of 2.8 eV substantially underestimates the reported experimental gap of $\sim 3.6$ eV Tsoi et al. (2006); Alawadhi et al. (2007). This well-known shortcoming of standard $G_{0}W_{0}$ for ZnO is due to a deficient LDA starting point Lim et al. (2012), and indicates the need for a more accurate starting point or self-consistent schemes. This work reaches a consensus on the value of the $G_{0}W_{0}$ band-gaps of prototype systems, and hence facilitates future work studying beyond-standard $GW$ schemes to improve the accuracy of $GW$ when using a poor mean-field starting point.

V Conclusions

In this work, we have revisited the $GW$ approximation for prototype systems with three representative plane-wave-based codes: Yambo, Abinit and BerkeleyGW. Within certain choices of approximations and a given set pseudopotentials, the converged $GW$ QP energies calculated with the different codes agree within 0.1 eV, addressing long-standing controversies on the $GW$ results for difficult systems such as ZnO and rutile.

More specifically, we have studied the validity of approximations within one-shot $G_{0}W_{0}$ which can give rise to disagreement in $GW$ results between different codes, in particular the treatment of the Coulomb divergence, convergence, plasmon-pole model approximations, and scheme for capturing the full frequency dependence of $\Sigma$. We have benchmarked different techniques to treat the Coulomb divergence, and identified several effective techniques, in particular an auxiliary-function method used in Abinit, the RIM in Yambo and the MC average in BerkeleyGW. The latter was implemented in Abinit in this work. We have provided new insights into the details of PPMs and their effect on $GW$ results, such as the treatment of unfulfilled PPM modes, which for some systems can lead to large deviations ($>0.5$ eV) from FF references. We have shown that specific PPM techniques, when treated at the same level in the different codes, provide results in complete agreement independently of the code. Beyond the PPM approximation we have also shown that the FF-CD method implemented in BerkeleyGW provides results in agreement with FF implementations in Abinit and Yambo. We highlight that QP energies predicted with the FF-CD method (in the complex plane) agree quantitatively with real-axis FF references, a numerical proof of the validity of the FF-CD.

In summary, our work provides a framework for users and developers to validate and document the precision of new applications and methodological improvements relating to $GW$ codes.

VI Acknowledgments

This work was supported by the Center for Computational Study of Excited State Phenomena in Energy Materials at the Lawrence Berkeley National Laboratory, which is funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division under Contract No. DE-AC02-05CH11231, as part of the Computational Materials Sciences Program. This work is also supported by the Molecular Foundry through the U.S. Department of Energy, Office of Basic Energy Sciences under the same contract number. We acknowledge the use of computational resources at the National Energy Research Scientific Computing Center (NERSC). We also acknowledge the use of HPC resources from GENCI-CCRT-TGCC (Grants No. 2014-096018). F.B. acknowledges the Enhanced Eurotalent program and the France Berkeley Fund for supporting his sabbatical leave in UC Berkeley. D.V., A.F. and A.M. acknowledge support from European Union H2020-INFRAEDI-2018-1 programme under grant agreement No. 824143 project “MaX - materials at the exascale”. AM acknowledges support from European Union H2020-INFRAIA-2015-1 programme under grant agreement No. 676598 project “Nanoscience Foundries and Fine Analysis - Europe”. O.K.O. and D.O’.R. acknowledge the support of Trinity College Dublin’s Studentship Award and School of Physics. Their work was supported by TCHPC (Research IT, Trinity College Dublin), where calculations were performed on the Lonsdale cluster funded through grants from Science Foundation Ireland, and on the Kelvin cluster funded through grants from the Irish Higher Education Authority through its PRTLI program. D.V. and A.F. acknowledge PRACE for awarding access to resource Marconi based in Italy at CINECA. M.J.v.S. and G.-M.R. are grateful to F.R.S.-FNRS for financial support through the PDR Grants T.1031.14 (HiT4FiT). They also thank the CÉCI facilities funded by F.R.S.-FNRS (Grant No. 2.5020.1) and Tier-1 supercomputer of the Fédération Wallonie-Bruxelles funded by the Walloon Region (Grant No. 1117545).

Appendix A The choice of pseudopotential for $GW$

In this appendix, we study the variation of the bandgap with respect to the choice of pseudopotential for TiO2 and ZnO. We emphasize that the validation of pseudopotentials for $GW$ requires all-electron references and is beyond the scope of the present manuscript. In Table 7 we show the $G_{0}W_{0}$ direct gap of rutile calculated with different choices of pseudopotentials. We use a DFT-PBE starting point from Abinit and consider norm-conserving PPs of the Fritz Haber Institute (FHI)Fuchs and Scheffler (1999), Optimized Norm-Conserving Vanderbilt (ONCV)Hamann (2013) and Hartwigsen-Goedecker-Hutter (HGH)Hartwigsen et al. (1998) kinds. The configuration of choice for Ti is [Ne]3s23p63d24s2 (including semicore states), and [He]2s22p6 for O. We only use PPs available in the literature (see Table 7). Note that the HGH and PD PPs contain non-local core corrections (NLCC), which are subtracted from $\Sigma$ when calculating the QP energies. In the table, we show the energy cutoff required to converge the DFT total energy per atom to 0.01 eV and the PP radii, which can be taken as a measure of the PP “hardness”. Here we use BerkeleyGW to compute the $G_{0}W_{0}$ direct gap of rutile using a set of under-converge parameters for $GW$: $N_{\scriptsize\textrm{eps.}}=N_{\scriptsize\textrm{sig.}}=2000$, $\epsilon_{\scriptsize\textrm{cut}}=20$ Ry, the MC avg. technique and a $\Gamma$-centered homogeneous grid of $6\times 6\times 10$ k-points. Importantly, the $GW$ gaps correponding to different PP types agree within 0.1 eV, indicating a small dependence of the gap of rutile with the choice of PPs used here.

We now study the sensitiveness of the $GW$ results with respect to the choice of pseudopotential for ZnO. In Table 8 we show the QP gap of ZnO calculated with $G_{0}W_{0}$@LDA using different PPs. The configuration of choice for Zn is [Ne]3s23p63d104s2 (including semicore states), and [He]2s22p6 for O. As in the TiO2 case, some of the HGH and PSP8 PPs considered here contain NLCCs. We also show the minimum kinetic energy cutoff for the plane-wave expansion to converge the DFT gap within $0.05$ eV, and the corresponding DFT-LDA and $GW$ gaps. Here we use BerkeleyGW, the FF-CD method with 20 imaginary frequencies, an uniform sampling of real frequencies spaced by 0.25 eV from 0 to 6 eV, the modified static-reminder method of Ref. Deslippe et al., 2013 and unconverged $GW$ parameters: $\epsilon_{\scriptsize\textrm{cut}}=30$ Ry, $N_{\scriptsize\textrm{sig.}}=N_{\scriptsize\textrm{eps.}}=500$ . For ZnO the $GW$ and DFT gaps change little, by up to 0.14 and 0.1 eV respectively, with the different choices of PPs. Therefore, the results for ZnO and TiO2 presented in this manuscript are negligibly affected by the choice of PPs.