GW100: Benchmarking G0W0 for Molecular Systems

: We present the GW 100 set. GW 100 is a benchmark set of the ionization potentials and electron a ﬃ nities of 100 molecules computed with the GW method using three independent GW codes and di ﬀ erent GW methodologies. The quasi-particle energies of the highest-occupied molecular orbitals (HOMO) and lowest-unoccupied molecular orbitals (LUMO) are calculated for the GW 100 set at the G 0 W 0 @PBE level using the software packages TURBOMOLE, FHI-aims, and Berkeley GW . The use of these three codes allows for a quantitative comparison of the type of basis set (plane wave or local orbital) and handling of unoccupied states, the treatment of core and valence electrons (all electron or pseudopotentials), the treatment of the frequency dependence of the self-energy (full frequency or more approximate plasmon-pole models), and the algorithm for solving the quasi-particle equation. Primary results include reference values for future benchmarks, best practices for convergence within a particular approach, and average error bars for the most common approximations.


INTRODUCTION
Computational spectroscopy is developing into a complementary approach to experimental spectroscopy.It facilitates the interpretation of experimental spectra and can predict properties of hitherto unexplored materials.In computational spectroscopy, as in any other theoretical discipline, the first step is the definition of the physical model.In theoretical physics and chemistry, this model is governed by a set of equations.Computational sciences solve these equations numerically, and the solutions should of course be independent from the computational settings.However, in reality this is not always the case.In practice, the equations are often complicated, and the numerical techniques introduce many, often interdependent, computational parameters.A thorough validation of these parameters can be very time-consuming.
To validate computational approaches, theoretical benchmarks are essential.In quantum chemistry, benchmark sets are well-established (e.g., G2/97, 1−3 GMTKN30, 4 ISO34, 5,6 S66 7 ).In solid state physics, a validation benchmark set for elementary solids has only recently been published for ground state properties calculated in density-functional theory (DFT). 8ccording to refs 8 and 9, several DFT codes differ by a surprising amount even for the computationally efficient semilocal Perdew−Burke−Ernzerhof (PBE) functional. 10he physical model we address in this article describes charged electronic excitations, and we apply it here to molecules.We focus on Hedin's GW approximation 11 where G is the single particle Green's function and W the screened Coulomb interaction.−53 In its simplest form (G 0 W 0 ), the GW approach is applied as a correction to the electronic spectrum of a noninteracting reference Hamiltonian, such as Kohn−Sham DFT or Hartree− Fock 11−14 (in the following denoted @reference).However, despite G 0 W 0 's more than 50 year history and, starting 30 years ago, its practical implementation within an electronic structure framework, 54−56 results from different codes and approximations have rarely been directly compared.In this article, we provide a thorough assessment of G 0 W 0 , with a particular starting point, for gas-phase molecules using three different codes, validating different computational implementations and elucidating best practices for convergence parameters, such as basis sets, treatment of the unoccupied subspace, and discretization meshes.
In this work, we make the first step and establish a consistent set of benchmarks of ionization energies and electron affinities of 100 molecules−the GW100 set.We present converged G 0 W 0 calculations based on the Perdew−Burke−Ernzerhof (PBE) 10 generalized gradient approximation to DFT.Our G 0 W 0 @PBE results for these molecules can serve as a reference for future G 0 W 0 implementations and calculations.We apply three different G 0 W 0 codes in this work: TURBOMOLE, 49,57 FHIaims, 40,58 and BerkeleyGW. 59The three codes differ in their choice of basis set (atom centered orbitals in TURBOMOLE and FHI-aims and plane waves in BerkeleyGW) and in their implementation.The validation process was crucial to remove conceptual and numerical inconsistencies from our implementations and to test the influence of all computational settings.In the end, all three codes agree on average within 0.2 eV for ionization energies and electron affinities.TURBOMOLE and FHI-aims are all-electron codes and use the same basis sets in this work.(The calculations perfomed in this work are all explicitly nonrelitivistic to exclude the effects of different relativistic approaches.The QZVP basis sets for fifth row elements used in this work are also fully all-electron and do not contain an effective core potential.)They agree to ∼1 meV.BerkeleyGW, employing a real-axis full-frequency method and pseudopotentials, leads to results that differ from TURBO-MOLE and FHI-aims by 200 meV on average.We consider this difference in residual discrepancy as acceptable for the time being because FHI-aims and TURBOMOLE use the exact same (local-orbital) basis set, and BerkeleyGW uses a plane wave basis set with pseudopotentials.The plane wave basis set makes BerkeleyGW also applicable to extended systems.Whether or not there is another dominating source remains to be investigated.
In the GW100 set, we also supply experimental ionization energies and electron affinities, where available.These are intended for future reference.Care has to be taken in the comparison to experimental values, because experimental data tends to carry uncertainties that are intrinsic to the measuring process or reflect external influences (e.g., defects or disorder) and other environmental parameters (e.g., temperature).These effects, and intrinsic effects such as the zero-point motion, are not included in our current theoretical approach.An assessment of the GW method as such is beyond the scope of this work and would require an extensive study of the starting point dependence of G 0 W 0 . 14,39,42,43,47,51,60,61Indeed, the dependence on the starting point is a well-known problem of perturbative GW calculations for molecules as well, which can be solved only through fully self-consistent GW calculations. 29,39,62,63The starting-point dependence for local and semilocal DFT functional starting points are usually small (<0.1 eV 49 ).Considering the whole range of starting-points from semilocal DFT to Hartree−Fock, the starting-point dependence can easily exceed 1.0 eV. 29,49,64In this work, we do not investigate the starting point dependence and only compare results obtained using PBE as a starting point.
Before moving on to discuss the actual results of the GW methods described above, we comment on the differences in computational cost between the G 0 W 0 method and other approaches that may be employed in the evaluation of the quasi-particle excitations of molecules.−68 The high numerical cost of G 0 W 0 calculations mostly stems from the computation of the polarizability and of the screened Coulomb interaction.The precise computational bottlenecks that limit the applicability of the G 0 W 0 method to large systems depends on the type of basis function that is adopted in the calculation.In plane wave-codes, the main bottleneck arises from the necessity to include a large vacuum region to avoid spurious interactions between the fictitious periodic replicas of the system.As a consequence, a large number of reciprocal-lattice vectors are required to numerically converge the calculations, which makes the computation of the polarizability very demanding.In localized-basis code, such as FHI-aims and TURBOMOLE, on the other hand, the most demanding computational operation is the evaluation and storage of the Coulomb integrals, which are typically handled within the framework of the resolution-of-the-identity.
Numerical implementations of the G 0 W 0 approach typically scale as N ( ) 3 or N ( ) 4 with N being the number of basis functions.The scaling with system size is the same.Despite their high numerical cost, G 0 W 0 calculations remain considerably cheaper than wave function-based approaches, such as coupled-cluster singles-doubles (CCSD), CCSD with triples (CCSDT), or configuration interaction (CI).In these methods, the main computational bottlenecks arise from the complex ansatz adopted for the many-body wave function, which requires the explicit inclusion of the excited Slater determinants.Correspondingly, the scaling of computational cost relative to system size is typically N ( ) 6 for CCSD, N ( ) 8 for CCSDT, and exponential for full CI.
The remainder of the article is structured as follows: We start by describing the test set that will be used in this paper in section 2. In section 3, we present the G 0 W 0 approaches used in this work and explain their similarities and differences.In section 4, the ionization energies and electron affinities of the GW100 are presented.In section 5, we discuss the different ways to treat the analytic structure of W in our three G 0 W 0 approaches.The main conclusions are summarized in section 6.

THE GW100 SET
The 100 molecules in the GW100 set include different elements and thereby cover a considerable range of ionization potentials (from ∼4 eV for Rb 2 to ∼25 eV for He).The selected molecules exhibit a spectrum of typical chemical bonding situations.For instance, for carbon, we include a variety of covalent bonds, such as C 2 H 6 , C 2 H 4 , C 2 H 2 , C 6 H 6 , CO, CO 2 , and C 4 .Special interest is devoted to bonds of metal atoms by including Cu 2 and Ag 2 , Li 2 , K 2 , Na 2 , and Rb 2 as well as small metallic clusters Na 4 and Na 6 .In contrast, alkaline metal halides are prototypes for ionic bonds with LiF being an extreme and KBr a more moderate case.The alkaline earth metal compounds MgF 2 and MgO are also ionic, the former with a vanishing and the latter with a large dipole moment.Including series of homologous like N 2 -P 2 -As 2 , F 2 -Cl 2 -Br 2 -I 2 , or CF 4 , CCl 4 , CBr 4 , or CI 4 facilitates the identification of trends within a group of elements.These trends can then be correlated with certain physical or chemical properties, such as the decreasing ionic bond character in the last example.Furthermore, we have also included several simple organic molecules, like alcohols, aldehydes, and nitrogenous bases, as well as the most typical test cases often appearing in benchmark sets such as water and carbon mono-and dioxide.
We have considered using a standard test set, for example, the G2 set.This was not done for the following reasons: First, some of the GW codes are restricted to closed-shell systems, and the G2 set contains many open-shell systems, for example, the alkaline atoms.Second, it is restricted to compounds of the first two periods, excluding the elements Li, Be, Na, Mg, and Al and thus covers neither bonds between metal atoms nor pronounced ionic bonds.Third, it contains only small molecules in which the degree of delocalization is rather small.In our work, we observe that none of the nine compounds that show a discrepancy between GW resulting in larger than 1 eV in the IP at GW-level is contained in the G2 set.The respective systems turn out to be strongly ionic alkaline (earth) hydrides, fluorides, chlorides, bromides, and oxides as well as C4 and O3.The worst case that is also present in the G2 set is FH with a maximal discrepancy of 0.8 eV.We would thus not expect any other G2 molecule to show a larger deviation.
The molecular geometries used in this work are mainly taken from experimental data.For some molecules, the final structure was obtained by optimizing a known morphology using DFT in the PBE approximation for the exchange-correlation functional using the def2-QZVP basis set.All molecular geometries are included in the Supporting Information.

COMPUTATIONAL METHODOLOGY
The objective of this work is to establish a set of converged and validated G 0 W 0 results for molecules that will serve as a benchmark for future work.Specifically, we calculate the G 0 W 0 self-energy 11 where G 0 σ denotes the one-particle causal Green's function for spin channel σ = ↑,↓ and W 0 is the screened Coulomb interaction in the random-phase approximation (RPA).Traditionally, one splits the self-energy into energy dependent correlation and energy independent exchange terms as G 0 σ is given in terms of the single-particle wave functions ψ nσ KS (r) and eigenvalues ϵ nσ KS of a reference Kohn−Sham DFT calculation. 69 where ϵ F is the Fermi level (chemical potential) and η a positive infinitesimal.W 0 follows from the noninteracting response function χ 0 that in a real-space representation assumes the following form where the second line is the sum-over states representation of Adler and Wiser, 70,71 and f m and f n are the occupation factors of states m and n, respectively.Given the noninteracting response function χ 0 , W 0 is expanded in powers of χ 0 where we have omitted the space and frequency variables for simplicity.The last equal sign in eq 6 introduces the reducible response function Given the G 0 W 0 self-energy, the quasi-particle (QP) energies ε n QP are computed by solving the diagonal QP equation in the basis of the single-particle states |nσ⟩, i.e.
In the above equation, v xc denotes the exchange-correlation potential of the preceding DFT calculation.Since we will restrict ourselves to closed shell molecules for the remainder of the paper, we drop the spin index σ.
In our comparative study, we use a DFT-PBE starting point.Once we have decided on this reference Hamiltonian, G 0 is uniquely defined.As a result, Σ in the G 0 W 0 approximation is also uniquely defined and all G 0 W 0 implementations should, in principle, produce the same G 0 W 0 @PBE results.However, G 0 W 0 implementations can differ in several aspects, that in practice can lead to deviations.The most critical aspects are listed as follows: (a).The Choice of Basis Set.In this work, we will compare GW results obtained with both local orbital (LO) and plane wave (PW) basis sets.In a local basis set, the size of the Hamiltonian is usually significantly smaller than in plane waves.For example, for ethene with the LO basis used here, approximately 350 functions were used; the analogous plane wave calculation used a factor of 225 more functions.However, for an LO basis, there is no unique recipe to systematically increase the basis set to approach the complete basis set limit.Conversely, a plane wave basis set is conceptually easier to converge: one main parameter, the kinetic energy cutoff, needs to be increased until convergence is achieved.However, plane waves require periodic boundary conditions, and molecules therefore have to be placed in a supercell, which is a large, periodically repeated box that is filled with vacuum.Although the size of the supercell usually convergences quickly for local and semi-local DFT functionals, the screened Coulomb interaction W 0 in G 0 W 0 can be long-ranged.A brute force supercell convergence is therefore computationally impractical. 72Instead, the Coulomb interaction is truncated, 73,74 which introduces, in principle, an additional convergence parameter.In BerkeleyGW, this truncation is automatically based on supercell size.To compute ionization energies and electron affinities, we correct for the shift in vacuum level that is present due to periodic boundary conditions.Taking the molecule at the center of the supercell, the vacuum correction is calculated as the electrostatic potential at the supercell edges.
−77 These energy cutoffs can reduce the computational cost, albeit at the expense of new convergence parameters. 78,79onversely, the use of a small number of basis functions in local orbital calculations can limit the number and the accuracy of empty states above the vacuum level.If in this situation a systematic extrapolation to the basis set limit is not feasible, this can lead to limitations in the accuracy of the response function given by eq 5 and G 0 and W 0 and, therefore, the self-energy itself.
G 0 W 0 has also been implemented in other basis sets, such as projector augmented waves (PAW), see refs  (b).The Treatment of Core Electrons and Valence Electrons in the Core Regions.In the FHI-aims and TURBOMOLE calculations in this work, we included all electrons explicitly at each step of the calculation.However, at the moment, a pure plane wave basis set requires pseudopotentials that remove the rapid oscillations of the electronic wave functions near the nuclei and hence drastically reduce the required energy cutoff.In a pseudopotential approach, the electrons are divided into core and valence electrons.−92 On the DFT level, pseudopotentials are typically derived from single-atom DFT calculations performed with the same functional as used later.
However, interactions between core and valence electrons that depend on the environment in a polyatomic system undermine transferability and lead to deviations between pseudopotentials and all-electron calculations.The use of pseudopotentials can lead to errors in GW calculations because of the neglect of core polarization effects, 14,93 deviations in the pseudo and all-electron wave functions near the nucleus, 94,95 and core−valence interactions. 94,95To shed more light on the quantitative role of core electrons and pseudopotentials, we directly compare all-electron, frozen core, and pseudopotential G 0 W 0 calculations in this work.
DFT-PBE pseudopotentials are used in our planewave G 0 W 0 calculations.All valence electrons of the same principal quantum number are treated on the G 0 W 0 level, whereas the description of the core electrons and the core−valence interaction remains on the DFT level.
(c).Treatment of the Frequency Dependence.All quantities in G 0 W 0 depend explicitly on a frequency (or time) argument (see section 3).Different G 0 W 0 implementations differ in their treatment of this frequency dependence.Because the poles in G 0 and all subsequent quantities lie close to the real frequency axis, all quantities in G 0 W 0 exhibit a pronounced fine structure on the real frequency axis, whose resolution requires fine frequency grids and a large number of frequency points.Different strategies are employed to avoid the computational bottleneck of dense frequency grids.In this work, we will compare several different ways: the fully analytic treatment in TURBOMOLE (TM-RI) and (TM-noRI), 49 an integration on the imaginary axis with subsequent analytic continuation to the real axis as implemented in FHI-aims (AIMS-2P) and (AIMS-P16), 40,58 and the BerkeleyGW implementation of the full frequency treatment on the real axis (BGW-FF) 96 and of a generalized plasmon pole model (BGW-GPP). 55Techniques employing contour deformation and approaches to circumvent the sum over empty states have also been reported in the literature, 97−101 but they are not considered in this work.
(d).Solution of the QP Equation.−104 Because not all GW codes search for all solutions, different G 0 W 0 implementations may give different answers for multiplesolution cases even though the underlying self-energies may be very similar.
In summary, the impact of these differing approaches (a−d) can affect the QP energies in a significant way, which has only in part been quantified by previous studies. 40,45,49,94,95Apart from code validation, a second main goal of our work is a quantitative comparison of different methods.This goal will be achieved by comparing the G 0 W 0 results (i.e., the QP energies) obtained from three different G 0 W 0 implementations: TUR-BOMOLE, FHI-aims, and BerkeleyGW.
In the following sections, we describe the conceptual and technical differences that distinguish the TURBOMOLE, FHIaims, and BerkeleyGW G 0 W 0 implementations.A more general description of the GW method and its application to molecules in particular can be found in various reviews [12][13][14]49,105,106 and is not the topic of this paper.We will also provide detailed convergence studies of the relevant computational parameters.
3.1.The Frequency Dependence of the G 0 W 0 Self-Energy.A pronounced difference between different G 0 W 0 implementations is the treatment of the frequency dependence of the self-energy in eq 1 and of intermediate quantities.In this work, we will compare an analytic treatment of the pole structure facilitated by the spectral representation of the response function with a plasmon-pole model and a numerical real as well as imaginary frequency treatment.In the next four sections, we describe the technical aspects of these approaches.The implications on the results of these different approaches will be discussed in the Results and Discussion sections.
3.1.1.Implementation of the Fully Analytic (FA) Spectral Representation in TURBOMOLE.The RPA response function (eqs 5 and 6) is calculated explicitly in its spectral representation.Because the Green's function G 0 has N occ + N unocc poles, the screened interaction W 0 exhibits 2 × N occ N unocc poles.As the exact pole positions of W 0 are inherited from G 0 and therefore known, we can evaluate the energy integral for Σ (eq 1) analytically.This gives for Σ 2 × (N occ + N unocc ) N occ N unocc poles.In the rest of this section, we summarize the most important technical details following ref 49, where we here focus on the nonmagnetic case ψ = ψ ↑ = ψ ↓ .
The implementation of G 0 W 0 in TURBOMOLE is based on the spectral representation of the reducible response function.
The pole positions, Ω m , are the (charge neutral) excitation energies, and ρ m (r) denotes transition densities.The ρ m are expanded in a basis of orbital products, where i,j,.. label occupied states, and a,b,.. label empty states.The vectors |X m ,Y m ⟩ are solutions of the eigenvalue problem ) , 0 m m m (10)   under the orthonormality constraint The operators contain the orbital rotation Hessians (14)   with . From the reducible response function, the screened Coulomb interaction W can be directly constructed by contracting with the Coulomb interaction v The self-energy can finally be obtained directly by performing the energy integral analytically because in this formalism the energy structure of both G and W is known.Performing the integral leads to a closed expression for the matrix elements of the self-energy.The real part of the diagonal matrix elements of Σ includes the exchange contribution while for the correlation contribution, we have where, as before, m runs over all density excitations, and η̃is an infinitesimal.
3.1.2.BerkeleyGW Implementation of the Full Frequency (FF) Integration along the Real Frequency Axis.The FF-BerkeleyGW approach evaluates the frequency-dependent selfenergy numerically along the real frequency axis.Similar to the TURBOMOLE implementation, G 0 , W 0 , and Σ retain the full pole structure.However, unlike in TURBOMOLE, the pole structure is represented on the frequency grid.
In FF-BerkeleyGW the frequency dependent self-energy is where W r/a are the retarded (r) and advanced (a) screened Coulomb matrix, which can be expressed as where the retarded/advanced dielectric function ϵ r/a has the form and the retarded/advanced reducible polarizability is defined as ∑ ∑ where Δϵ i,j = ϵ j − ϵ i ≥ 0. The integral in eq 18 is evaluated numerically on the interval [0,ω high cut ], beyond which the integrand is negligibly small.We further reduce computational cost by dividing the interval [0,ω high cut ] into two intervals: [0,ω low cut ] and [ω low cut ,ω high cut ] treated with two different integration schemes.The low frequency cutoff ω low cut is chosen to ensure that all poles of the numerator in the integrand lie below ω low cut ; therefore, the integrand decays smoothly beyond this point.A uniform fine frequency step Δω is used on [0,ω low cut ], whereas a smaller number of quadrature points can be used to perform a standard numerical integration of eq 18 within the interval [ω low cut ,ω high cut ].Because the integrand contains a number of singularities in the interval [0,ω low cut ], we use the special integration scheme of ref 107 to perform the numerical integration accurately on this interval.
For all 18 of the molecules studied with FF-BGW, we use η = 0.2 eV and Δω = 0.2 eV.The left panel of Figure 1 shows a typical example (ethylene molecule) of the convergence with the dependence of the quasi-particle energy on the parameters η and Δω for predicted QP-HOMO energy.Because our aim here is to test the convergence behavior for η and Δω, we used a low E ϵ cut of 5.0 Ry (68 eV).(In principle, we use eV as the energy unit; however, in cases where the input parameters of a specific code are supplied in another unit, we also provide the quantity in this unit.)Figure 1 shows that for a given η, Δω = η is sufficient for an accuracy of 0.01 eV, and that η = 0.2 eV is sufficient.The right panel of Figure 1 shows an example for a molecule with a large slope of the self-energy.Here, a convergence within 0.2 eV is reached.
In a plane wave basis, the number of basis functions is much larger than in the localized basis set used in TURBOMOLE and FHI-aims.Therefore, intermediate quantities, such as the dielectric function and the self-energy, require larger matrices.To keep the calculation tractable, we only computed the dielectric function for a reduced set of plane waves with plane wave cutoff energy E ϵ cut .In addition, the number of unoccupied states that enter the Green's function and the self-energy is controlled by the energy of the highest unoccupied state, E max .Because E ϵ cut and E max are interdependent, 37 careful convergence studies are required.Here, we converged these parameters within GPP-GW as will be described below.We converged the dielectric function by increasing both E max and E ϵ cut until the screened exchange component of the self-energy changed by less than 0.2 eV and extrapolated the Coulomb-hole term of the self-energy by means of the static completion method. 108ithin this approach, a term is added to the Coulomb-hole component of Σ, which corrects for the truncation of the number of unoccupied states. 108n Table 1, the values of the parameters introduced in this section are listed for 18 molecules studied with FF-BGW.For all FF-BGW calculations, values of η = 0.2 eV, Δω = 0.2 eV, and E ϵ cut = 24 Ry (326 eV) have been used.

FHI-aims Implementation of the Imaginary Frequency Treatment Including Analytic Continuation (AC).
In imaginary frequency implementations, G 0 , W 0 , and Σ retain their full pole structure but on the imaginary frequency axis.This representation requires significantly fewer frequency points due to the smooth behavior of all quantities on the imaginary frequency axis.The self-energy is then analytically continued to the real frequency axis.The "quality" of the selfenergy on the real-frequency axis depends on the type of analytic continuation.
FHI-aims provides two different analytic continuation models: (i) a "two-pole fit" 109 and (ii) a Padéapproximation (see ref 110).In the "two-pole fit", each matrix element of the self-energy is fitted to the following expression on the imaginary axis where the dependence of a and b on i and j has been omitted for simplicity.In the Padéapproximation, Σ ij is given by where N is the total number of parameters employed in the Padéexpansion.The ionization energies and electron affinities in Tables 2 and 5 have been calculated using N = 16 (which is equivalent to a sum of 8 poles).Σ on the real-frequency axis is then obtained by replacing iω with ω in eqs 21 or 22, and the G 0 W 0 quasi-particle energies are obtained by solving eq 7.
3.1.4.Generalized Plasmon Pole (GPP) in the BerkeleyGW Implementation.Within our GPP implementation, 59 the expression for the self-energy, Σ, is written as the sum of two terms, termed the screened-exchange (SX) and Coulomb-hole (CH).Here, and where j runs over occupied states, n″ runs over both occupied and unoccupied states, and Ω, ω̃, λ, and Φ are the effective bare plasma frequency, GPP mode frequency, amplitude, and phase of the renormalized Ω , respectively, defined in reciprocal space as Here, ρ is the electron charge density in reciprocal space, and ω p 2 = 4πρ(0)e 2 /m is the classical plasma frequency.The convolution of G 0 and W 0 is performed analytically.Finally, Σ contains the same number of poles as G 0 in this setting.
As mentioned above, GW calculations within any basis require thorough convergence studies because the dielectric function and self-energy converge slowly with respect to the energy of the highest unoccupied state (E max ) included in the calculation as well as the G-vector cutoff for the dielectric matrix (E ϵ cut ).Moreover, these parameters are interdependent. 37e tested the convergence of the quasi-particle energies by varying E max such that the unoccupied states reached 30, 45, 90, 120, and 150 eV above the vacuum level while simultaneously varying E ϵ cut from 82, 163, 327, and 408 to 544 eV.We note that E max defines the number of unoccupied states used to construct the dielectric function and the self-energy, whereas E ϵ cut affects the self-energy only through its effect on the dielectric function.
Because of the computational cost of these convergence studies, we systematically checked the convergence of the calculated ionization energies and electron affinities of only 35 randomly chosen molecules and applied consistent convergence parameters to all 100 molecules.An example of the convergence of the ionization energy is shown for magnesium oxide and ozone in Figure 2. We determined that, at E max = 90 eV and E ϵ cut = 408 eV, the error was less than 0.2 eV for all but the ozone molecule (see Figure 3) when compared to the highest convergence criteria tested (E max = 150 eV and E ϵ cut = 544 eV).Therefore, this was set as the convergence criteria for all 100 molecules.When computationally feasible (for 10 additional molecules), we checked that going beyond this convergence criteria did not alter predicted energies for the remaining molecules by more than 0.2 eV.The ozone molecule was exceptionally difficult to converge and required a higher number of unoccupied states.The reported value of ionization energy and electron affinity of ozone was E max = 120 eV and E ϵ cut = 544 eV.
In our pseudopotentials for transition metal atoms, semicore d-states are explicitly treated as valence.Here, the density that is used to construct the plasma frequency does not include the semicore, but inclusion of semicore electrons is necessary for describing the nodal structure of the valence wave function. 111ecause the planewave cutoff is very high with the use of semicore states, we performed a limited convergence study for Cu 2 of the influence of dielectric function cutoffs and number of bands.We estimate an error of ∼0.3 eV for the predicted ionization energies and electron affinities of these transition metal-containing compounds.
3.2.Basis Sets.Types of Local Basis Sets.The G 0 W 0 QP energies in TURBOMOLE and FHI-aims are calculated using the TURBOMOLE def2 basis sets of contracted Gaussian orbitals. 112The QZVP basis sets for fifth row elements used in this work are also fully all-electron and do not contain an effective core potential.In the TURBOMOLE calculations, the contracted Gaussians are treated analytically, exploiting the standard properties of Gaussian functions.With FHI-aims, in contrast, the Gaussian orbitals are treated numerically to be compliant with the numeric atom-centered orbital (NAO) technology of FHI-aims. 40,58NAOs are of the form where u i (r) are the numerically tabulated radial functions, and Y lm (Ω) is the spherical harmonics.
In the Supporting Information, we provide a convergence study for the most critical numerical parameters in TURBOMOLE and FHI-aims for four representative molecules.The quantities that only depend on the occupied KS-  reference states, the KS-energies, the (matrix elements of the) exchange-correlation potentials, and the exchange part of the self-energy Σ x are converged at the meV level and agree between FHI-aims and TURBOMOLE at this level.The correlation part of the self-energy Σ c also depends on the unoccupied states and introduces additional convergence parameters that depend on the method that is used to calculate it.Their convergence will therefore be discussed separately.Because the ground state quantities obtained with TURBO-MOLE and FHI-aims are practically identical, we will conduct a basis-set convergence study for the ground state only once using results obtained with TURBOMOLE.
Convergence of the KS States in the Local Basis Sets (DFT).The convergence of the KS-HOMO with respect to the size of the basis set is presented in Figure 4.The KS eigenvalues calculated with def2-SVP, def2-TZVP, and def2-QZVP basis sets are compared to the "complete" basis set extrapolation.The extrapolation is obtained from the def2-TZVP and def2-QZVP results by a linear regression against the inverse of the total number of basis functions.This extrapolation technique was applied previously in ref 49 and described and validated there in more detail.The KS energies converge systematically with increasing basis set size (details are provided in the Supporting Information).The molecules containing fifth row elements are excluded from the studies because there are no allelectron SVP and TZVP basis sets available for these elements.
Most molecules are seen to converge from below.At the QZVP level, the KS energies of 66 (88) out of the 100 molecules are converged to within 50 meV (100 meV).We note that the def2 basis sets are optimized for Hartree−Fock total energies.There, convergence for KS eigenvalues might thus be slower.
Convergence of the QP States in the Local Basis Sets (G 0 W 0 ).The convergence of the QP energies with respect to the size of the basis sets is illustrated in Figures 5 and 6.In the upper part of Figure 5 we report the difference of the QP-HOMO energies to the basis set extrapolated results for the def2-SVP, def2-TZVP, and def2-QZVP basis sets.The QP-LUMOs are shown in the lower part of Figure 5.For the QP-HOMO levels, the same data is replotted in Figure 6 as a scatter plot.All results have been calculated with FHI-aims.TURBOMOLE gives identical results (see section 4).As expected, the QP energies converge slower with basis size than the Kohn−Sham energies.For the QP-LUMO levels, we observe a similar pattern (see Figure 5).However, the  deviations are approximately a factor of 2 larger for each basis set.
When comparing the QP results between the local orbital codes FHI-aims and TURBOMOLE and the plane wave code BerkeleyGW in the Results section, we will always use the extrapolated values.These will be referred to as "EXTRA".For the QP energies, we estimate the error in the extrapolation by comparing to a second extrapolation scheme, i.e., extrapolating against 1/C n 3 where C n is the cardinal number of the basis set (2 for SVP, 3 for TZVP, etc.).The mean absolute error is 0.04 for QP-HOMO and 0.1 for QP-LUMO.In the Supporting Information, a systematic comparison against other basis sets (including Dunning's correlation consistent basis sets 113 ) is given for four typical molecules.In the Results section, the estimated errors are provided together with the extrapolated results.
Auxiliary Basis Sets.The G 0 W 0 implementations in TURBOMOLE and FHI-aims make use of the resolution-ofidentity (RI) technique to compute four-center Coulomb integrals of the type To avoid the numerical cost associated with the computation and storage of the (ij|kl) matrix, an auxiliary basis set {P μ } is introduced to expand the product of basis function pairs as where C ij μ are the expansion coefficients, and N aux is the number of auxiliary basis functions P μ .The G 0 W 0 equations can be conveniently rewritten employing eq 30 as described in detail in refs 40, 49, and 114.In our experience, RI speeds up DFT calculations by an order of magnitude in TURBOMOLE.For G 0 W 0 , the computational effort is reduced by an order of magnitude for calculations with ∼100 basis functions.In addition, the scaling reduces from N 5 to N 4 or better (where N is the number of atoms) up to 1000 basis functions. 49For G 0 W 0 in FHI-aims, RI is essential. 40The scaling of G 0 W 0 in FHI-aims is always N 4 or better.
The auxiliary basis functions in TURBOMOLE are supplied in a database. 115,116They are designed such that the RI-induced error lies in the meV regime for DFT total energy calculations. 116An important difference for FHI-aims is that the auxiliary basis functions in FHI-aims are constructed on the fly 40,114 and are not predefined as in TURBOMOLE.
We define the G 0 W 0 RI error as the difference between the QP energies calculated with and without applying RI in all steps of the calculation.This comparison is shown in Figure 7 for TURBOMOLE using a subset of GW100.We observe that G 0 W 0 is more sensitive to the RI approximation than DFT with local or semi-local functionals.In an earlier study, we found for a smaller set of molecules that the G 0 W 0 RI error for the QP-HOMO energies is below 0.1 eV. 49The same trend is observed here (see Figure 7).Only very small systems, such as helium and hydrogen, tend to have a larger G 0 W 0 RI error of 0.24 and 0.13 eV, respectively.In the FHI-aims calculations, the parameters for constructing the auxiliary basis sets are chosen such that the quasi-particle energies agree with the RI-free TURBOMOLE values to better than 1 meV compared for all systems (see also Table 2).[All input and output files of the FHI-aims calculations are available at https://NOMAD-Repository.eu DOI: 10.17172/NOMAD/2015.11.03-1.] Convergence in the Plane-Wave Basis.The BerkeleyGW 59 package computes the dielectric function and self-energy within a plane-wave basis set.The input DFT-PBE eigenvalues are computed with the Quantum Espresso package 117 with a planewave wave function cutoff defined such that the total DFT-PBE energy is converged to <1 meV/atom.The wave function cutoffs for all 100 molecules are provided in the Supporting Information.Typical values between 50 (680) and 120 Ry (1633 eV) are sufficient, but in some cases, even 300 Ry (4082 eV) are necessary.The molecules are placed in a large supercell that is twice the size necessary to contain 99.9% of the charge density.To avoid spurious interactions between periodic images at the G 0 W 0 step, the Coulomb potential is truncated at half of the unit cell length. 59he KS eigenvalues computed with Quantum Espresso agree well with TURBOMOLE and FHI-aims (see Table S6 in the  The RI error for the G 0 W 0 calculations is defined as the difference between the QP-energies calculated with and without applying the RI approach at all steps of the calculation.Supporting Information).In most cases, the deviation is well below 0.1 eV, the mean absolute deviation is 0.048 eV with the extrapolated and 0.045 eV with the QZVP results.The largest discrepancies, on the order of 0.1 eV, occur almost exclusively in Fluor-containing molecules.
As alluded to in section 3, we obtain the QP energies by solving the QP equation (eq 7, repeated here for convenience) We observe that Σ x − v xc also differs by less than 50 meV for the GW100 molecules between the three codes (see Table S6 in the Supporting Information).
The Response Function in the Local Basis.Because of the compactness of the basis in TURBOMOLE and FHI-aims, the response function can be treated in the full Hilbert space of the def2-QZVP basis for all molecules in GW100.Even for the smallest molecules, the def2-QZVP basis describes a homogeneous excitation spectrum up to 350 eV.For all molecules with more than 4 electrons, the excitation spectrum ranges well beyond 800 eV.We tested that at least states up to 400 eV are required in the construction of Σ (i.e., in the sum over m in eq 17) for all systems.The QP-HOMO energies are then converged to within 0.05 eV.For states with energies higher than 400 eV, the denominator in eq 17 becomes large enough to suppress further contributions.
We mention this assessment here to estimate the contribution of higher lying states in G 0 W 0 .In practice, we always include all states of the Hilbert space spanned by the chosen basis set in our TURBOMOLE and FHI-aims calculations.
3.3.Treatment of Core Electrons.In the FHI-aims and TURBOMOLE calculations, all electrons have been taken into account explicitly at each step of the calculation.They are included fully in the KS and the G 0 W 0 calculations and thus also take part in the screening.
3.3.1.Pseudopotentials.In the BerkeleyGW calculations, we employed Troullier−Martins norm-conserving pseudopotentials. 118For all atoms except the transition metals, the pseudopotentials were taken from version 0.2.5 of the Quantum Espresso pseudopotential library. 117For Ti, Cu, and Ag, we found it necessary to include semicore states in the pseudopotential to properly describe the valence orbitals near the core. 119For Cu and Ag, we use the pseudopotentials of Hutter and co-workers 120,121 in which 19 valence electrons are treated explicitly.For Ti, we generated the pseudopotential using the FHI98pp package, 87 treating 12 electrons as valence.The plane-wave cutoff is set such that total DFT energies are converged to 10 meV/atom and the DFT HOMO energies are converged to 50 meV for all molecules.As shown in Figure 8, the HOMO energies are converged to <5 meV for the majority of the molecules (79 molecule).The error is computed as the difference between the converged cutoff energy and a cutoff of 120 Ry (1633 eV).For cases where the converged cutoff energy is determined to be greater than 120 Ry (1633 eV), the error is estimated as the difference between this value and a planewave cutoff, which is increased by 10 Ry (136 eV).The core radii cutoff for the pseudopotentials and the wave function cutoff for all 100 molecules are listed in the Supporting Information.This same planewave cutoff is used at the GW step.
3.4.The Quasi-Particle Equation.The final technical step of a G 0 W 0 calculation is the solution of the quasi-particle equation.Once the G 0 W 0 self-energy is obtained, the quasi-particle energies ε n QP are calculated by solving the diagonal quasi-particle equation (eq 7, repeated here for convenience) Here, ϵ n KS are the Kohn−Sham eigenvalues (computed in PBE in this work).Eq 32 is generally solved iteratively due to the interdependence of the self-energy Σ(ε n QP ) and the quasiparticle energy ε n QP .In most cases, ε n QP falls in a region in which the self-energy has no poles, making it featureless and almost constant.Then, the solution is unique in the region of interest, and eq 32 may even be linearized such that a single evaluation of Σ at the KS energy is sufficient and there is no iteration process.In BerkeleyGW GPP calculations, the quasi-particle equation is always linearized.This is justified because the GPP self-energy is smooth near the quasi-particle energy.
In some cases, however, the initial KS energy is close to a pole of Σ, and as a result, eq 32 has more than one solution.These solutions are relatively close in energy, which implies that the correction to the KS eigenvalue is not unique.A schematic example is shown in Figure 9.In practice, most available G 0 W 0 codes search for only a single solution of the quasi-particle equation.Which solution is found depends on the initialization and the type of the iterative procedure.In general, different codes might find different solutions to the nonlinear quasi-particle equation even though the self-energy is similar.We expect all solutions to be physically relevant in principle.However, which one is actually physically relevant depends, e.g., on the quasi-particle weight (i.e., the pole residue), and furthermore may also depend on the physical observable one would like to study.
Generically, one expects those solutions with the largest quasi-particle peak Z/(ω − ε n QP ), the so-called Z-factor Figure 8. Error in the computed DFT HOMO eigenvalues due to the chosen planewave cutoff.As noted in the text, the error is calculated as the orbital energy difference for a calculation that employs the chosen cutoff and 120 Ry (1633 eV) cutoff.If the converged cutoff is greater than or equal to 120 Ry (1633 eV), the error is determined as the difference between the converged value and 10 Ry (136 eV) higher than this value.
to be most important.Especially very close to the poles, the slope of Σ at the quasi-particle energy is large, Z is small, and the weight of the quasi-particle peak in the spectral function is reduced.The remaining weight is shifted to the other solutions, in the case of multiple solutions, or the incoherent background.Solutions that result from intersections with the almost vertical lines when a pole in Σ changes sign are spurious in a sense.For the molecules considered here, these poles are very sharp, as we will demonstrate in the Results section.The slope is thus almost infinite and the corresponding weight will go to zero.A lesser Padéapproximant or large damping η, as often applied in a full frequency treatment, will broaden the pole and can give rise to a nonvanishing spurious Z-factor.At this point, we will refrain from further analysis of the role of the different solutions for physical observables.In this work, we ascribe the solution with the highest energy to the QP-HOMO.
One may conclude now that states whose derivative of the self-energy is large at the quasi-particle energy, i.e., that have a small Z-factor, might lend themselves to multisolution behavior.If a small Z-factor is detected, it might indeed be necessary to test whether the solution of the quasi-particle equation actually corresponds to the highest in energy intersection between the Δ(ω) and Σ(ω).A small Z-factor, however, is not always an indicator of multisolution behavior.We will see examples of this below.

RESULTS
In the following, we present our numerical assessment of the G 0 W 0 implementations in TURBOMOLE, FHI-aims, and BerkeleyGW.We have chosen the energies of the highest occupied molecular orbital (QP-HOMO, ε H QP ) and lowest unoccupied molecular orbital (LUMO, ε L QP ) as our main observables in this comparison.Both have a well-defined physical meaning: ε H QP corresponds to the first vertical ionization energy, and ε L QP is the electron affinity.For most molecules in the GW100 set, the quasi-particle equation yields a unique solution for QP-HOMO and QP-LUMO.However, four molecules exhibit the aforementioned multisolution behavior described in section 3.4.Section 4.1.1reports all solutions for these four molecules in detail.In the final subsection of the results section, we will then make a comparison to other G 0 W 0 results that are available in the literature. 2 presents the G 0 W 0 @PBE QP-HOMO energies for the molecules of the GW100 set.(For a subset of the molecules considered in this work, the G 0 W 0 @ PBE ionization energies−evaluated with the three codes−have been reported in previous publications. 37,40,49The results presented in Table 2 show small numerical differences for some of these molecules.The deviations of these values to previously published data are generally smaller than 0.1 eV and can be mainly attributed to the different (larger) basis sets employed in the present study.Furthermore, in the previous TURBO-MOLE calculations, we added an exchange-correlation kernel to the RPA response function and the quasi-particle equation was linearized.The only molecule previously calculated with BerkeleyGW is benzene, for which there is no difference between our current and the previous calculation.)For brevity, we have introduced the following abbreviations in this section: AIMS-2P for a two pole fit in FHI-aims and AIMS-P16 and AIMS-P128 for a 16 and 128 parameter Padéfit in FHI-aims, respectively, BGW-GPP and BGW-FF for the generalized plasmon model and the full-frequency treatment in Berke-leyGW, respectively, and TM-RI and TM-noRI for the RI and the RI-free treatment in TURBOMOLE.respectively.EXTRA denotes extrapolated local orbital results obtained by extrapolating the def2-TZVP and QZVP values calculated using FHI-aims (see section 3.2).The TM-noRI and the BGW-FF calculations are computationally very demanding.They have therefore only been performed for subsets of the GW100 set (see Tables 2 and 5).

Ionization Energies. Table
The absolute values of the differences of the approaches used in this work are reported in Figure 10 for all molecules considered in this work.The FHI-aims and TURBOMOLE results are compared at the QZVP level, and the BerkeleyGW results are compared to the extrapolated results.The AIMS-P16 and TM-noRI QP energies (green shading) generally differ by less than 1−2 meV.There are, however, some molecules for which we observe a larger discrepancy.In these cases, we observe a QP weight Z in the range between 0.6 and 0.8.AIMS-P16 is slightly less accurate in these cases.For the systems for which we observe the multisolution behavior (ozone (index 81), boron nitride (index 65), beryllium oxide (index 84), and magnesium oxide (index 85)), the quasi-particle equation for the QP-HOMO is solved close to a pole of the self-energy.As alluded to in section 3.4, this leads to multiple solutions of the quasi-particle equation even though the underlying self-energies show only minor numerical differences.These cases will be discussed separately in section 4.1.1.
As is also seen in Table 2 and Figure 10, slightly larger discrepancies are observed for the other implementations (see red and yellow shading in Figure 10).AIMS-2P and TM-RI yield values that are generally within 0.1 eV of the AIMS-P16 and TM-noRI results.Similarly, the BGW-FF QP energies agree with AIMS-P16 values within 0.2 eV and BGW-GPP within 0.5 eV.
Table 3 condenses the information given in Figure 10 by reporting the mean deviation (MD) and mean absolute deviation (MAD) of the QP-HOMO energies obtained from the different G 0 W 0 calculations.The mean deviations reported in Table 3 show that TM-RI tends to underestimate QP-HOMO, whereas the generalized plasmon pole model leads to larger QP-HOMO values, overestimating the full-frequencydetermined QP-HOMO for all systems.
The most important conclusion that can be drawn from Tables 2 and 3 is that the AIMS-P16 and TM-noRI results agree to within 3 meV; the numerical aspects in these considerably different implementations are under control.It further shows that the analytic continuation for quasi-particle energies is very accurate, provided that a Padéfit is used and the number of parameters is converged explicitly.We will therefore use AIMS-P16 and AIMS-P128 to extrapolate the QP-HOMO energies to the complete basis set limit, which is also shown in Table 2 in the column EXTRA.This column presents our main result, converged benchmark numbers for the GW100 set.
There is a difference of ∼0.5 eV between the ionization energies from BGW-GPP and TM or AIMS as shown in Table 3.The comparison with BGW-FF shows that the majority of these differences arise from the GPP approximation.For the 18 molecules we tested explicitly (discounting possible deviations arising from the multiple solution behavior), BGW-FF and EXTRA agree to better than 0.2 eV.
4.1.1.Multiple Solutions of the Quasi-Particle Equation.The quasi-particle energies listed in Table 2 have been obtained from the solution of the nonlinear quasi-particle (eq 32).As alluded to before, the solution for the QP-HOMO (and QP-LUMO) is unique for the majority of the systems considered here.However, for the QP-HOMO of ozone, boron nitride, magnesium oxide, and beryllium oxide, we find multiple solutions.Figure 11 illustrates the behavior for ozone for the different calculations.An analogous plot for magnesium oxide is shown in section 5.1, and plots for boron nitride and beryllium oxide are reported in the Supporting Information.The solutions of the QP equation are the intersections of the red line with the self-energy curves.Table 4 reports all solutions we find for the four multisolution molecules.
It should be noted that, technically speaking, each solution of the quasi-particle equation is potentially physically relevant.However, at least two mechanisms can be identified, whereas in practice usually not more than one solution needs to be considered.First, in cases with strong broadening, secondary intersections are suppressed, see, for example, the blue trace in Figure 11.Second, even if the broadening is weak, the slope of the self-energy at the intersection point (i.e., the quasi-particle weight) will be different, favoring in general one point against all others.Only in intermediate situations, as appears to be the case with ozone, for example, could two or more solutions survive.
The absolute value of the difference between the various calculated electron affinities are reported in Figure 12.The analysis of the deviations in Table 6 and Figure 12 shows trends similar to those observed for ε H QP in section 4.1.Again, the agreement between the TM-noRI and AIMS-P16 electron affinities is on the order of a few meV.The TM-RI G 0 W 0 implementation yields electron affinities that are generally within 14 meV from the TM-noRI and AIMS-P16 values.One may conclude from this that QP energies of occupied states are more sensitive to the RI approximation than those of unoccupied states.
In the comparison between the extrapolated and BGW data, we observe a significant difference between systems with a bound and unbound LUMO state.In the case of a bound LUMO state, the agreement between the EXTRA values and BGW results is similar to that observed for the occupied states.For the unbound LUMOs, the agreement can be far off.Here, the local orbital description of the unbound states clearly is not well-converged.It is noted, however, that the convergence is the worst for small systems.Most molecules with more than four none-hydrogen atoms show a small deviation also for unbound LUMOs.
4.3.Literature Comparison.4.3.1.Comparison with Plane-Wave Results.A subset of the GW100 molecules (C 2 H 4 O is present in both sets; in GW100, this is acetaldehyde, whereas in Pham et al.'s set, this is ethylene oxide) has previously been calculated by Pham et al. using a plane-wave basis set. 38,45The diagonal quasi-particle equation was solved as in the present work.Moreover, the analytic structure of the selfenergy was taken into account consistently.These results are hence fully comparable to the results of the present paper.For seven out of 12 of the molecules in the intersecting subset, we observe an agreement within ∼0.1 eV of our extrapolated results.The molecules containing fluorine show a somewhat larger discrepancy and so do ozone and pyridine (up to 0.5 eV).For ozone, the slope of the self-energy around the quasi-particle energy is steep, which implies a low quasi-particle weight (see section 4.1.1)and furthermore exhibits multisolution behavior as discussed previously and in section 4.1.1.In pyridine, the self-energy also exhibits a relatively large slope of 0.35.However, no multisolution behavior aggravates the solution of the quasi-particle equation.
DNA bases have been studied by Qian et al. 31 They used a plane-wave basis in combination with a Wannier-type optimized basis for the response function.The diagonal quasi-particle equation was solved as in the present work, and the self-energy was calculated using the analytic continuation.They estimated that their results are accurate to within 0.1 eV.Our extrapolated results for G, A, C, T lie approximately 0.2 eV lower in energy (0.4 eV for uracil).There are two possible causes for the observed deviation.First, the comparison of extrapolated results to none-extrapolated results, which tends to lower the extrapolated results.The second cause is the comparison of all-electron to froze-core results; increasing the number of valence electrons increases the screening.This reduces the absolute size of the (in general positive) contribution of the correlation part of the self-energy and hence lowers the QP energies.Both effects thus point in the direction of the observed deviation.

Comparison with Local-Orbital Results
. Recently, several local-orbital G 0 W 0 implementations have emerged.We have included the G 0 W 0 @PBE values reported with these codes 29,30,32,33,40,41,44,47−50 in Table 7.We often observe numerical differences that are larger than 0.1 eV and, for some systems, on the order of 1 eV.For example, the G 0 W 0 @ PBE ionization energy of the CO 2 molecule has already been studied in several other works. 29,40,47,49The reported values span a range of almost 1 eV (12.8−13.6 eV, the experimental ionization energy being 13.78 eV 122 ).For other systems, e.g., N 2 and NH 3 , the published G 0 W 0 @PBE ionization energies are also distributed in a similar energy range.
We attribute the large spread to the different basis function choices (Gaussian, molecular orbital, etc.) and sizes (triple-ζ, quadruple-ζ, etc.) and the fact that no extrapolation to the complete basis set limit has been performed.The different treatments of the frequency dependence also play a role.The results reported in this work allow us to estimate the impact of these approximations.For instance, Figure 10 illustrates that the error of a two-pole analytic continuation is mostly smaller than 0.12 eV.Conversely, Figure 6 indicates that the basis set size has a considerable impact on the ionization energies.Even for the relatively large QZVP basis, the QP-HOMO energies are underestimated by ∼0.12 eV on average.
Another factor that may affect the comparison with the published G 0 W 0 ionization energies is the way in which the quasi-particle equation is solved.The quasi-particle energies may be determined either from the iterative solution of eq 32 or, alternatively, from the linearized quasi-particle equation (see, for example, ref 136.)  that the ionization energies obtained from these two approaches may differ by up to 0.5 eV.
In summary, we believe that deviations up to a few hundred meV can be explained by the technical differences in the G 0 W 0 implementations discussed above.We mention that the results of Rostgaard et al. 29 deviate in 4 out of 22 cases by ∼1 eV or more (up to 2 eV), which seems hard to reconcile with our general findings.
4.4.Comparison to Experimental Values.In this work, the focus is on methodology; we establish a G 0 W 0 @PBE benchmark that holds (with certain bounds) irrespective of how the G 0 W 0 equations have been implemented.Residual discrepancies to experimental data persist, which are typically several hundred meV (see Table 8).Our work strongly suggests that on the theory side these discrepancies reflect (i) the neglect of self-consistency and (ii) vertex corrections.
For some classes of molecules, an informed choice of the starting point, e.g., B3LYP, PBE0, or static COHSEX instead of PBE, could lead to a reduction of the discrepancies to experiment by factors of order unity (see, for example, refs 43 and 137).However, corrections that are applicable to broad system classes and that can be systematically improved upon can probably not be achieved in this way.Our table also shows a comparison between the simplifying GPP model and experimental data.The comparison is favorable, at first sight, which at this point is not fully understood and a matter of ongoing investigation. 63

DISCUSSION
As discussed in section 3.2, both the KS energies and the exchange part of the self-energy are well converged.Moreover, we have already discussed in section 3.2 the errors introduced by the finite basis set and additional computational approximations, such as RI and the freezing of core states.In this section, we discuss the differences between the QP energies reported in Tables 2 and 5 that originate from the correlation part of the GW self-energy.We start by discussing the appearance of multiple solutions of the QP equation and their impact on the QP energies.Then, we discuss the aspects that relate to different approaches to treat the frequency dependence of the self-energy.5.1.Multiple Solutions of the QP Equation.For four molecules, ozone, boron nitride, beryllium oxide, and magnesium oxide, we found multiple solutions with comparable energy for the QP-HOMO level.To illustrate this behavior, we showed the correlation part of the G 0 W 0 @PBE self-energy for ozone in Figure 11.In Figure 13, we show that of magnesium oxide.Similar plots for boron nitride and beryllium oxide are provided in the Supporting Information.For boron nitride, Δ(ω) crosses one pole of Σ(ω), which results in three solutions.For ozone, beryllium oxide, and magnesium oxide, we obtain five solutions.
Inspection of Figures 11 and 13 confirms that the multisolution behavior typically emerges when the intersection is close to poles in the self-energy.For BrK, LiH, KH, copper cyanide, and Cu 2 , we also observe large self-energy derivatives (small Z-factors), but the poles are not close enough to the solution of the quasi-particle equation to give rise to multiple solutions.An example is show in Figure 14 for copper cyanide.The slope of the self-energy is large but rather constant between the KS energy and the QP energy.Any Newton− Raphson-like algorithm to iterate the QP equation will converge to a unique QP energy when started at the KS energy.This QP energy is, however, close to a pole, hence small differences in the self-energy calculated by the different approaches do result in changes in the final value.This behavior differs from the case of magnesium oxide, for example, where the first iteration step takes the algorithm to the other side of the first pole.
Figures 11 and 13, as well as Table 4, illustrate that the occurrence of multiple solutions strongly depends on the frequency treatment of the G 0 W 0 implementation.In general, the GPP approach and the two pole analytic continuation gives rise to a self-energy that has too few poles or poles that are considerably broadened.Additionally, the GPP approximation shifts the self-energy poles away from the quasi-particle energies, resulting in higher predicted G 0 W 0 quasi-particle energies.The result is a slowly varying self-energy with a unique solution to the quasi-particle equation in the energy region of interest.For many systems (including many of the molecules in GW100), the true poles of Σ are sufficiently far away from the quasi-particle energies or several poles merge into a broad pole with a gently sloped self-energy in the energy region of interest.In these cases, the GPP model and the two pole analytic continuation can be expected to agree well with the analytic and numerical full frequency methods (see also section 5.2).
In the case of multiple solutions, the solution that is reported by the code depends on the iterative solver that is employed for the QP equation.The solver in TURBOMOLE is constructed to always find the solution with the largest quasi-particle weight, i.e., the smallest derivative.In FHI-aims, a simple iteration is used, which sometimes yields a different solution.In BerkeleyGW, the QP equation is linearized in the case of GPP and solved by graphically finding the point where ϵ KS + Σ − V XC is equivalent to ϵ QP in the case of FF.
5.2.Pade, Two-Pole, and Plasmon Pole Approximations.For most systems, the solutions of the diagonal QP equation (eq 32) are found in a region in which the self-energy is rather flat, i.e., ∂Σ/∂ε is small and the quasi-particle weight is close to unity.In these cases, the self-energy is only evaluated in an energy region in which there are no poles nearby.For this purpose, a simple approximation for the energy dependence of Σ, such as the two-pole analytic continuation, may provide quasi-particle energies converged within 0.2 eV, as the data reported in Tables 2 and 5 indicate.This is further illustrated in Figure 15, where we report the difference between the AIMS-2P and TM-noRI ionization energies (left panel) and between BGW-GPP and EXTRA (right panel), as a function of Σ′ = ∂Σ/∂ε For small values of Σ′ (−0.5<Σ′ < 0), the deviation between AIMS-2P and TM-noRI is consistently below 0.2 eV.In this case, the two-pole analytic continuation tends to overestimate the ionization energies by ∼0.1 eV on average.However, for Σ′ < −0.5, AIMS-2P systematically underestimates the ionization energies by 0.4−0.5 eV as compared to TM-noRI. 138 Molecules that present large derivatives ∂Σ/∂ϵ for the QP-HOMO level are, for instance, BrK, LiH, KH, Cu 2 , NCCu, O 3 , BN, BeO, and MgO.Here, the quasi-particle equation has solutions close to a pole of the self-energy, as illustrated in Figures 11 and 13.As a consequence, the resolution of the pole-structure of the self-energy, as achieved by different methods, now matters more.For O 3 , BeO, and MgO, AIMS-P16 provides enough fitting parameters.Conversely, for BN and CNCu, it is necessary to increase the number of Padeṕ arameters up to 128 to describe the poles of the self-energy sufficiently well, as illustrated in Figures 11 and 13 and Table 4.
Very smooth self-energies result from plasmon-pole models for the response function.Besides the "flattening" of the correlation part of the self-energy, we also observe an overall decreased correlation part of the self-energy (see Figures 11  and 13) that is not present in the 2-pole self-energy.
5.3.Numerical Full Frequency Integration.In contrast to imaginary frequency treatments and plasmon pole models, a full frequency treatment can in principle be converged to the exact analytic G 0 W 0 self-energy by increasing the density of the frequency grid.(The Padéapproximation could, in principle from a mathematical point of view, also be converged fully to the analytic self-energy.However, in practice, it is very hard to fit more than 64 poles to the smoothly varying self-energy on the imaginary axis.)We indeed see that the FF results reproduce the pole-positions in the analytically treated selfenergy accurately (see Figures 11 and 13).Even in MgO, which features many low intensity poles, we observe this agreement for many poles positions.In contrast, the pole-broadening in FF is much more pronounced.This is caused by the finite mesh size of the numerical integration, Δω, and the imaginary shift η (see eq 18).The full frequency approach ideally requires Δω ≪ η.In practice, however, Δω ≲ η, which introduces considerable additional broadening.To simulate this effect in TURBOMOLE, we ran calculations with η = 0.2 eV and obtained self-energy matrix elements almost identical with FF (not shown).
The agreement of the pole positions and the self-energy curves calculated at equal broadening η suggests that both the plane-wave basis and the localized basis would, in a fully converged calculation, produce the same results.Decreasing η in eq 18 would, however, necessitate a smaller Δω, which would increase the number of frequency points to a computationally prohibitive amount.The deviations seen in Figure 16 are therefore caused mostly by dicretization errors.
In principle, the accuracy of the FF approach can systematically be improved by decreasing Δω and η, respecting the condition Δω≪ η.However, fully converging the FF approach to reproduce the complex pole structure of Σ that arises in molecular systems may often be computationally prohibitive for actual production calculations.

CONCLUSIONS
In this work, we have benchmarked the G 0 W 0 method for molecular systems.We have devised the GW100 benchmark set of 100 closed shell molecules containing a variety of elements and chemical bonding types.Using G 0 W 0 @PBE, we report converged QP-HOMO and QP-LUMO energies that can serve as reference for future GW development work.For completeness, we also included the experimental ionization energies and affinities where available.
We have compared the three G 0 W 0 implementations in FHIaims, BerkeleyGW, and TURBOMOLE, as well as several approximations and numerical approaches that are available in these codes.
At convergence, the QP-HOMO and QP-LUMO levels agree on the order of 200 meV.AIMS-P16 and TM-noRI even agree to a precision of 1 meV.
In the process of this work, we have identified two crucial aspects that control the accuracy of the G 0 W 0 quasi-particle energies: the size of the basis set and the treatment of the frequency dependence.The complex pole structure of the response function of molecules cannot accurately be described by commonly applied approximate approaches, such as plasmon-pole models or an analytic continuation with only 2 poles.In contrast, an analytic continuation with high order Padéapproximants and numerical full frequency treatments facilitates an accurate ab initio description but care needs to be to taken to converge the numerical parameters (the number of pole parameters and the frequency grid on the imaginary axis in the Padéapproach and the frequency grid in the FF method).
With the derivative of the self-energy, we have also identified a rule of thumb criterion that allows us to detect problematic cases that may require special attention.As long as the derivative at the QP energies stays small, |∂Σ/∂ ω| ≪ 1, there are no poles close by and the detailed pole structure of the selfenergy is of lesser importance.If the absolute value of the derivative becomes of order unity, poles in the self-energy are   situated closer to the quasi-particle energies.This could be an indication for the emergence of multiple solutions.As it turns out, our test set includes four such molecules.For these molecules, the different solvers for the iterative solution of the quasi-particle equation that are implemented in the three codes converge on different solutions.

Figure 2 .
Figure 2. Convergence behavior of the BGW-GPP ionization energy of ozone and magnesium oxide.The left panel shows convergence with respect to the dielectric function cutoff with the number of unoccupied states fixed such that it spans 90 eV above the vacuum level, and the right panel shows the convergence with respect to number of states for a dielectric function cutoff of 408 eV.The energy is given as a difference with respect to the highest convergence criteria.

Figure 3 .
Figure 3. Number of molecules for which the estimated convergence error associated with the chosen E max and E ϵ cut is within a given interval (given in eV).The estimated error is the difference between the highest convergence parameter tested and the parameter used for all 100 molecules.

Figure 4 .
Figure 4. Histogram of the deviation of the KS-HOMO energies in the GW100 set calculated with the def2-SVP (left), def2-TZVP (center), and def2-QZVP (right) basis sets from the extrapolated complete basis set limit using TURBOMOLE.The raw data can be found in the Supporting Information.The extrapolated results are obtained by linear extrapolation of the def2-TZVP and def2-QZVP results as a function of the inverse of the number of basis functions.

Figure 5 .
Figure 5. Histogram of the deviation of the QP-HOMO and QP-LUMO energies of the GW100 set calculated with the def2-SVP (left), def2-TZVP (center), and def2-QZVP (right) basis from the extrapolated complete basis set limit using TURBOMOLE.The extrapolated results are obtained by linear extrapolation of the def2-TZVP and def2-QZVP results as a function of the inverse of the number of basis functions.

Figure 6 .
Figure 6.Deviation of the QP-HOMO energies from the extrapolated complete basis set limit for the def2-SVP, def2-TZVP and def2-QZVP basis sets as a function of the inverse number of basis functions.The results have been obtained with TURBOMOLE.

Figure 7 .
Figure 7. G 0 W 0 RI error of the QP-HOMO energies in the TURBOMOLE calculations compared to the reciprocal of the number of basis functions (1/NBF).The RI error for the G 0 W 0 calculations is defined as the difference between the QP-energies calculated with and without applying the RI approach at all steps of the calculation.

Figure 9 .
Figure 9. Schematic of a graphical solution of the quasi-particle equation (eq 32).All intersections of the red line with the correlation part of the selfenergy (black line) are solutions of the quasi-particle equation.The left panel shows the most common situation with a clear single solution near the KS starting point; the right panel shows the situation that can lead to multiple solutions.

Figure 10 .
Figure 10.Absolute difference between QP-HOMO energies calculated with the different G 0 W 0 implementations.The FHI-aims and TURBOMOLE results are compared at the QZVP level.The BerkeleyGW results are compared to the extrapolated values.The averages of the absolute deviations are shown as horizontal lines.

Figure 11 .
Figure 11.Comparison of the energy-dependent correlation part of the self-energy Σ c (ε) calculated with the three different codes using different procedures for ozone."TM-5" and "TM-3" indicate TURBOMOLE results calculated with imaginary shifts of η = 10 −3 and 10 −5 Hartree, respectively.The intersections of these curves with the (red) line ω − ϵ KS + V xc − Σ x correspond to the solutions of the QP equation.

Figure 12 .
Figure 12.Absolute difference between QP-LUMO energies calculated with the different G 0 W 0 implementations.The FHI-aims and TURBOMOLE results are compared at the QZVP level.The BerkeleyGW results are compared to the extrapolated values.The averages of the absolute deviations are shown as horizontal lines.

Figure 13 .
Figure 13.Comparison of the energy-dependent correlation part of the self-energy Σ c (ε) calculated with the three different codes using different procedures for magnesium oxide."TM-5" and "TM-3" indicate TURBOMOLE results calculated with imaginary shifts of η = 1e−3 and 1e−5 Hartree, respectively.The intersections of these curves with the (red) line ω − ϵ KS + V xc − Σ x correspond to the solutions of the QP equation.

Figure 14 .
Figure 14.Energy-dependent correlation part of the self-energy Σ c (ε) of copper cyanide calculated using TURBOMOLE.

Figure 16 .
Figure 16.Deviation between the extrapolated quasi-particle energies and BerkeleyGW evaluated in the GPP (left) and with a full frequency treatment (right).

Table 1 .
Computational Parameters for 18 Molecules Computed within FF-BGW a a Columns 2, 3, and 5 are fixed at the converged value from the BGW-GPP calculation.b Number of valence bands c Number of conduction bands

Table 2 .
Ionization Potentials (the Negative of the QP-HOMO Energies) of the GW100 Set Calculated with G 0 W 0 @PBE using TURBOMOLE, FHI-aims, and BerkeleyGW a

Table 4 .
. Looking at our CO 2 example, it is linearized in refs 29, 47, and 49 and solved in ref40.We find Solutions of the Quasi-Particle Equation for the QP-HOMO Level of Ozone, Boron Nitride, Beryllium Oxide, and Magnesium Oxide for the Different Calculations a The numbers in boldface are the solutions as reported by the specific code (as reported solution).In Table2, we always report the solution with the highest energy. a

Table 5 .
LUMO Energies Analagous to Table2 a

Table 5 . continued
The experimental electron affinities marked with * are laser photoelectron spectroscopy values.We note that the QP-LUMOs of the single atoms in the GW100 set (He−Xe) exhibit such a strong dependence on the basis set that a controlled extrapolation is not possible within the SVP, TZVP, or QZVP series.As before, no extrapolated values are presented for molecules containing fifth row elements because SVP and TZVP all-electron basis sets are not available. a

Table 6 .
Mean Deviation (MD) and Mean Absolute Deviation (MAD) between the QP-LUMO Energies Presented in Table2

Table 7 .
G 0 W 0 @PBE Literature Values for QP-HOMO Energies Compared to the BGW-GPP and Extrapolated Data from this Study a a Using either a local-orbital (LO) or plane-wave (PW) basis.

Table 8 .
Mean Deviation (MD) and Mean Absolute Deviation (MAD) between the QP-HOMO Energies presented in Table 2 and the Experimental Ionization Energies with the Second Half of the Table Comparing Only Those Systems for which the Experimental Ionization Energy is Vertical