Practical GW scheme for electronic structure of 3 d -transition-metal monoxide anions: ScO − , TiO − , CuO − , and ZnO −
PPractical GW scheme for electronic structure of 3 d -transition-metal monoxide anions:ScO − , TiO − , CuO − , and ZnO − Young-Moo Byun and Serdar ¨O˘g¨ut ∗ Department of Physics, University of Illinois at Chicago, Chicago, IL 60607, USA (Dated: September 4, 2019)The GW approximation to many-body perturbation theory is a reliable tool for describing chargedelectronic excitations, and it has been successfully applied to a wide range of extended systems forseveral decades using a plane-wave basis. However, the GW approximation has been used to testlimited spectral properties of a limited set of finite systems (e.g. frontier orbital energies of closed-shell sp molecules) only for about a decade using a local-orbital basis. Here, we calculate thequasiparticle spectra of closed- and open-shell molecular anions with partially and completely filled3 d shells (shallow and deep 3 d states, respectively), ScO − , TiO − , CuO − , and ZnO − , using variouslevels of GW theory, and compare them to experiments to evaluate the performance of the GW approximation on the electronic structure of small molecules containing 3 d transition metals. We findthat the G -only eigenvalue self-consistent GW scheme with W fixed to the PBE level ( G n W @PBE),which gives the best compromise between accuracy and efficiency for solids, also gives good resultsfor both localized ( d ) and delocalized ( sp ) states of 3 d -transition-metal oxide molecules. The successof G n W @PBE in predicting electronic excitations in these systems reasonably well is likely due tothe fortuitous cancellation effect between the overscreening of the Coulomb interaction by PBEand the underscreening by the neglect of vertex corrections. Together with the absence of theself-consistent field convergence error (e.g. spin contamination in open-shell systems) and the GW multi-solution issue, the G n W @PBE scheme gives the possibility to predict the electronic structureof complex real systems (e.g. molecule-solid and sp - d hybrid systems) accurately and efficiently. I. INTRODUCTION
It is a challenging task to accurately determine theelectronic structure of an interacting many-electron sys-tem. In experiment, electron removal and addition en-ergies of both extended and finite systems are measuredby direct and inverse photoelectron spectroscopy (PESand IPES, respectively). In theory, it is well known thatthe GW approximation to many-body perturbation the-ory (MBPT) describes bandgaps and band structures ofsolids more accurately than local and semi-local approxi-mations to density-functional theory (DFT). However,less is known about the performance of the GW approx-imation on the electronic structure of atoms, molecules,and clusters. Especially, GW calculations for the quasi-particle (QP) spectra of open-shell molecules containing3 d transition metals are scarce. There are a few reasonsfor it.First, it is easier to test only frontier orbital energiessuch as the ionization energy (IE) and the electron affin-ity (EA) than the full QP spectrum (all orbital energies).There are mainly two ways to calculate IE and EA ofmolecules. On one hand, IE (EA) can be obtained fromDFT, HF (Hartree-Fock), MP2 (second-order Møller-Plesset perturbation theory), RPA (random-phase ap-proximation), or CCSD(T) (coupled-cluster singles anddoubles plus perturbative triples) total energy differencesbetween a neutral and a cation (anion) within the so-called ∆SCF (self-consistent field) method. Generally,the ∆SCF method gives good results for frontier or-bital energies of finite systems, but cannot be appliedto extended systems. Also, it is not straightforward forthe ∆SCF method to access the full QP spectrum. On the other hand, IE and EA can be obtained from GW eigenvalues for the HOMO (highest occupied molecularorbital) and the LUMO (lowest occupied molecular or-bital), respectively. Due to the simplicity of the ∆SCFmethod, many studies have evaluated the performance ofthe GW approximation on molecules by comparing GW IE and EA to ∆SCF ones, but that approach does notutilize the full power of the GW approximation, which isthe ability to provide the QP spectrum for both finite andextended systems. For example, Bethe-Salpeter equation(BSE) calculations for optical excitations require moreorbital energies than IE and EA as input. Second, it is easier to test closed-shell systems thanopen-shell ones. Most of quantum chemistry-based GW implementations for finite systems, such as MOLGW, FIESTA, TURBOMOLE, FHI-AIMS, and CP2K, use local-orbital basis sets such as Gaussian basis sets. GW calculations require mean-field self-consistent calcu-lations, such as restricted and unresticted Hartree-Fockor Kohn-Sham (RHF or RKS and UHF or UKS, respec-tively) calculations for closed- and open-shell systems,respectively. The problem is that unlike RHF and RKSself-consistent calculations, UHF and UKS ones are notguaranteed to converge, as their convergence stronglydepends on the initial guess wavefunctions. This is es-pecially the case for spin-unrestricted calculations per-formed with hybrid exchange-correlation (xc) function-als, which include a fraction of exact exchange (EXX),and HF on 3 d -transition-metal-containing molecules dueto the near-degeneracy of energy levels. Partially dueto this SCF convergence issue, most existing studies haveused only closed-shell systems to assess the performanceof the GW approximation on finite systems. For exam- a r X i v : . [ c ond - m a t . m t r l - s c i ] S e p ple, Refs. 13–16 used the so-called GW
100 benchmarkset, which is composed of only closed-shell molecules.Last, it is easier to test sp -electron systems than d -electron ones. Fundamentally, it is more difficult to accu-rately predict the electronic structure of d systems (espe-cially, 3 d systems) than sp ones because of the strong lo-calization, and thus the strong correlation, of d electrons.For example, it is challenging for GW to accurately re-produce the experimental bandgap and d -band positionof bulk ZnO at the same time. Practically, it is com-putationally more demanding to tackle systems with d -electrons than those with only sp -electrons. For example, d elements have more basis functions than sp ones, whichincreases the computational effort, and transition-metal-containing molecules, especially with partially filled d shells and low multiplicity states, aggravate the above-mentioned SCF convergence issue, which increases thehuman effort by making it necessary to manually exploremany minima with similar energies using many initialguess wavefunctions. The GW approximation is unique, but due to its highcomputational costs, there are various GW schemes andvariants. Generally, there are two approaches. One ap-proach is to vary the self-consistency level in the GW approximation. The GW self-consistent levels from thelowest to the highest include the perturbative non-self-consistent (one-shot) GW ( G W ) scheme, the eigen-value self-consistent GW (ev GW ) scheme (with twotypes G n W and G n W n , which update eigenvalues onlyin G and in both G and W , respectively), the QP self-consistent GW (QS GW ) scheme using a static and Her-mitian approximation to the GW self-energy, and thefully self-consistent GW (SC GW ) scheme. Generally, asthe GW self-consistency level increases, the GW approx-imation depends less on the mean-field starting point andbecomes more conserving with respect to particle num-ber, momentum, and energy. However, the higher GW self-consistency level does not necessarily give more ac-curate QP energies because vertex corrections are miss-ing in the GW approximation. For example, SC GW and QS GW systematically overestimate the bandgapsof solids, displaying worse performance than ev GW ,which currently provides the best balance between accu-racy and efficiency for solids. The other approach is to vary the amount of EXX inthe GW starting point to reduce the self-interaction er-ror by (semi-)local xc functionals. Typically, the G W scheme chooses this approach to obtain good resultsat low computational costs. However, the predictivepower of this approach is questionable, since the optimalamount of EXX in the GW starting point is stronglysystem-dependent. For example, for extended systems,the reported values for the optimal amount of EXX arenarrowly spread between 0% and 25%, , while for fi-nite systems, they are widely spread between 25% and100%. The purpose of this work is to evaluate the perfor-mance of the GW approximation on the electronic struc- ture of small oxide molecules containing 3 d transitionmetals. To this end, we calculate the QP spectra ofclosed- and open-shell molecular anions with partiallyand completely filled 3 d shells, ScO − , TiO − , CuO − , andZnO − , using various levels of GW theory. There are afew reasons why we chose these molecular systems: (i)their anion PES data is available, (ii) CuO − andZnO − are molecular analogs to bulk Cu O and ZnO,respectively, which are challenging systems for the GW method, and (iii) shallow and deep 3 d states are mea-sured in TiO − and CuO − , respectively.This article is organized as follows: First, we givea brief introduction to the GW approximation and itsimplementation in the framework of quantum chem-istry. Second, we present various convergence test re-sults and show that care should be taken to obtain reli-able and reproducible QP energies of finite systems fromGaussian-based GW implementations. Third, we as-sess various GW schemes, focusing on ionization energiesand 3 d -electron binding energies, and conclude that the G n W @PBE scheme gives the best performance among GW schemes considered in this work in terms of accuracyand efficiency. Last, we discuss the origin of seeminglyconflicting GW results for finite systems in the literature. II. THEORETICAL BACKGROUND
In this section, we briefly review the GW approxima-tion and its implementation using local-orbital basis sets.This section contains only a minimal number of equa-tions, which will be needed later. More details can befound in Refs. 6, 8, 13, and 30. Generally, we followthe notation in the MOLGW implementation paper forconsistency: (i) Hartree atomic units are used in all equa-tions, (ii) The complex conjugate notation is not used forwavefunctions, because they are real in finite systems,(iii) State indices i and j run over only occupied states, a and b run over only empty (virtual) states, and m and n run over all states, (iv) The response function is referedto as the polarizability instead of the susceptibility, and(v) χ is used for the polarizability instead of P and Π. A. GW Approximation
In Hedin’s GW approximation, the non-local, dynam-ical, and non-Hermitian self-energy Σ σ at frequency ω isgiven byΣ σ ( r , r (cid:48) , ω ) = i π (cid:90) dω (cid:48) e iηω (cid:48) G σ ( r , r (cid:48) , ω + ω (cid:48) ) W ( r (cid:48) , r , ω (cid:48) ) , (1)where σ is the spin channel ( ↑ or ↓ ), G σ is the time-ordered one-particle Green’s function, W is the dynam-ically screened Coulomb interaction, and η is a positiveinfinitesimal.The self-energy in Eq. (1) can be calculated from firstprinciples by solving the coupled Hedin’s equations in or-der. One starts by constructing the one-particle Green’sfunction using the one-electron eigenvalues (cid:15) σm and corre-sponding wavefunctions ϕ σm ( r ) obtained from the Hartreeor mean-field approximation: G σ ( r , r (cid:48) , ω ) = (cid:88) i ϕ σi ( r ) ϕ σi ( r (cid:48) ) ω − (cid:15) σi − iη + (cid:88) a ϕ σa ( r ) ϕ σa ( r (cid:48) ) ω − (cid:15) σa + iη , (2)where i runs over occupied states and a runs over emptystates. Note that G σ in Eq. (2) is not the interact-ing (dressed) Green’s function, but the non-interacting(bare) one, which are conventionally denoted by G σ and G σ , respectively. In this work, we use the subscript 0 todistinguish the non-self-consistent GW scheme from theself-consistent one.Using the one-particle Green’s function in Eq. (2), onecan successively obtain the non-interacting (irreducible)polarizability χ and the interacting (reducible) polariz-ability χ = χ [1 − vχ ] − within the RPA, the screenedCoulomb interaction, and the self-energy: χ = − i (cid:88) σ G σ G σ , (3) W = v + vχ W = v + vχ v + vχ vχ v + ... = v + vχv, (4)Σ σ = iG σ W = iG σ ( v + vχv ) = Σ σ x + Σ σ c , (5)where v denotes the bare (unscreened) Coulomb inter-action v ( r , r (cid:48) ) = 1 / | r − r (cid:48) | , Σ x is the exchange part ofthe self-energy, and Σ c is the correlation part of the self-energy. Note that in Eqs. (3), (4), and (5), space and fre-quency variables ( r , r (cid:48) , ω ) are omitted for simplicity, and χ ( ω ), χ ( ω ), Σ σ ( ω ), and Σ σ c ( ω ) are dynamic, whereas v and Σ σ x are static. Note also that W is obtained withoutusing the dielectric matrix.Using the real part (Re) of the self-energy in Eq. (5)and the first-order perturbation theory, one can obtainthe (diagonal) QP equation: (cid:15) G W ,σm = (cid:15) σm + (cid:104) ϕ σm | ReΣ σ ( (cid:15) G W ,σm ) − v σ xc | ϕ σm (cid:105) , (6)where (cid:15) G W ,σm are the perturbative one-shot GW QPenergies and v σ xc is the xc potential. Experimentally, (cid:15) G W ,σm correspond to vertical IEs and EAs in PES andIPES, respectively. Theoretically, (cid:15) G W ,σm correspond tothe positions of poles of the Green’s function in the spec-tral (Lehmann) representation and thereby to the posi-tions of QP peaks and plasmon satellites in the corre-sponding spectral function A σ : A σmm ( r , r (cid:48) , ω ) = 1 π | Im G σmm ( r , r (cid:48) , ω ) | , (7)where A σmm are the diagonal elements of the spectralfunction, G σmm are the diagonal elements of the Green’s function, and Im represents the imaginary part. Notethat A σmm give the local density of states. G σ in Eq. (7) is the interacting Green’s function,whereas G σ in Eq. (2) is the non-interacting one. In otherwords, by plugging G σ in Eq. (2) into Eq. (7) after replac-ing (cid:15) σm , where m = i or a , by (cid:15) σm + (cid:104) ϕ σm | Σ σ ( ω ) − v σ xc | ϕ σm (cid:105) ,one can find that A σmm have Lorentzian peaks at ω = (cid:15) σm + (cid:104) ϕ σm | ReΣ σ ( ω ) − v σ xc | ϕ σm (cid:105) , (8)which shows that solving the QP equation in Eq. (6) andlocating the peak positions in the spectral function inEq. (7) are equivalent ways of obtaining the QP energies.The QP equation in Eq. (6) is non-linear, because Σ σ depends on (cid:15) G W ,σm , so it should be solved numerically.Additionally, Hedin equations are coupled, because W and Σ σ depend on G σ , so they should be solved self-consistently. Multiple ways to numerically solve the non-linear QP equation and to iteratively solve the coupledHedin equations will be discussed later. B. Self-Consistent Field Method
In order to obtain the ingredients for the one-particleGreen’s function in Eq. (2) using local-orbital basis sets,molecular orbitals (MOs) and corresponing MO energiesare used as one-electron wavefunctions and correspondingeigenvalues. MOs are expanded as a linear combinationof atomic orbitals (AOs) φ µ : ϕ σm ( r ) = (cid:88) µ C σµm φ µ ( r ) , (9)where C σµm are MO expansion coefficients. In MOLGW,atom-centered (contracted) Gaussian orbitals are used asAOs.The MO coefficients in Eq. (9) and MO energiescan be obtained by solving the generalized Kohn-Sham(gKS) equation [the Hartree-Fock–Kohn-Sham scheme for (semi-)local functionals, hybrid functionals, and HF]: H σ C σ = SC σ (cid:15) σ , (10)where C σ is a matrix of MO coefficients, (cid:15) σ is a diagonalmatrix of MO energies, S is the AO overlap matrix withelements: S µν = (cid:90) d r φ µ ( r ) φ ν ( r ) , (11)and H σ is the Hamiltonian matrix with elements: H σµν = T µν + V ext ,µν + J µν − αK σµν +(1 − α ) V PBE ,σ x ,µν + V PBE ,σ c ,µν , (12)where T , V ext , J , and K σ are the kinetic energy, externalpotential energy, Hartree, and Fock exchange terms, re-spectively, V σ x and V σ c are the exchange and correlationpotentials, respectively, and α is the fraction of EXX inhybrid functionals that will be introduced later.We briefly explain only a few terms in the Hamiltonianmatrix in Eq. (12), which will be needed later. The ma-trix elements of the Hartree term in Eq. (12) are givenby J µν = (cid:88) λτ ( µν | λτ ) (cid:88) σ D σλτ , (13)where ( µν | λτ ) are the 4-center two-electron Coulomb re-pulsion integrals:( µν | λτ ) = (cid:90) (cid:90) d r d r (cid:48) φ µ ( r ) φ ν ( r ) 1 | r − r (cid:48) | φ λ ( r (cid:48) ) φ τ ( r (cid:48) ) , (14)and D σ is the density matrix with elements: D σµν = (cid:88) m f σm C σµm C σνm , (15)where f σ is the occupation number (0 or 1). The matrixelements of the Fock exchange term in Eq. (12) are givenby K σµν = (cid:88) λτ D σλτ ( µλ | τ ν ) . (16)The exchange and correlation potentials in Eq. (12) de-pend on the density ρ σ (and the density gradient ∇ ρ σ ): ρ σ ( r ) = (cid:88) µν D σµν φ µ ( r ) φ ν ( r ) . (17)The gKS equation in Eq. (10) (the restrictedRoothaan-Hall or unrestricted Pople-Nesbet equations)should be solved using the SCF method, because J , K σ , V σ x , and V σ c in Eq. (12) depend on the density matrix inEq. (15), as shown in Eqs. (13), (16), and (17). C. GW Self-Energy
In order to obtain the ingredients for the interactingpolarizability in Eq. (4), one should solve the Casidaequation in matrix form: (cid:18)
A B − B − A (cid:19) (cid:18) X s Y s (cid:19) = Ω s (cid:18) X s Y s (cid:19) , (18)where A and B are the resonant and coupling matri-ces, respectively, and Ω s and ( X s , Y s ) are the eivenvalues(the neutral two-particle excitation energies) and corre-sponding eigenvectors, respectively. The matrix elementsin A and B are given by A jbσ (cid:48) iaσ = ( (cid:15) σa − (cid:15) σi ) δ ij δ ab δ σσ (cid:48) + ( iaσ | jbσ (cid:48) ) + f jbσ (cid:48) xc ,iaσ , (19) B jbσ (cid:48) iaσ = ( iaσ | bjσ (cid:48) ) + f bjσ (cid:48) xc ,iaσ , (20)where i and j are for occupied states, a and b arefor empty states, f xc is the time-dependent density-functional theory (TDDFT) xc kernel, and ( iaσ | jbσ (cid:48) ) are the 4-orbital two-electron Coulomb repulsion integrals:( iaσ | jbσ (cid:48) ) = (cid:90) (cid:90) d r d r (cid:48) ϕ σi ( r ) ϕ σa ( r ) 1 | r − r (cid:48) | ϕ σ (cid:48) j ( r (cid:48) ) ϕ σ (cid:48) b ( r (cid:48) ) . (21)In this work, we used the RPA by setting f xc = 0. Notethat Ref. 8 showed that TDDFT and RPA polarizabilitiesmake a difference of ∼ G W @PBE QP HOMOenergy. Note also that the dimension of the Casida ma-trix in Eq. (18) scales as O ( N ) with N being the sys-tem size, so building and completely diagonalizing theCasida matrix scale as O ( N ) and O ( N ), respectively.The MO integrals in Eq. (21) are transformed from theAO integrals in Eq. (14) through the AO-MO integraltransformation:( iaσ | jbσ (cid:48) ) = (cid:88) µνλτ C σµi C σνa C σ (cid:48) λj C σ (cid:48) τb ( µν | λτ ) , (22)which scales as O ( N ). Note that this integral transfor-mation is a bottleneck in Gaussian-based GW and MP2calculations.Diagonalizing the Casida matrix in Eq. (18) yieldseigenvalues Ω s and eigenvectors ( X s , Y s ). In MOLGW,the diagonalization is performed without using theTamm-Dancoff approximation (TDA), which sets B tozero, but efficiently using the so-called beyond-TDAmethod. . Using Ω s and ( X s , Y s ), one can constructthe spectral representation of the interacting polarizabil-ity χ ( ω ). From χ ( ω ) and Eq. (4), one can obtain thespectral representation of the screened Coulomb interac-tion W ( ω ). Using W ( ω ) and Eq. (5), and analyt-ically performing the convolution of G σ ( ω ) and W ( ω ) inthe frequency domain, one can obtain the exchange andcorrelation parts of the GW self-energy Σ σ x and Σ σ c ( ω ),respectively, whose diagonal matrix elements are givenby (cid:104) ϕ σm | Σ σ x | ϕ σm (cid:105) = − (cid:88) i ( miσ | imσ ) , (23) (cid:104) ϕ σm | Σ σ c ( ω ) | ϕ σm (cid:105) = (cid:88) is w smiσ w smiσ ω − (cid:15) σi + Ω s − iη , + (cid:88) as w smaσ w smaσ ω − (cid:15) σa − Ω s + iη , (24)where i runs over occupied states, a runs over emptystates, s runs over all excitations, and w smnσ are given by w smnσ = (cid:88) iaσ (cid:48) ( mnσ | iaσ (cid:48) )( X siaσ (cid:48) + Y siaσ (cid:48) ) . (25)Note that unlike the plasmon-pole approximation(PPA), the analytic continuation method, and thecontour deformation technique, the fully analyticmethod employed in RGWBS, TURBOMOLE, andMOLGW gives the exact GW self-energy at all frequencypoints because it does not rely on any approximation andnumerical parameter. D. Spin Contamination
In unresticted HF and KS calculations for open-shellsystems, the expectation value of the total angular mo-mentum (cid:104) S (cid:105) is given by (cid:104) S (cid:105) = S ( S + 1) + N ↓ − N ↑ (cid:88) i N ↓ (cid:88) j |(cid:104) ϕ ↑ i | ϕ ↓ j (cid:105)| , (26)where N ↑ and N ↓ are the numbers of ↑ - and ↓ -spin elec-trons, respectively, and S is ( N ↑ − N ↓ )/2 with N ↑ > N ↓ .The last two terms on the right side of Eq. (26) are calledthe spin contamination, which is non-negative. Thespin contamination becomes large when a ground stateis mixed with (contaminated by) excited states.In restricted calculations for closed-shell systems, theSCF cycle always converges to a global minimum and thespin contamination is zero for all (semi-)local and hybridfunctionals as well as HF. In unrestricted calculations foropen-shell systems, the SCF convergence and the spincontamination depend on EXX amount and basis size.For (semi-)local functionals, the SCF cycle almost alwaysconverges to a global minimum and the spin contamina-tion is small [generally smaller than ∼
10% of S ( S + 1)].For hybrid functionals and HF, there is a chance (whichincreases with EXX amount and basis size) that the SCFcycle fails, does not converge, or converges to local min-ima or the spin contamination is large.There are a few points to note about the spin contami-nation. First, the spin contamination is just an indicatorfor the SCF convergence error, therefore, a small spincontamination does not guarantee the correct SCF con-vergence. Second, that the spin contamination generallyraises, but sometimes lowers the gKS total energy, so thelowest gKS total energy does not guarantee the correctSCF convergence, either. Last, the spin contaminationand the SCF cycle are independent of each other. Forexample, the SCF cycle can converge quickly with largespin contamination or slowly with small spin contamina-tion. E. Auxiliary Basis Sets and Multi-threadParallelization
In Gaussian-based GW , the 4-center integrals ( µν | λτ )in Eq. (14), which scale as O ( N ), are a common bottle-neck in gKS and GW parts in terms of compute timeand memory storage. One way to reduce the bottle-neck is the resolution-of-identity (RI) approximation (thedensity-fitting approximation), which expands the prod-uct of basis functions φ µ ( r ) φ ν ( r ) as a linear combinationof auxiliary basis functions φ P ( r ). There are twotypes of the RI approximation: RI-V using a Coulombmetric and RI-SVS using an overlap metric. For example,FIESTA uses both RI-V and RI-SVS, whereas MOLGWuses only RI-V, which is known to be superior to RI-SVS. Within RI-V, the 4-center integrals ( µν | λτ ) in Eq. (14)approximate to( µν | λτ ) ≈ (cid:88) P Q ( µν | P )( P | Q ) − ( Q | λτ ) , (27)where P and Q run over auxiliary basis functions, ( µν | P )and ( Q | λτ ) are the 3-center integrals, and ( P | Q ) are the2-center integrals.RI can be applied to both gKS [ J and K σ in Eqs. (13)and (16)] and GW [ A , B , Σ σ x , and Σ σ c ( ω ) in Eqs. (19),(20), (22), (23), (24), and (25)] parts. In this work, werefer to RI applied to one (both) of them as a half (full)RI method. For example, FIESTA uses a half RI method,whereas MOLGW uses a full RI method. In this work,we observed that a full RI method in MOLGW reducesboth compute time and memory storage by about thenumber of basis functions (by ∼
100 times as shown inTable I).RI is an approximation, so it causes an error. Thereare mixed results for the RI error in the literature, rang-ing from ∼ ∼ GW ), xc functionals (PBE vs HF), and basis sets areused to evaluate the quality of RI. Another way to reduce the bottleneck without caus-ing an error is the parallelization. We parallelizedthe 4-center integrals in Eqs. (13), (16), (19), (20),(22), (23), (24), and (25) [as well as other bottlenecks,such as the integral transformation in Eq. (22) andthe correlation part of the self-energy in Eq. (24)] us-ing Open Multi-Processing (OpenMP), which consumesmuch less memory than Message Passing Interface byusing shared-memory threads. The performance gain byour OpenMP parallelization is shown in supplementarymaterial. We also optimized our OpenMP implementa-tion to reduce Non-Uniform Memory Access (NUMA)effects in modern multi-core processors by enhancing thememory bandwidth and reducing the memory latency.Our OpenMP implementation in MOLGW 1.F has re-cently been merged into MOLGW 2.A. F. G W Quasiparticle Energy
In this work, we used three methods to obtain (cid:15) G W m ,as it is practically impossible to obtain unique, correct,and accurate (cid:15) G W m for all energy levels of all molecularsystems using a single method, which will be discussed indetail later. Note that in the following, the spin channel σ and the real part Re are omitted for simplicity.The first method is to linearize the non-linear QP equa-tion in Eq. (6): (cid:15) G W m ≈ (cid:15) m + Z linear m (cid:104) ϕ m | Σ( (cid:15) m ) − v xc | ϕ m (cid:105) ≡ (cid:15) G W , linear m , (28)where (cid:15) G W , linear m is the perturbative one-shot QP energyobtained from the linearization, and Z linear m is the QPrenormalization factor for the linearization: Z linear m ( (cid:15) m ) = 11 − ∂∂ω (cid:104) ϕ m | Σ( ω ) | ϕ m (cid:105)| ω = (cid:15) m , (29)where the derivative of the self-energy is obtained fromthe finite difference method using two frequency pointsat (cid:15) m ± ∆ ω with ∆ ω being the frequency grid spacing,which is set to 0.001 Ha in this work.There are a few points to note about the linearizationmethod. First, one can choose different frequency pointsfor the finite difference method (e.g. (cid:15) m ± . (cid:15) m ± . G W method, different frequency points give similar re-sults for (cid:15) G W , linear m , because PPA makes (cid:104) ϕ m | Σ c ( ω ) | ϕ m (cid:105) in Eq. (24) smooth around (cid:15) m by drastically reducingthe number of self-energy poles. However, in the full-frequency G W method, different frequency points cangive very different results for (cid:15) G W , linear m when the fi-nite difference method fails due to weak self-energy polesaround (cid:15) m , which will be discussed later. Second, thederivative of the self-energy can be evaluated analyti-cally, but it is not used in this work. Last, Z linear m inEq. (29) is slightly different from that in Ref. 13, Z m ( (cid:15) G W m ) = 11 − ∂∂ω (cid:104) ϕ m | Σ( ω ) | ϕ m (cid:105)| ω = (cid:15) G0W0 m . (30)The derivative of the self-energy is evalutated at ω = (cid:15) m and ω = (cid:15) G W m in Z linear m and Z m , respectively. Gener-ally, Z m is smaller than Z linear m because the self-energyhas a steeper slope at ω = (cid:15) G W m than at ω = (cid:15) m . Z m represents the QP weight (the pole residue of the Green’sfunction), which equals the area under the Lorentzian QPpeak and depends on the spectral weight transfer fromthe QP peak to plasmon satellites and the incoherentbackground.The second method is to graphically solve the non-linear QP equation in Eq. (6) using the secant (quasi-Newton) method. In this work, we refer to (cid:15) G W m ob-tained from the graphical solution as (cid:15) G W , graph m . Notethat the above linearization method corresponds to thefirst step of the secant method. Note also that when thenon-linear QP equation in Eq. (6) has multiple solutions,the secant method gives only one of them, which dependsstrongly on the choice of η .The last method is to use the position of the QP peakwith the highest spectral weight in the spectral function A mm in Eq. (7). In this work, we refer to (cid:15) G W m obtainedfrom the spectral function as (cid:15) G W , spect m and define Z spect m by replacing Z linear m and (cid:15) G W , linear m in Eq. (28) by Z spect m and (cid:15) G W , spect m , respectively. Note that we searched forthe QP peak at 0 < Z spect m < Z m in Eq. (30) because Z m represents thespectral weight, as explained above. G. G n W and G n W n Quasiparticle Energy
As introduced in Section I, there are various levels ofself-consistency in the GW approximation (from the low-est to the highest): G W , G n W , G n W n , QS GW , andSC GW . In this work, we used G n W and G n W n forsimplicity, efficiency, and stability. G n W updates onlygKS eigenvalues in G σ [ (cid:15) σi and (cid:15) σa in Eq. (24)], whereas G n W n updates gKS eigenvalues in G σ and A [ (cid:15) σi and (cid:15) σa in Eq. (19)] as well as Casida eigenvalues in W [Ω s inEq. (24)]. Therefore, G n W n is computationally more ex-pensive than G n W by the time to build and completelydiagonalize the RPA Casida matrix in Eq. (18). Notethat G n W n can be viewed as a diagonal approximationto QS GW .In this work, we obtained G n W and G n W n QP en-ergies ( (cid:15) G n W m and (cid:15) G n W n m , respectively) by iterating therecurrence relations ( n ≥ (cid:15) evGW , m = (cid:15) m + Z evGW (cid:104) ϕ m | Σ( (cid:15) m ) − v xc | ϕ m (cid:105) , (31) (cid:15) evGW,2 m = (cid:15) evGW,1 m + Z evGW (cid:104) ϕ m | Σ( (cid:15) evGW,1 m ) − Σ( (cid:15) m ) | ϕ m (cid:105) , (32) (cid:15) evGW ,nm = (cid:15) evGW ,n − m + Z evGW (cid:104) ϕ m | Σ( (cid:15) evGW ,n − m ) − Σ( (cid:15) evGW ,n − m ) | ϕ m (cid:105) , (33)where (cid:15) evGW ,nm is (cid:15) G n W m or (cid:15) G n W n m , and Z evGW = 1. Sum-ming up Eqs. (31), (32), and (33), we get (cid:15) evGW ,nm = (cid:15) m + (cid:104) ϕ m | Σ( (cid:15) evGW ,n − m ) − v xc | ϕ m (cid:105) , (34)which we refer to as the ev GW QP equation in this work.When the ev GW convergence is reached ( (cid:15) evGW ,nm = (cid:15) evGW ,n − m = (cid:15) evGW , ∞ m , where (cid:15) evGW , ∞ m are convergedev GW QP energies), the ev GW QP equation in Eq. (34)becomes similar to the G W QP equation in Eq. (6).Whereas most GW codes use 0 < Z evGW < MOLGW uses Z evGW = 1. Even though we implementedev GW with 0 < Z < GW with Z = 1 in this work for a few reasons. First, Eq. (33)shows that converged ev GW QP energies ( (cid:15) evGW , ∞ m ) isindependent of whether 0 < Z < Z = 1. Sec-ond, Z = 1 gives a unique solution that satisfies theQP equation in Eq. (34), which allows us to avoid the GW multi-solution issue from the graphical-solution andspectral-function methods and the ∼ GW with 0 < Z < GW variant that updates only a few states nearHOMO and LUMO and rigidly shifts all the other statesfor efficiency, but we updated all eigenvalues in thiswork for accuracy. Last, ev GW with Z = 1 has no vari-ant and does not need a QP equation solver, but ev GW with 0 < Z < GW variants with 0 < Z < TABLE I. Number of occupied and empty states for the ↑ -spin channel ( N ↑ occ and N ↑ emp , respectively) used in GW cal-culations. FC means the frozen-core approximation. AE andECP mean all electron and effective core potential, respec-tively. CN means the cardinal number. N ↑ occ N ↑ emp Anion Potential FC CN=2 CN=3 CN=4 CN=5ScO − AE Yes 9 67 124 205 314ScO − AE No 15 67 124 205 314TiO − AE Yes 10 66 123 204 313TiO − AE No 16 66 123 204 313CuO − AE Yes 13 63 120 201 310CuO − AE No 19 63 120 201 310CuO − ECP Yes 13 63 120 201 310CuO − ECP No 14 63 120 201 310ZnO − AE Yes 14 62 119 200 309ZnO − AE No 20 62 119 200 309ZnO − ECP Yes 14 62 119 200 309ZnO − ECP No 15 62 119 200 309
QP equation solvers give different solutions (especially,for states far away from HOMO and LUMO). Z evGW = 1 and the efficiency comparison of ev GW and G W are discusssed in supplementary material. III. COMPUTATIONAL DETAILS AND TESTRESULTSA. Computational Details
Our gKS calculations were carried out using bothMOLGW and NWChem in order to cross-check the re-sults and to ascertain the correct SCF convergence. For GW calculations, we used only MOLGW. MOs were ex-panded using augmented Dunning correlation-consistentGaussian basis sets, aug-cc-pV n Z ( n = D, T, Q, and5), which are designed to smoothly converge with basissize. Augmentation using diffuse functions is essential inground-state calculations for anions and in excited-statecalculations for both neutrals and anions. Without aug-mentation, gKS and GW eigenvalues for empty statesconverge very slowly with basis size. In the following,the cardinal number (CN = 2, 3, 4, and 5) is used to rep-resent the approximate size of diverse basis sets employedin the literature and this work. For example, CN=4means def2-QZVP in Ref. 13, and aug-cc-pVQZ in thiswork. Table I summarizes the exact size of CN=2,3,4,5basis sets used in this work. To determine the optimizedbond lengths of TMO anions, we used NWChem withPBE and CN=3. We obtained bond lengths of 1.695,1.642, 1.697, and 1.765 ˚Afor ScO − , TiO − , CuO − , andZnO − , respectively.In order to study the starting-point dependency of the GW approximation, we used global hybrid functionals: E PBE α,σ xc = αE HF ,σ x + (1 − α ) E PBE ,σ x + E PBE ,σ c , (35)where E HF ,σ x , E PBE ,σ x , and E PBE ,σ c are Fock exact ex-change, PBE exchange, and PBE correlation energies, re-spectively. We refer to the hybrid functionals in Eq. (35)as PBE α functionals in this work. While we tested otherfunctionals such as B3LYP, HSE06, BHLYP, and HF, wediscuss only PBE α (0 . ≤ α ≤ .
00) results because theEXX amount in the starting point has a stronger effecton GW results than other factors such as range sepa-ration (to screen the Coulomb interaction) and correla-tion type. As shown in supplementary material, HSE06,PBE0, and BHLYP α ( α =0.25) [where PBE is replaced byLYP in Eq. (35)] give similar GW results. Note that thetype of the correlation functional is not important (e.g.PBE vs LYP), but the existence of it is. As shown in sup-plementary material, PBE α ( α =1.00) and HF can makea large difference ( ∼ GW results for some states. B. gKS Test Results
1. Effective Core Potentials
Unlike Sc and Ti, Cu and Zn have two choices of basissets: AE (All Electron) and ECP (Effective Core Poten-tial). ECP allows to remove core electrons and includerelativistic effects. We first tested scalar relativistic ef-fects by comparing AE and ECP GW binding energies.We did not include spin-orbit coupling because (i) spin-orbit ECP is not implemented in MOLGW, and (ii) spin-orbit effects are very small in Cu and Zn, which are rel-atively light elements. The test results are presentedin supplementary material. We found that (i) ECP andAE GW IEs differ by 0.01–0.15 eV, depending on thesubtle competition between direct and indirect relativis-tic effects ( s and p orbital contraction and stabilizationand d and f orbital expansion and destabilization, re-spectively), which is consistent with Ref. 48, and (ii)ECP GW d -electron binding energies are smaller thanAE ones by 0.16–0.66 eV due to indirect effects, whichis consistent with Ref. 49. We next tested the efficiencyof ECP with respect to AE. We found that ECP is moreefficient than AE because the absence of core states notonly makes the basis size smaller, which benefits bothgKS and GW parts, but also makes the SCF cycle fasterand more stable, which benefits only the gKS part. Inthis work, we present mainly AE results not because AEis superior to ECP, but because scalar relativistic effectsmake smaller changes than large ( ∼ d -electron binding energies, where scalar rel-ativistic effects are considerable ( ∼
10% of the experimen-tal 3 d -electron binding energy). Note that Ref. 13 usedonly AE even though the GW
100 benchmark set containsAg , Cu and NCCu molecules.
2. RI for gKS
We did not use RI in this work because our goal is toassess the range of applicability of the GW approxima-tion as accurately as possible using small molecules, butRI is unavoidable for the practical GW study of largemolecules. Thus, we evaluated the quality of RI for bothAE and ECP by comparing RI and no-RI gKS eigenval-ues and total spins. The evaluation results are presentedin supplementary material. We found that CN=5 RIECP causes a large random error in gKS results (e.g. ∼ ∼ − and ZnO − , respectively, ingKS-PBE IEs), which decreases with the EXX amount.It is important to note that unlike the SCF convergenceerror, which occurs only in open-shell systems with non-zero EXX amounts, this gKS RI error occurs in bothclosed- and open-shell systems with all EXX amounts. Itis difficult to detect the gKS RI error because all SCF cy-cles with different convergence parameters smoothly con-verge to the same local minimum with no or small spincontamination. Therefore, we conclude that RI should beused only after the potential gKS RI error is thoroughlytested.Note that because we did not use RI and the 4-centerintegrals are computed at each SCF step, a single gKScalculation is as expensive as a single GW calculation inthis work, which is consistent with Ref. 9. Note also thatwe discussed the effect of RI on GW results in supple-mentary material.
3. SCF Convergence Tests
It is not straightforward to obtain the correct mean-field input for GW calculations, because successful SCFconvergence could come from both correct convergenceto a global minimum as well as wrong convergence tosome local minima. This is a particularly critical issuein gKS calculations on open-shell systems involving non-zero EXX and large basis. Many minima with similartotal energies and total spins due to nearly degenerateenergy levels in 3 d transtion metals make it more diffi-cult to obtain correct SCF convergence. For closed-shell systems or open-shell systems with (semi-)local xcfunctionals, the SCF cycle is generally guaranteed toconverge to a global minimum. However, when EXXis used for open-shell systems, wrong SCF convergenceoccurs frequently and randomly, which makes manual,time-consuming, and error-prone SCF convergence testsmandatory.In order to obtain the correct mean-field input, weperformed three-step SCF convergence tests. First, weused 12 and 96 sets of SCF convergence parameters forMOLGW and NWChem, respectively. Second, we man-ually searched for correctly converged SCF results usingmultiple indicators: gKS total energy, total spin (cid:104) S (cid:105) inEq. (26), the number of total SCF cycles, a trend overbasis size (CN=2,3,4,5), and a trend over EXX amount (by manually choosing gKS total energies and total spinsthat vary smoothly with basis size and EXX amount).Last, we cross-checked all MOLGW and NWChem gKSresults. Our SCF convergence test results are presentedin supplementary material.Note that because of our heavy SCF convergence tests, total gKS calculations are more expensive than total GW calculations in this work. C. GW Test Results
1. Complete Basis Set Limit
Like MP2, RPA, and CCSD(T) correlation energies, GW QP energies converge slowly with basis size. Accord-ingly, one should extrapolate GW QP energies obtainedfrom different basis sizes to the complete basis set (CBS)limit to avoid the incomplete basis set error of ∼ We, therefore, tested the effect of fitting function type,EXX amount, and basis size on the GW CBS limit.Two fitting functions are most widely used for the CBSlimit, which we refer to as standard fitting functions inthis work: E m = a + bN BF , (36) E m = a + b CN , (37)where E m are correlation or m th QP energies, a and b arefitting parameters, N BF is the number of basis functions(see Table I), and CN is the cardinal number. In Eqs. (36)and (37), a gives the correlation or QP energy in the CBSlimit. Note that there are various non-standard fittingfunctions used in the literature. Fig. 1 compares CBS results obtained from two stan-dard fitting functions in Eqs. (36) and (37) (as well asone non-standard one used in Refs. 6 and 54). We seethat different fitting functions always give different GW CBS limits, deviating from each other by up to ∼ GW CBS limit. We observe that the incomplete basis set er-ror increases with the EXX amount. CN=2 occasionallyand randomly causes a significant error ( ∼ GW CBS limit, which is commonly observed in the lit-erature.
Based on these test results, we conclude thatit is important to check whether extrapolation is usedor not, whether CN=2 is used or not for extrapolation,and which fitting function is used when analyzing andcomparing GW results. For example, Ref. 13 reportedthat IEs obtained from Gaussian- and planewave (PW)-based GW implementations with and without extrapola-tion, respectively, differ by ∼ GW IEs with extrapolationreduces the difference to ∼ GW CBS resultsusing the fitting function in Eq. (36) with CN=2,3,4,5. -1.4-1.2-1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0 0.005 0.01 0.015
CuO - PBE α ( α =1.00) E ne r g y [ e V ] -1.4-1.2-1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 1 2 3 4 5 6 CuO - PBE α ( α =1.00) E ne r g y [ e V ] CN GW HOMOFit (a + b/CN )Fit (a + b/CN )CBS limit (a + b/CN )CBS limit (a + b/CN ) FIG. 1. (Color online) Effect of the fitting function type on the GW complete basis set (CBS) limit. The calculations arepresented for the HOMO of CuO − at the G W @PBE α ( α =1.00) level of theory. NBF and CN represent the number of basisfunctions and the cardinal number, respectively. G W @PBE α ( α =1.00) HOMO energies are obtained from gKS-PBE α ( α =1.00)HOMO-1, which corresponds to gKS-PBE HOMO (see text). -0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 0.005 0.01 0.015 ScO - PBE α ( α =0.00) E ne r g y [ e V ] -1.8-1.7-1.6-1.5-1.4-1.3-1.2 0 0.005 0.01 0.015 ScO - PBE α ( α =1.00) E ne r g y [ e V ] FIG. 2. (Color online) Effect of the EXX amount on the GW complete basis set (CBS) limit. The calculations are presented forthe HOMO of ScO − at the G W @PBE (left) and G W @PBE α ( α =1.00) (right) levels of theory. NBF represents the numberof basis functions. We chose Eq. (36) not because it is superior to Eq. (37),but because it can be used by both Gaussian- and PW-based GW implementations.
2. Number of Empty States
By enabling MOLGW 1.F to support the largest avail-able basis set (CN=5), we also tested the effect of CN=5on the GW CBS limit. The test results are presentedin supplementary material. Here, we briefly mentiona couple of trends. In most cases, CN=5 has a small( ∼
10 meV) effect on the GW CBS limit, since CN=4,5 GW QP energies are very similar. However, in somecases, CN=5 has an appreciable ( ∼ GW CBS limit by reducing the effect of the large ran-dom CN=2 error on the GW CBS limit. In other words,CN=5 barely improves the accuracy of the GW CBS limit, but mostly acts as a bumper for the CN=2 er-ror. Moreover, CN=5 calculations are expensive (due tothe large number of empty states and the slow SCF con-vergence speed) and error-prone (due to the high chanceof SCF convergence and gKS RI errors). Therefore, weconclude that it is more beneficial to obtain the GW CBS limit from CN=3,4 than from CN=2,3,4,5. UsingCN=4 ( ∼
100 empty states per atom, as shown in Ta-ble I) instead of CN=5 as the largest basis set for the GW CBS limit tremendously reduces the computationalcosts. This conclusion is consistent with Ref. 13, whichused only CN=3,4 for extrapolation, and gives a roughestimate for two important and inter-dependent conver-gence parameters for Σ σ c ( ω ), the dimension of the dielec-trix matrix and the number of empty states, in sum-over-states PW GW implementations. The above conclusion holds only for occupied states.The effect of CN=5 on the GW CBS limit for empty0states is discussed in supplementary material. The effectof the number of occupied states on GW results usingthe frozen-core (FC) approximation, which reduces thenumber of occupied states used in the construction of G and W and thus speeds up GW calculations, is alsodiscussed in supplementary material. G W Quasiparticle Energy
A full-frequency G W method used in this work pro-duces complicated self-energy pole (and spectral-functionpeak) structures at non-frontier orbitals, so it is notstraightforward to automatically obtain correct and ac-curate (cid:15) G W m by using a single value of η and a single QPequation solver. Thus, we used three values of η (0.001,0.002, and 0.005 Ha with ∆ ω = 0.001 Ha) and three QPequation solvers mentioned in Section II F (linearization,graphical-solution, and spectral-function methods).Our analysis of a total of nine solutions for (cid:15) G W m shows that the graphical-solution method using η =0.001 Ha and the linearization method randomly give in-correct solutions. Therefore, we obtained (cid:15) G W m fromthe graphical-solution or spectral-function method using η = 0.002 or 0.005 Ha. When (cid:15) G W , graph m and (cid:15) G W , spect m are different, we manually selected a correct solution byanalyzing Σ c ( ω ) and A ( ω ).A large distance between (cid:15) m and (cid:15) G W m and multipleself-energy poles between them makes it difficult to ob-tain correct and accurate (cid:15) G W m for non-frontier orbitals,which is especially the case for (semi-)local xc function-als. Fig. 3 shows various examples of successes andfailures of three QP equation solvers for G W @PBE.In the following, we analyze each example individually.Note that we used η = 0.0001, 0.0002, and 0.0005 Hawith ∆ ω = 0.0001 Ha in a few of the following exam-ples to demonstrate the danger of small η values for thespectral-function method. Such small η values are notrecommended, as they significantly increase disk storagerequirements while barely improving accuracy.First, the top two left panels of Fig. 3 show a gen-eral example, in which all three methods succeed. Wesee a few general trends. Typically, all three methodsgive correct solutions at m = HOMO and LUMO, whichhave a simple pole structure in Σ c ( ω ). Graphical-solutionand spectral-function methods generally give multiple so-lutions, whereas the linearization method always givesa unique solution, at which two straight lines inter-sect. Generally, graphical-solution and spectral-functionmethods give identical (correct and accurate) solutions(in this case, at ω = − .
48 eV), whereas the lineariza-tion method gives a different (correct but inaccurate) so-lution (in this case, at ω = − .
69 eV) due to an intrinsicerror of ∼ η (0.0002 Ha) sharpensa weak self-energy pole at ω = − . c ( ω ), but thesharpened pole does not cause an error in the graphical-solution method because it is not between (cid:15) G W m and (cid:15) m . The very small η heightens the weak peak C in A ( ω ), but the heightened peak cannot cause an errorin the spectral-function method because it is still lowerthan other peaks A and B. In other words, the spectral-function method depends more weakly on the choice of η than the graphical-solution method. The very small η has little effect on the linearization method, because (i)the linearization method in this work depends only onΣ c ( (cid:15) m ± ∆ ω ), and (ii) ω = − . (cid:15) m to affect the finite difference method.Second, the top two right panels of Fig. 3 show a spe-cial example, in which the graphical-solution method cangive an incorrect solution. We see a few special trends.Generally, some of the three methods give incorrect solu-tions at m = HOMO- n and LUMO+ n ( n = 1, 2, 3, ...),which have a complicated pole structure in Σ c ( ω ). A verysmall η (0.0002 Ha) sharpens a weak pole at ω = − . c ( ω ), and this sharpened pole causes a large error of0.4 eV in the graphical-solution method because it is be-tween (cid:15) G W m and (cid:15) m and sharpened enough to make thesecant method fail by causing it to find an incorrect inter-section point. The very small η heightens a weak peak Bin A ( ω ), but the heightened peak cannot cause an errorin the spectral-function method because it is still lowerthan the other peak A. The very small η has little ef-fect on the linearization method because ω = − . (cid:15) m .Third, the bottom two left panels of Fig. 3 show aspecial example, in which the linearization method cangive an incorrect solution. We see a few special trends.A very small η (0.0002 Ha) sharpens a weak pole at ω ≈ (cid:15) m in Σ c ( ω ), but the sharpened pole does not cause anerror in the graphical-solution method even though it isbetween (cid:15) G W m and (cid:15) m , because the pole is not sharpenedenough ( η = 0.0001 Ha, on the other hand, causes a largeerror of 1.3 eV). The very small η heightens a weak peakat ω ≈ (cid:15) m in A ( ω ), but the heightened peak cannot causean error in the spectral-function method because it isstill lower than the other peak at ω ≈ (cid:15) G W m . The verysmall η sharpens a weak pole at ω ≈ (cid:15) m in Σ c ( ω ), andthe sharpened pole causes a large error of 2.1 eV in thelinearization method because it is too close to (cid:15) m , makingthe finite difference method fail by causing a large errorin the slope of the tangent line at ω = (cid:15) m .Last, the bottom two right panels of Fig. 3 show a spe-cial example, in which the spectral-function method cangive an incorrect solution. We see a few special trends.A very small η (0.0002 Ha) sharpens a weak pole at ω = − . c ( ω ), but the sharpened pole does not cause an error in the graphical-solution method becauseit is not between (cid:15) G W m and (cid:15) m . Two peaks A and B (at ω = − . − . A ( ω ) have similarspectral weights (and peak heights), so it is not straight-forward to unambiguously determine which one is a QPpeak or a satellite. Spectral weights (practically, peakheights) of the two peaks depend on the basis size: thepeak B (A) is higher than the peak A (B) for CN=2,3,4(CN=5). We chose peak B as the QP peak, because (i) itis consistent with the solution from the graphical-solution1 -10-5 0 5 10 -3 -2 -1 0 1 2 3 ScO - HOMOE KS E QPlin E QPgra Σ c [ e V ] ω [eV] Σ c ( η = 0.0002 Ha) Σ c ( η = 0.0020 Ha) ω - E KS + v xc - Σ x -10 0 10 20 -6 -5 -4 -3 -2 -1 0 ScO - HOMO-1E KS E QPlin E QPgra E QPgra Σ c [ e V ] ω [eV] Σ c ( η = 0.0002 Ha) Σ c ( η = 0.0020 Ha) ω - E KS + v xc - Σ x ScO - HOMOE KS E QPspe
AB C A [ e V - ] ω [eV] η = 0.0002 Ha η = 0.0020 Ha ScO - HOMO-1E KS E QPspe E QPspe
A B A [ e V - ] ω [eV] η = 0.0002 Ha η = 0.0020 Ha -10 0 10 20 -4 -3 -2 -1 0 1 2 CuO - HOMO-1E KS E QPlin E QPgra Σ c [ e V ] ω [eV] Σ c ( η = 0.0002 Ha) Σ c ( η = 0.0020 Ha) ω - E KS + v xc - Σ x -10 0 10 20 30 -8 -6 -4 -2 0 2 4 CuO - HOMO-2E KS E QPlin E QPgra Σ c [ e V ] ω [eV] Σ c ( η = 0.0020 Ha) Σ c ( η = 0.0050 Ha) ω - E KS + v xc - Σ x CuO - HOMO-1E KS E QPspe A [ e V - ] ω [eV] η = 0.0002 Ha η = 0.0020 Ha CuO - HOMO-2E KS E QPspe E QPspe
A B A [ e V - ] ω [eV]CN=4, η = 0.0050 HaCN=5, η = 0.0050 Ha FIG. 3. (Color online) Comparison of three QP equation solvers: linearization, graphical-solution, and spectral-functionmethods. (Top left) All the three methods give correct solutions. (Top right) The graphical-solution method can give incorrectsolutions. (Bottom left) The linearization method can give incorrect solutions. (Bottom right) The spectral-function methodcan give incorrect solutions. E KS represents the gKS-PBE eigenvalue. E QPlin , E QPgra , and E QPspe represent G W @PBE QP energiesobtained from linearization, graphical-solution, and spectral-function methods, respectively. The green straight line is a tangentto the red curve at ω = E KS . Except for A ( ω ) at the bottom right, all results are obtained from CN=2 no-RI AE. G W @PBE α (0 . ≤ α ≤ .
00) QP HOMO-2 energies of CuO − using the solution from the peak Bvary smoothly with α , as shown in Fig. 6]. The choiceof η and CN has little effect on the linearization method,but (cid:15) G W , linear m = − . G W @PBE binding energy.There are several points to note about the above ex-amples: (i) we chose simple examples, in which only onemethod can give an incorrect solution, for demonstrationpurposes; multiple methods can give incorrect solutionssimultaneously, as shown in supplementary material, (ii)not only a very small η , but also a very large η (e.g. ∼ n , where n = 5, 6, ...) have muchmore complicated pole [peak] structures in Σ c ( ω ) [ A ( ω )]than those in Fig. 3, so it is very difficult to choose cor-rect and accurate (cid:15) G W m for deep states not only auto-matically, but also manually.We conclude this section by summarizing several guide-lines to obtain a reliable and reproducible G W @PBEQP spectrum. First, one should try multiple η (and ∆ ω )values. There is no single general η value that workswell for all QP equation solvers, molecular systems, andmolecular orbitals. In other words, while η is typicallyviewed as a convergence parameter (the smaller η , themore accurate GW QP energy), it is practically an ad-justable parameter, which should be not too small or toolarge (e.g. 0 and ∼ η depends on | (cid:15) G W m - (cid:15) m | ,which generally decreases with the amount of EXX andincreases with the depth of the m th state. For exam-ple, when calculating G W @PBE HOMO and LOMO(lowest occupied molecular orbital) energies, one may try ∼ ∼ η .Second, we recommend using multiple QP equationsolvers. As shown in Fig. 3, the G W @PBE QP spec-trum automatically obtained from a single QP equationsolver can contain a large ( ∼ G W QP ener-gies (by ∼ GW results withoutsmall ( ∼ GW results with-out large ( ∼ GW multi-solution issues, respectively.Fourth, one should be fully aware of the large randomerrors that the linearization method, which is the mostwidely used QP equation solver, can cause. Ref. 15 sug-gests the linearization method as a preferable method fora fair comparison of G W @PBE IE (and EA) from dif-ferent GW implementations, because it gives a uniquesolution and thus is free of the GW multi-solution issue.The idea works well for the IE, but it does not perform aswell for the QP spectrum. For HOMO (and LUMO), thelinearization method generally succeeds and systemically overestimates the IE only by ∼ ∼ G W @PBE with respect to experiment. How-ever, for deep states, it randomly succeeds or fails, asshown in the bottom left of Fig. 3, and randomly overes-timates or underestimates G W @PBE binding energiesby ∼ G W @PBE QPspectrum.Last, one should be aware that different ways to handlethe GW multi-solution issue are found in the literature.For example, Ref. 13, which suggests that all solutionsare physically relevant, manually searched for an actu-ally relevant one by varying η , whereas Ref. 40 avoidedthe issue by automatically selecting the solution with thelargest Z m in Eq. (30) and using only η = 0. In this work,we adopted a combined approach. In other words, weautomatically chose the solution with the highest spec-tral weight, which is identical to the solution with thelargest Z m , as explained in Section II F, when one solu-tion is clearly more relevant than others, but we man-ually selected one solution that gives smoothly varying GW binding energies with a change in G W startingpoint and ev GW self-consistency level without causingunphysical kinks when multiple solutions are equally rel-evant [e.g. two solutions at peaks A and B in the bot-tom two right panels of Fig. 3 give similar Z m ( ∼ m = HOMO-2 of CuO − due to similar slopes of the self-energy (approximately, − . η , we adopted the approach of Ref. 13 instead of thatof Ref. 40 because η = 0 frequently causes the secantmethod in the graphical-solution method to find an in-correct intersection point and makes (cid:104) ϕ σm | Σ σ c ( ω ) | ϕ σm (cid:105) inEq. (24) diverge at ω = (cid:15) σi − Ω s and ω = (cid:15) σa +Ω s . The cu-mulant expansion, which describes plasmon satellitesbetter than the GW approximation, may allow us to ad-dress the GW multi-solution issue when the QP picturebreaks down, but it is beyond the scope of this work. G n W and G n W n Quasiparticle Energy
In this work, we used only η = 0 .
001 Ha for ev GW because it is small enough to obtain the convergence ofev GW QP energies with respect to η within ∼ GW QP energies withrespect to the iteration number are shown in supplemen-tary material. QS GW and our ev GW are quasiparticle-only GW schemes with no spectral weight transfer ( Z =1), and G n W n is a diagonal approximation to QS GW .Therefore, we compared the convergence behaviors ofQS GW and our ev GW and found a couple of similari-ties and differences between them.3 TABLE II. TM 3 d character in molecular orbitals ofTMO anions, obtained from NWChem gKS-PBE and gKS-PBE α ( α =1.00) results with CN=2 no-RI AE using the Mul-liken population analysis. The gKS-PBE orbital order is usedfor gKS-PBE α ( α =1.00) molecular orbitals (see text). Boldnumbers are used to highlight entire TM 3 d character.ScO − TiO − CuO − ZnO − ↑ ↓ ↑ ↓ PBEHOMO 0.06 α ( α =1.00)HOMO 0.00 First, the ev GW convergence is reached after onlya few iterations, which is consistent with the litera-ture. Due to the fast and stable convergence, amixing scheme is not used in our ev GW . Unlike ev GW ,QS GW generally needs ∼ Second, the orbital character affects the starting-pointdependency of ev GW . For example, we observed thatev GW QP energies for HOMO of CuO − depend morestrongly on the EXX amount in the ev GW starting pointthan those for HOMO of ScO − . We attribute this to dif-ferent amounts of 3 d character in HOMOs of ScO − andCuO − (6% and 23%, respectively, as shown in Table II).In other words, as the 3 d character in MO increases, thestarting-point dependency of ev GW increases. We alsoobserved that G n W n has a weaker (stronger) starting-point dependency for HOMO of ScO − (CuO − ) than G n W . Our observations for ScO − are consistent withRef. 41, which studied the ev GW starting-point depen-dency using small water clusters and concluded that asthe ev GW self-consistency level increases from G n W and G n W n , the ev GW starting-point dependency de-creases. However, our observations for CuO − are notconsistent with this conclusion. This is likely becauseCuO − has strong 3 d character in HOMO, whereas ScO − and small water clusters do not. This orbital-character-dependent starting-point dependency of ev GW may berelated to conflicting results for QS GW in the litera-ture: Ref. 59 showed the starting-point independency ofQS GW using a small sp − bonded molecule (CH ), whileRef. 60 showed the strong starting-point dependency ofQS GW using a d solid ( α -Fe O ). ZnO - ScO - TiO - CuO - ☎ (cid:0) ☎ (cid:0) HOMOHOMO-1
HOMO-2HOMO-3
HOMO-4
FIG. 4. (Color online) Contour plots of molecular orbitalsof TMO anions, obtained from MOLGW gKS-PBE resultswith CN=2 no-RI AE using VESTA . Red boxes are usedto highlight entire TM 3 d character. IV. RESULTS AND DISCUSSION
In this section, we compare our GW calculations toanion PES experiments, focusing especially on thefirst IE, the lowest 3 d -electron binding energy, and the or-bital order. We present our results from two approachesseperately: First, we discuss non-self-consistent GW withdifferent starting-points (namely, G W @PBE α calcula-tions as α is varied in steps of 0.25 from 0 to 1), and then,we discuss eigenvalue self-consistent GW ( G n W and G n W n ) with PBE starting point. We only briefly discussour GW results for the starting-point–self-consistencyhybrid approach, because (i) fundamentally, we foundthat the hybrid approach does not give any better re-sults than the two separate approaches, and (ii) prac-tically, the hybrid approach inherits disadvantages fromboth approaches.ScO − , TiO − , CuO − , and ZnO − are similar but differ-ent systems in several aspects. First, ScO − and CuO − are closed-shell systems, whereas TiO − and ZnO − areopen-shell systems. Second, ScO − and TiO − have par-tially filled 3 d shells, while CuO − and ZnO − have com-pletely filled 3 d shells. Third, TiO − has a shallow 3 d state, but CuO − and ZnO − have deep 3 d states. Fourth,3 d -electron photodetachment transitions are observed inTiO − and CuO − , but not in ScO − and ZnO − . Fifth,CuO − has strong 3 d character in HOMO, but ScO − ,TiO − , and ZnO − have weak 3 d character in HOMO.Last, two-electron transitions are observed in ScO − , butnot in TiO − , CuO − , and ZnO − . Due to these similaritiesand differences, TMO anions are an ideal set of systemsfor assessment of the performance of GW schemes.Table II shows the amount of TM 3 d character inall molecular orbitals considered in this work, obtainedfrom CN=2 gKS-PBE and gKS-PBE α ( α =1.00) using theMulliken population analysis. We see that gKS-PBE ↑ -HOMO of TiO − and gKS-PBE HOMO-2 of CuO − have4 TABLE III. Experimental (PES) and calculated (AE GW ) electron binding energies of TMO anions (in eV). MAE representsthe mean absolute error. Bold numbers are used to highlight 3 d -electron binding energies. G W @PBE α G n W @ G n W n @State Exp. α =0.00 α =0.25 α =0.50 α =0.75 α =1.00 PBE PBE OthersScO − ( Σ + )HOMO Σ + a g , 1.26 j , 1.19 k HOMO-1 ∆ 3.10 a g , 2.78 h , 3.31 i HOMO-2 Π 3.40 a g , 3.24 h , 3.44 i MAE 0.84 e e e e e e e TiO − ( ∆) ↑ -HOMO Σ + b g , 2.37 h , 2.34 i ↑ -HOMO-1 ∆ 1.73 b g , 1.72 i ↓ -HOMO ∆ 1.30 b g , 1.18 j , 1.14 n MAE 1.31 0.33 0.60 1.01 1.37 0.25 0.32CuO − ( Σ + )HOMO Π 1.78 c g , 1.52 j , 0.46 o HOMO-1 Σ + c g , 2.86 h , 1.60 o , 2.78 p , 2.47 q , 2.81 r HOMO-2 f g , 4.01 h , 4.58 p , 4.50 q HOMO-3 2.89 4.49 5.06 5.04 4.67 4.96 6.60HOMO-4 3.70 4.49 4.63 4.53 4.24 5.16 6.39MAE 1.69 0.47 0.27 0.40 0.75 0.29 1.37ZnO − ( Σ + ) ↑ -HOMO Σ + d g , 2.33 j , 2.29 l , 2.10 m , 1.06 o ↑ -HOMO-1 Π 2.71 d g , 1.43 o ↑ -HOMO-2 3.24 4.04 4.89 5.64 6.90 4.77 5.66 3.50 o ↓ -HOMO Π 2.40 d g , 1.20 o ↓ -HOMO-1 Σ + d g , 2.89 o MAE 1.25 0.33 0.29 0.64 1.19 0.12 0.78 a Ref. 24 b Ref. 25 c Ref. 26 d Ref. 27 e HOMO-1 and HOMO-2 are not included because our GW calculations cannot account for two-electron transitions (see text). f We chose this value from the Z band in the PES spectrum of CuO − (see text). g Ref. 62 using 6-3111+G* basis sets h Ref. 63 i Ref. 64 using the multi-reference configuration interaction (MRCI) method j Ref. 10 using the B3LYP functional k Ref. 11 using the B3LYP functional l Ref. 65 using the B3LYP functional m Ref. 21 using the G W @PBE0 method n Ref. 66 using the B3LYP functional o Ref. 56 using the G W @PBE method p Ref. 67 using the CCSD(T) method q Ref. 68 r Ref. 69 using the single and double excitation configuration interaction (SDCI) method entirely TM 3 d character. Fig. 4 shows the contour plotsof all molecular orbitals considered in this work, obtainedfrom CN=2 gKS-PBE. It is clearly seen that ↑ -HOMOof TiO − and HOMO-2 of CuO − are strongly localizedon Ti and Cu, respectivley. Table III summarizes our GW calculations with comparison to PES experimentsand existing calculations in the literature.Throughout this work, we use only gKS-PBE TM3 d character except when we discuss the subtle com- petition between direct and indirect relativistic effects,because the EXX amount has a small effect on TM3 d character. Also, throughout this work, we useonly the gKS-PBE orbital order to avoid confusion.The orbital order depends strongly on the amount ofEXX (e.g. PBE vs HF) and the level of theory (e.g.DFT vs GW ). For example, gKS-PBE HOMO ofCuO − corresponds to gKS-PBE α ( α =1.00) HOMO-1 and G W @PBE α ( α =1.00) HOMO of CuO − , as shown in5supplementary material and Fig. 6. In this work, gKS-PBE and G W @PBE were found to have the same or-bital order. A. G W Starting Points
Figs. 5 and 6 show PES and G W @PBE α (0 . ≤ α ≤ .
00) QP spectra of ScO − , TiO − , CuO − , and ZnO − . InPES spectra, vertial dashed and solid lines represent ex-perimental sp - and d -electron binding energies, respec-tively. In GW spectra, oblique dashed and solid linestrack calculated sp - and d -electron binding energies, re-spectively. In Figs. 5 and 6, we find a few general trendscommon in all TMO anions considered in this work: (i)no G W @PBE α (0 . ≤ α ≤ .
00) results are in per-fect agreement with experiment, (ii) G W @PBE under-estimates the IE of TMO anions by ∼ sp molecules( ∼ and (iii) G W @PBE α (0 . ≤ α ≤ .
50) reduces it to ∼
1. ScO − Scandium is the first transition metal and has onlyone 3 d electron. DFT and CCSD(T) calculations inRefs. 10 and 11 confirmed the ground state of ScO − as Σ + (8 σ π σ ), correcting the wrongly assumed state ∆ − (8 σ π σ δ ) in Ref. 24. There is no 3 d peak orband in the PES spectrum of ScO − , and the top three va-lence molecular orbitals have weak Sc 3 d character (6%,20%, and 18%, respectively), as shown in Table II.In the left panel of Fig. 5, we see that for HOMO-1 and HOMO-2 of ScO − , G W @PBE binding energiesslightly overestimate PES ones (by 0.20 and 0.02 eVfor HOMO-1 and HOMO-2, respectively), whereas G W @PBE α (0 . ≤ α ≤ .
00) binding energies signifi-cantly overestimate PES ones by ∼ G W @PBE α ( α =0.25) and G W @PBE α ( α =1.00)binding energies are greater than PES ones by 1.67and 2.72 eV, respectively]. This seems to suggest thatfor HOMO-1 and HOMO-2 of ScO − , G W @PBE per-forms better than G W @PBE α (0 . ≤ α ≤ . − are likely due to two-electron transitionsfrom 8 σ π σ ( Σ + ScO − ) to 8 σ π σ ( B Σ + ScO)and to 8 σ π δ ( A (cid:48) ∆ ScO) states, respectively, which GW calculations for quasiparticle excitations cannot ac-count for. In other words, the seemingly excellent agree-ment between G W @PBE and PES binding energies forHOMO-1 and HOMO-2 of ScO − is accidental. There-fore, we exclude HOMO-1 and HOMO-2 of ScO − fromour evaluation of the performance of GW schemes in thefollowing. We also see that as α increases, G W @PBE α IE al-ways increases, but this happens at different rates: As α increases from 0.00 to 0.25, it increases rapidly, whereasas α increases from 0.25 to 1.00, it increases slowly. Theweak sensitivity of G W @PBE α (0 . ≤ α ≤ .
00) IE toa change in α gives a large margin for an optimal amountof EXX: 25%–100%.
2. TiO − Titanium is the second transition metal and has two3 d electrons. Several theoretical studies in Refs. 10, 62,64, 66, and 70 confirmed 9 σ δ ( ∆) as the ground-stateelectron configuration of TiO − , correcting the wronglyassigned configuration 9 σ δ ( Σ + ) in Ref. 25. UnlikeScO − , which has an empty δ shell, TiO − has one 3 d electron in the δ shell. The transition of the 3 d electronfrom 9 σ δ ( ∆ TiO − ) to 9 σ ( Σ + TiO) states producesthe third peak in the PES spectrum of TiO − at 2.0 eV.In the G W @PBE α QP spectrum of TiO − , ↑ -HOMO isof entirely Ti 3 d character, as shown in Table II.The right panel of Fig. 5 clearly shows that G W @PBE α (0 . ≤ α ≤ .
00) binding energy for ↑ -HOMO of TiO − is much more sensitive to a change in α than those for other occupied molecular orbitals withmainly sp character, as shown in Table II. The orbital-character-dependent sensitivity of G W @PBE α bind-ing energy to a change in α causes a couple of prob-lems. First, G W @PBE underestimates the IE and the3 d -electron binding energy of TiO − non-uniformly (by0.99 and 1.74 eV, respectively), leading to the wrong or-bital order. In other words, G W does not correct thewrong orbital order produced by PBE. Second, the G W starting-point approach does not give accurate results forboth the IE and the 3 d -electron binding energy of TiO − at the same time. For example, G W @PBE α ( α =0.50)gives a better result for the IE of TiO − by 0.29 eV, but aworse result for the 3 d -electron binding energy of TiO − by 1.46 eV, than G W @PBE α ( α =0.25). This type ofbehavior is not uncommon in GW predictions for transi-tion metal oxides; for example, no existing GW schemecan accurately reproduce both the bandgap and the d -band position in the band structure of bulk ZnO at thesame time. The increase in α from 0 to 1 has a similar effect on G W @PBE α IE of both ScO − and TiO − : For bothScO − and TiO − , G W @PBE α (0 . ≤ α ≤ .
00) re-duces the underestimation of IE by G W @PBE from ∼ ∼ G W @PBE α ( α =0.25) reducesthe difference in IE between PES and G W @PBE from0.84 eV to 0.20 eV and from 0.99 eV to 0.30 eV, respec-tively]. However, unlike ScO − , the strong sensitivity of G W @PBE α d -electron binding energy of TiO − to achange in α gives a small margin for an optimal amountof EXX: ∼ ScO - exp G W @PBE α ( α =0.00) G W @PBE α ( α =0.25) G W @PBE α ( α =0.50) G W @PBE α ( α =0.75)-2 0 2 4 6 8G W @PBE α ( α =1.00) Binding Energy [eV] TiO - exp G W @PBE α ( α =0.00) updown G W @PBE α ( α =0.25) updown G W @PBE α ( α =0.50) updown G W @PBE α ( α =0.75) updown-2 0 2 4 6 8G W @PBE α ( α =1.00) Binding Energy [eV] updown FIG. 5. (Color online) Effect of the EXX amount in the G W starting point on the electronic structure of ScO − and TiO − .A Gaussian distribution function with a smearing width of 0.1 eV is used to broaden the spectra.
3. CuO − Copper is the 11th transition metal and has ten 3 d elec-trons. DFT calculations in Ref. 10 confirmed the groundstate of CuO − as Σ + with the electron configurationof 3 d pσ pπ . There are three bands (named as X,Y, and Z in Ref. 26) in the PES spectrum of CuO − , asshown in the top of the left panel of Fig. 6. Ref. 26 sug-gested that the photodetachment transition of 3 d elec-trons (3 dδ dπ dσ ) from Σ + CuO − d pσ pπ toZ CuO 3 d pσ pπ states produces the broad Z bandin the PES spectrum of CuO − at ∼ G W @PBE α QP spectrum of CuO − ,HOMO-2 is of entirely Cu 3 d character, as shown in Ta-ble II.In the left panel of Fig. 6, we see that G W @PBE α (0 . ≤ α ≤ .
50) binding energy for HOMO-2 of CuO − is more sensitive to a change in α than those for other occupied molecular orbitals withweaker Cu 3 d character than HOMO-2, as shown inTable II, and G W @PBE α (0 . ≤ α ≤ .
75) givesgood results for the IE and the 3 d -electron bindingenergy (corresponding to HOMO and HOMO-2, respec-tively) of CuO − at the same time. Scalar relativisticeffects in ECP reduce G W @PBE α ( α =0.50) and G W @PBE α ( α =0.75) binding energies for HOMO-2of CuO − by 0.31 and 0.24 eV, as shown in supplemen-tary material, without changing the conclusion that G W @PBE α (0 . ≤ α ≤ .
75) gives good results forthe 3 d -electron binding energy of CuO − . Like TiO − , theorbital-character-dependent sensitivity of G W @PBE α binding energy to a change in α causes G W @PBEto underestimate the IE and the 3 d -electron bindingenergy of CuO − non-uniformly (by 1.38 and 2.32 eV,respectively). G W @PBE α (0 . ≤ α ≤ .
00) bindingenergies for all valence molecular orbitals considered inthis work are weakly sensitive to a change in α . This7 CuO - exp X Y Z G W @PBE α ( α =0.00) G W @PBE α ( α =0.25) G W @PBE α ( α =0.50) G W @PBE α ( α =0.75)-2 0 2 4 6 8G W @PBE α ( α =1.00) Binding Energy [eV] ZnO - exp G W @PBE α ( α =0.00) updown G W @PBE α ( α =0.25) updown G W @PBE α ( α =0.50) updown G W @PBE α ( α =0.75) updown-2 0 2 4 6 8G W @PBE α ( α =1.00) Binding Energy [eV] updown FIG. 6. (Color online) Effect of the EXX amount in the G W starting point on the electronic structure of CuO − and ZnO − .A Gaussian distribution function with a smearing width of 0.1 eV is used to broaden the spectra. trend suggests that PBE α (0 . ≤ α ≤ .
00) orbitalswith large amounts of EXX are good for localizedstates of CuO − with strong 3 d character [i.e. forCuO − , PBE α (0 . ≤ α ≤ .
00) wavefunctions areclose to QP ones], and is consistent with the relativelygood performance of HF on molecules with weakscreening.
4. ZnO − Zinc is the 12th transition metal and has ten 3 d elec-trons. Zinc is rather distinct from other first row transi-tion metals due to its closed-shell electron configuration.In other words, zinc is more similar to alkaline earth met-als than other transition metals because Zn 3 d electronsgenerally do not participate in bonding. DFT calcula-tions in Ref. 10 confirmed the ground-state electron con-figuration of ZnO − as Σ + σ σ π δ . There are fourbands in the PES spectrum of ZnO − , as shown in the top of the right panel of Fig. 6. The photodetachment of 3 d electrons is not measured in the PES experiment due toinsufficient photon energy of 4.66 eV. Unlike CuO − , allvalence molecular orbitals in the G W @PBE α QP spec-trum of ZnO − have weak Zn 3 d character, as shown inTable II.The right panel of Fig. 6 shows that unlike CuO − , G W @PBE underestimates electron binding energiesfor all valence molecular orbitals of ZnO − uniformly(e.g. by 1.18, 1.35, 1.23, and 1.25 eV for ↑ -HOMO, ↑ -HOMO-1, ↓ -HOMO, and ↓ -HOMO-1, respectively) pos-sibly because all valence molecular orbitals have simi-lar Zn 3 d character and thus their G W @PBE α bind-ing energies have similar sensitivity to a change in α . G W @PBE α ( α =0.50) gives good results for the IE andthe orbital order of ZnO − at the same time. UnlikeCuO − , G W @PBE α (0 . ≤ α ≤ .
00) binding energiesfor all valence molecular orbitals, except for ↑ -HOMOand ↓ -HOMO-1, are strongly sensitive to a change in α .This trend can be explained in terms of the importance8of spin-splitting in open-shell molecules (as has been dis-cussed in Ref. 54 for CuO − molecule). B. ev GW Self-Consistency Levels
Fig. 7 shows PES and ev GW QP spectra of ScO − ,TiO − , CuO − , and ZnO − . In the following, we analyze G n W @PBE and G n W n @PBE results individually.In Fig. 7, we see that as the ev GW self-consistencylevel increases from G W to G n W n , GW bindingenergies always increase, but this occurs at differentrates: As the ev GW self-consistency level increasesfrom G W to G n W , GW binding energies increaserapidly (e.g. the IE increases by 0.70, 0.75, 1.79, and1.20 eV for ScO − , TiO − , CuO − , and ZnO − , respectiv-ley), while as it increases from G n W to G n W n , theyincrease slowly (e.g. the IE increases by 0.14, 0.15,and 0.46 eV for ScO − , TiO − , and ZnO − , respectiv-ley) except for CuO − (0.98 eV), which will be dis-cussed later. G W @PBE always underestimates elec-tron binding energies, whereas G n W n @PBE generallyoverestimates them. G n W @PBE binding energies arealways in between G W @PBE and G n W n @PBE onesand generally close to experiment. In other words, G W @PBE and G n W n @PBE act as lower and up-per bounds for G n W @PBE, generally producing over-and under-screenings, respectively. This trend of theev GW self-consistency approach in electronic structureof molecules is also observed in band structure of solids. We also see that the ev GW self-consistency has astrong effect on GW binding energies for molecular or-bitals with strong 3 d character (e.g. ↑ -HOMO of TiO − and HOMO-2 of CuO − ). For example, G n W @PBE re-duces the underestimation errors of G W @PBE in theIE and the 3 d -electron binding energy of TiO − with re-spect to experiment from 0.99 and 1.74 eV to 0.24 and0.17 eV, respectively. As a result, G n W @PBE correctsthe wrong G W @PBE orbital order in TiO − . Anotherexample is that G n W @PBE gives small ( ∼ − , which are uniformly underestimated by G W @PBE by ∼ d char-acter. For ZnO − , G W @PBE and G n W @PBE yieldmean absolute errors (MAEs) of 1.25 and 0.12 eV, re-spectively, as shown in Table III.CuO − exhibits particularly large differences between G n W @PBE and G n W n @PBE binding energies com-pared to other TMO anions. This trend is not associ-ated with scalar relativistic effects in ECP, which reduce G n W @PBE and G n W n @PBE binding energies by sim-ilar amounts (e.g. by 0.57 and 0.66 eV, respectively, forHOMO-2 of CuO − , as shown in supplementary mate-rial). We attribute this trend to strong 3 d character inmolecular orbitals of CuO − . For example, CuO − has alarger difference between G n W @PBE and G n W n @PBEIEs than ScO − (0.98 and 0.14 eV, respectively) possiblybecause CuO − has stronger 3 d character in HOMO than ScO − (23% and 6%, respectively, as shown in Table II). C. Comparison of G W starting-point and ev GW self-consistency approaches From our results presented so far, it appears thatboth G W starting-point and ev GW self-consistent ap-proaches can, in principle, be good GW methods for fi-nite systems: both G W @PBE α (0 . ≤ α ≤ .
50) and G n W @PBE can reduce the large and orbital-character-dependent non-uniform errors for electron binding en-ergies of TMO anions produced by G W @PBE withrespect to experiment from ∼ ∼ G W @PBE α ( α =0.25) and G n W @PBEgive satisfactory results for the bandgap and the d -electron binding energy of solids, and drew the con-clusions that (i) for accuracy, one can choose either G W @PBE α ( α =0.25) or G n W @PBE because theygive similar results, but (ii) for efficiency, one may wantto choose G W @PBE α ( α =0.25) over G n W @PBE be-cause the former is computationally cheaper than the lat-ter. However, in the case of molecular systems, we arguethat G n W @PBE has several practical advantages over G W @PBE α .First, G n W @PBE does not contain system-dependentadjustable parameters. Unlike extended systems, thereis no unique amount of EXX for the G W startingpoint, which works well for all finite systems. Forexample, we showed in Section IV A that 25% EXXis optimal for ScO − and TiO − , whereas 50% EXX isoptimal for CuO − and ZnO − . Also, it appears thatatoms and small molecules require more amount of EXXthan clusters and large molecules. Second, G n W @PBEis transferable between finite and extended systems. G n W @PBE works well for both molecules and solids(e.g. ZnO − anion and bulk ZnO, respectively). Thisgreatly extends the range of applicability of the GW method. For example, G n W @PBE may be applicableto solid-molecule hybrid systems such as molecular junc-tions and molecules adsorbed on solid surfaces. Also, G n W @PBE may be used for the study of quantum sizeeffects in clusters because it is independent of the clustersize. Third, G n W @PBE is easy to use and reliable. Un-like PBE α (0 . < α ≤ . G W , ev GW with Z = 1is immune to the GW multi-solution issue. Therefore, G n W @PBE does not need manual, time-consuming, anderror-prone tests to address the two issues, which are ex-plained in detail in Section III.Furthermore, G n W @PBE has a few desirable prop-erties. One of them comes from the PBE part. PBEcauses the smallest incomplete basis set error, as shownin Fig. 2, allowing one to use smaller basis sets for theCBS limit, which makes G n W @PBE cheaper. Two de-sirable properties come from the GW part. G n W (aswell as G n W n ) gives faster and more stable GW con-9 ScO - exp G W @PBE α ( α =0.00) G n W @PBE α ( α =0.00)-2 0 2 4 6 8G n W n @PBE α ( α =0.00) Binding Energy [eV] TiO - exp G W @PBE α ( α =0.00) updown G n W @PBE α ( α =0.00) updown-2 0 2 4 6 8G n W n @PBE α ( α =0.00) Binding Energy [eV] updown CuO - exp X Y Z G W @PBE α ( α =0.00) G n W @PBE α ( α =0.00)-2 0 2 4 6 8G n W n @PBE α ( α =0.00) Binding Energy [eV] ZnO - exp G W @PBE α ( α =0.00) updown G n W @PBE α ( α =0.00) updown-2 0 2 4 6 8G n W n @PBE α ( α =0.00) Binding Energy [eV] updown FIG. 7. (Color online) Effect of the ev GW self-consistency level on the electronic structure of ScO − , TiO − , CuO − , and ZnO − .A Gaussian distribution function with a smearing width of 0.1 eV is used to broaden the spectra. vergence than QS GW and depends more weakly on thechoice of η (e.g. we used a single value of η for ev GW inthis work), as discussed in Section III C 4. Also, G n W is cheaper than G n W n , as pointed out in Ref. 17 anddiscussed in Section II G. In fact, G n W is the cheapestself-consistent GW scheme.One may argue that G W @PBE α should be a choiceof GW methods because it is computationally more effi-cient than G n W @PBE by the number of self-consistent G n W iterations. However, as discussed in Section II G and supplementary material, this is not the case sincethe compute time difference between G W @PBE α and G n W @PBE does not depend only on the number of G n W iterations; there are other factors such as the num-ber of eigenvalues to update for (cid:15) G n W m , the number offrequency points to use for Σ c ( ω ), the number of ∆ ω and η values to test for (cid:15) G W m , and the number of initial guesswavefunctions to test for gKS calculations. Some factorscan cancel each other out; for example, G n W @PBE re-quires a few G n W iterations, but one typically needs to0test a few η values for G W @PBE α . In other words,when all factors are taken into account, the total com-pute time to obtain reliable and reproducible QP spectraat G W @PBE α and G n W @PBE levels of theory canbe comparable, as is especially the case for open-shellsystems. D. Comparison with results in the literature
Some of our results for the performance of G W starting-point and ev GW self-consistency approaches inthis work may seem to be at odds with some of the resultsin the literature. In this section, we discuss the origin ofthe apparent differences between them.We begin with the G W starting-point approach. Ta-ble IV summarizes a few selected results for the optimalamount of EXX in the G W starting point out of nu-merous results, such as Refs. 14 and 73, in the literature.Interestingly, we see that there is a wide range of EXXamounts from 25% to 100%, and Refs. 4 and 23 obtaineddifferent results (50% and 100%, respectively) from thesame set of molecules. It seems that 75% and 100% aretoo large compared to our results: 25–50%. One mayguess that the large difference is due to implementationdifferences such as basis type (e.g. Gaussian vs PW) andfrequency integration type (e.g. analytical vs numerical).However, Refs. 15 and 16 showed that such implementa-tion differences have little effect on G W IE ( ∼ G W results than implementation differences.One factor is the choice of system and property. Asshown in Section IV A, G W @PBE α (0 . ≤ α ≤ . sp systems are slightly different (by ∼ G W studies used the IE of sp -bondedsystems to determine the optimal amount of EXX in the G W starting point. The other factor is that the choiceof QP equation solver and CBS extrapolation method.As shown in Section III C 1 and Section III C 3, the lin-earization method and the CBS extrapolation method(e.g. whether to extrapolate or not and which fittingfunction and basis set to use for extrapolation) can causea difference in G W IE on the order of ∼ G W startingpoint, and thus is likely to produce the wide range ofamounts that exist in the literature.Next, we move on to the ev GW self-consistency ap-proach, and discuss the origin of apparently conflict-ing ev GW results for IE and starting-point dependency.First, Ref. 7 reported that the G n W n approach with alocal-density approximation starting point gives good re-sults for the IE of large sp molecules, whereas we foundin Section IV B that G n W @PBE gives satisfactory re-sults for the electronic structure (including the IE) ofsmall 3 d molecules. A comparison of ev GW implemen-tations in Ref. 7 and this work is provided in supplemen-tary material. We believe that the main origin of the different results is the orbital-character-dependent sensi-tivity of ev GW binding energy to a change in ev GW self-consistency level. As shown in Section IV B, G n W @PBEand G n W n @PBE binding energies are slightly differ-ent for delocalized HOMO with weak 3 d character by ∼ d character by ∼ ∼ GW IE,as shown in Section III C. Accordingly, they are mostlikely not the reason for the large (0.98 eV) differencebetween G n W @PBE and G n W n @PBE IEs of CuO − .Second, Ref. 41 reported that in small water clusters,as the ev GW self-consistency level increases, the ev GW starting-point dependency decreases, whereas we foundin Section III C 4 that in TMO anions, G n W n sometimesdepends more strongly on the starting point than G n W .As mentioned in Section III C 4, we believe that the or-bital character influences the ev GW starting-point de-pendency: for molecular orbitals with strong (weak) 3 d character, G n W n depends more strongly (weakly) on thestarting point than G n W . Overall, without molecularorbitals with strong 3 d character (e.g. HOMO of CuO − ),our ev GW results for IE and starting-point dependencyin this work are consistent with those in Refs. 7 and 41.To verify our idea about the origin of the seeminglydifferent results between this work and the literature, weperformed a simple test: (i) we chose ScO − and TiO − asour analogs of sp molecules in the literature because theirvalence molecular orbitals have weak transition-metalcharacter, except for ↑ -HOMO of TiO − with entirely Ti3 d character, (ii) we applied 15 different starting-point–self-consistency hybrid GW schemes ( G W , G n W , and G n W n ; 0%, 25%, 50%, 75%, and 100% EXX) to them,and (iii) we searched for GW schemes that give a rea-sonably small error of less than 0.5 eV in the IE and the3 d -electron binding energy with respect to experiment.Fig. 8 shows the results of the test. We see that GW IEs of ScO − and TiO − (red dashed lines) depend weaklyon the starting point and the self-consistency level, giv-ing a large margin for the choice of GW schemes. 14 GW schemes out of 15 ( G W @PBE is an exception asexpected) give a small error (less than 0.5 eV), whichexplains why there are a large number of different good GW schemes for the IE of sp molecules in the literature.We also see that the GW d -electron binding energy ofTiO − (green solid lines) depends strongly on the startingpoint and the self-consistency level, yielding a small mar-gin for the choice of GW schemes. Only two GW schemes( G W @PBE0 and G n W @PBE) out of 15 give a smallerror (less than 0.5 eV), which is why we obtained a smallnumber of good GW schemes for the electronic structureof d molecules in this work. Overall, we confirm thatevaluation results for the performance of GW schemesdepend strongly on the choice of system and property(e.g. the IE with mainly sp character vs the electronicstructure containing d states).1 TABLE IV. Optimal amount of EXX in the G W starting point for gas-phase small molecules (highlighted in bold).Reference K¨orbel et al. Bruneval et al. Kaplan et al. Rostgaard et al. This workCode FIESTA MOLGW TURBOMOLE GPAW MOLGWOptimal EXX
25% 50% 75% 100% 25–50%
Tested EXX 25 & 100% 0, 20, 25 & 50% 0, 25 & 75% 0 & 100% 0, 25, 50, 75 & 100%System 39 closed-shell 3,4,5 d
34 closed-shell sp a
29 closed-shell sp
34 closed-shell sp a d & 9 closed-shell sp System size 2–7 atoms 2–8 atoms 2–18 atoms 2–8 atoms 2 atomsProperty HOMO & LUMO HOMO HOMO & HOMO HOMO- n ( n = 0, 1, ...)HOMO- n ( n = 0, 1, ...) b (focusing on 3 d MO)Reference data Experiment ∆SCF c QS GW d Experiment Experiment ω integration Contour deformation e Fully analytic Fully analytic Fully analytic Fully analyticQP equation Linearization Linearization Spectral function Linearization f & Graphical solution &Spectral function g Spectral function η Not available Not available 0.001 eV 0 eV 0.002 or 0.005 HaCBS limit Not used Not used Not used Not used Used [employing Eq. (36)](CN=4 only) (CN=4 only) (CN=3 only) (CN=2 only) (CN=2,3,4,5)Potential ECP AE AE PAW h AERI Used Not used Used Not applicable i Not used a The same set of molecules is used. b For naphthalene only c Ref. 4 showed that ∆SCF using CCSD(T) with CN=4 causes an error of ∼ sp molecules with respect toexperiment (the largest being 0.67 eV for NaCl). d Ref. 55 showed that QS GW with CN=5 causes a mean absolute error of 0.18 eV in the IE of the first row atoms with respect toexperiment (the largest being ∼ e Refs. 13 and 34 showed that the contour deformation technique produces almost the same GW self-energy as the fully analytic methodfor frontier and non-frontier orbitals, respectively. f For 0% EXX g For 100% EXX h Projector-Augmented Wave i GPAW uses augmented Wannier basis sets, whereas FIESTA, MOLGW, and TURBOMOLE use Gaussian basis sets. -4-3-2-1 0 1 2 3 4 0 0.25 0.5 0.75 1
ScO - ≤ α opt ≤ E c a l b i nd - E e x pb i nd [ e V ] α HOMO G W @PBE α -4-3-2-1 0 1 2 3 4 0 0.25 0.5 0.75 1 ScO - ≤ α opt ≤ E c a l b i nd - E e x pb i nd [ e V ] α HOMO G n W @PBE α -4-3-2-1 0 1 2 3 4 0 0.25 0.5 0.75 1 ScO - ≤ α opt ≤ E c a l b i nd - E e x pb i nd [ e V ] α HOMO G n W n @PBE α -4-3-2-1 0 1 2 3 4 0 0.25 0.5 0.75 1 TiO - α opt = 0.25 E c a l b i nd - E e x pb i nd [ e V ] α down-HOMO G W @PBE α up-HOMO G W @PBE α up-HOMO-1 G W @PBE α -4-3-2-1 0 1 2 3 4 0 0.25 0.5 0.75 1 TiO - α opt = 0.00 E c a l b i nd - E e x pb i nd [ e V ] α down-HOMO G n W @PBE α up-HOMO G n W @PBE α up-HOMO-1 G n W @PBE α -4-3-2-1 0 1 2 3 4 0 0.25 0.5 0.75 1 TiO - E c a l b i nd - E e x pb i nd [ e V ] α down-HOMO G n W n @PBE α up-HOMO G n W n @PBE α up-HOMO-1 G n W n @PBE α FIG. 8. (Color online) Effect of the GW starting point and the ev GW self-consistency level on the electronic structure of ScO − and TiO − . E expbind and E calbind represent experimental and calculated electron binding energies, respectively. Dashed and solidlines track sp - and d -electron binding energies, respectively. α opt represents an optimal fraction of EXX in the GW startingpoint. V. SUMMARY AND CONCLUSIONS
In summary, we calculated the electronic structure ofclosed- and open-shell molecular anions with partiallyand completely filled 3 d shells (shallow and deep 3 d states, respectively) using various GW schemes and com-pared calculated GW QP spectra to anion PES experi-ments to evaluate the performance of the GW approxi-mation on both localized and delocalized states of smallmolecules containing 3 d transition metals.We found that the perturbative one-shot G W @PBEscheme, which is the most widely used GW scheme forextended systems, has a couple of problems for finite sys-tems. Fundamentally, G W @PBE underestimates theIE and the 3 d -electron binding energy by ∼ ∼ ∼ GW binding energies, G W @PBEsometimes gives the incorrect orbital order. Practically, G W @PBE suffers from the GW multi-solution issuedue to the large distance between QP and gKS-PBEeigenvalues and the complicated pole (peak) structurein the self-energy (the spectral function).We found that the G W starting-point approach, G W @PBE α , can improve G W @PBE at the expenseof introducing a couple of problems. The G W starting-point approach can give good results for the IE and the3 d -electron binding energy at the same time, and thus,correct the wrong orbital order produced by PBE. Also,the G W starting-point approach can mitigate the GW multi-solution issue by reducing the distance between QPand gKS eigenvalues. However, the optimal amount ofEXX in the G W starting point depends strongly on theamount of 3 d character in molecular orbitals, leading tothe strong sensitivity of 3 d -electron binding energy to achange in the EXX amount. Thus, the optimal amount ofEXX is strongly system- and property-dependent. Moreimportantly, G W @PBE α suffers from the SCF con-vergence issue in open-shell systems, which is absent in G W @PBE.We found that the eigenvalue self-consistency ap-proaches, G n W @PBE and G n W n @PBE, can improve G W @PBE, too. Especially, G n W @PBE gives as goodresults for the IE and the 3 d -electron binding energy as G W @PBE α without suffering from GW multi-solutionand SCF convergence issues.We recommend G n W @PBE because of its practical advantages: (i) G n W @PBE is transferable, because itgives satisfactorily accurate results for both finite andextended systems, for both closed- and open-shell sys-tems, and for both localized and delocalized states, (ii) G n W @PBE is predictive, because it does not need anysystem- and property-dependent parameters, and (iii) G n W is efficient and easy to use, because it does notrequire computational and human efforts to address SCFconvergence and GW multi-solution issuesWe attribute the good performance of G n W @PBE to the fortuitous cancellation effect: the overscreeningof the Coulomb interaction due to the over-delocalizedPBE wavefunction is cancelled by the underscreening dueto the neglect of vertex corrections. In other words,for G W applied to finite systems, PBE is a “bad”starting point in the sense that it causes a large ( ∼ G n W appliedto finite and extended systems, PBE is a “good” startingpoint in the sense that it accidentally produces the over-screening just as much as vertex corrections do, which ismissing in self-consistent GW schemes.Our results in this work – (i) G W @PBE α (0 . ≤ α ≤ .
50) and G n W @PBE give good QP energies for molecu-lar orbitals with both weak and strong 3 d character, and(ii) the ev GW starting-point dependency is more relatedto the orbital character than the self-consistency level –may seem to disagree with some results in the literature,but this is not the case. The origin of the seeming dis-agreement is that except for G W @PBE α (0 . ≤ α ≤ . ∼ sp character, which is ac-cidentally comparable to individual or combined errorsfrom multiple sources, such as the incomplete basis set,the linearization method in G W , and the insufficientnumber of eigenvalues to update in ev GW . G n W @PBE is not a conserving and starting-point-independent GW scheme. It is not the most accurateor efficient GW scheme, either. However, G n W @PBEgives satisfactory and reliable results for a wide rangeof systems, such as solids with strong screening andmolecules with weak screening, at moderate computa-tional and minimal human efforts, and thus is ideal forautomated mass GW and BSE calculations for high-throughput screening and machine learning. Furtherstudies on the performance of more diverse GW schemeson larger and more complex systems containing a broaderrange of transition metals are needed to extend the rangeof applicability of the GW approximation. VI. SUPPLEMENTARY MATERIAL
See supplementary material for more details, results,and discussion.
ACKNOWLEDGMENTS
This work was supported by the U.S. Department ofEnergy Grant No. de-sc0017824. The computation forthis work was performed on the high performance com-puting infrastructure provided by Research ComputingSupport Services at the University of Missouri-Columbia.This research also used resources of the National EnergyResearch Scientific Computing Center, a DOE Office ofScience User Facility supported by the Office of Science3of the U.S. Department of Energy under Contract No.DE-AC02-05CH11231. We also would like to thank Bin Shi and Meisam Rezaei for useful discussions in the ear-lier stage of this work. ∗ [email protected] F. Bruneval and M. Gatti, “Quasiparticle self-consistentgw method for the spectral properties of complex materi-als,” in
First Principles Approaches to Spectroscopic Prop-erties of Complex Materials , edited by C. Di Valentin,S. Botti, and M. Cococcioni (Springer Berlin Heidelberg,Berlin, Heidelberg, 2014) pp. 99–135. L. Reining, “The gw approximation: content, suc-cesses and limitations,” Wiley Interdisciplinary Reviews:Computational Molecular Science , e1344 (2018),https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1344. X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko,A. Sanfilippo, K. Reuter, and M. Scheffler, “Resolution-of-identity approach to hartreefock, hybrid density func-tionals, rpa, mp2 and gw with numeric atom-centered or-bital basis functions,” New Journal of Physics , 053020(2012). F. Bruneval and M. A. L. Marques, “Benchmark-ing the starting points of the gw approximationfor molecules,” Journal of Chemical Theory andComputation , 324–329 (2013), pMID: 26589035,https://doi.org/10.1021/ct300835h. E. Pavarini, E. Koch, J. van den Brink, and G. Sawatzky,
Quantum Materials: Experiments and Theory , Modelingand Simulation, Vol. 6 (Forschungszentrum Jlich, Jlich,2016) p. 420 p. F. Bruneval, T. Rangel, S. M. Hamed, M. Shao, C. Yang,and J. B. Neaton, “molgw 1: Many-body perturbation the-ory software for atoms, molecules, and clusters,” ComputerPhysics Communications , 149 – 161 (2016). X. Blase, C. Attaccalite, and V. Olevano, “First-principles GW calculations for fullerenes, porphyrins, phtalocyanine,and other molecules of interest for organic photovoltaicapplications,” Phys. Rev. B , 115103 (2011). M. J. van Setten, F. Weigend, and F. Evers, “Thegw-method for quantum chemistry applications: The-ory and implementation,” Journal of Chemical Theoryand Computation , 232–246 (2013), pMID: 26589026,https://doi.org/10.1021/ct300648t. J. Wilhelm, M. Del Ben, and J. Hutter, “Gw inthe gaussian and plane waves scheme with applica-tion to linear acenes,” Journal of Chemical Theory andComputation , 3623–3635 (2016), pMID: 27348184,https://doi.org/10.1021/acs.jctc.6b00380. G. L. Gutsev, B. K. Rao, and P. Jena, “Electronicstructure of the 3d metal monoxide anions,” The Jour-nal of Physical Chemistry A , 5374–5379 (2000),https://doi.org/10.1021/jp000384f. J. M. Gonzales, R. A. King, and H. F. Schaefer,“Analyses of the sco − and sco − photoelectron spectra,”The Journal of Chemical Physics , 567–572 (2000),https://doi.org/10.1063/1.481832. Kudin, Konstantin N. and Scuseria, Gustavo E., “Converg-ing self-consistent field equations in quantum chemistry -recent achievements and remaining challenges,” ESAIM:M2AN , 281–296 (2007). M. J. van Setten, F. Caruso, S. Sharifzadeh, X. Ren,M. Scheffler, F. Liu, J. Lischner, L. Lin, J. R. Deslippe,S. G. Louie, C. Yang, F. Weigend, J. B. Neaton, F. Ev-ers, and P. Rinke, “Gw100: Benchmarking g0w0 formolecular systems,” Journal of Chemical Theory andComputation , 5665–5687 (2015), pMID: 26642984,https://doi.org/10.1021/acs.jctc.5b00453. F. Caruso, M. Dauth, M. J. van Setten, andP. Rinke, “Benchmark of gw approaches for thegw100 test set,” Journal of Chemical Theory andComputation , 5076–5087 (2016), pMID: 27631585,https://doi.org/10.1021/acs.jctc.6b00774. E. Maggio, P. Liu, M. J. van Setten, and G. Kresse,“Gw100: A plane wave perspective for smallmolecules,” Journal of Chemical Theory and Com-putation , 635–648 (2017), pMID: 28094981,https://doi.org/10.1021/acs.jctc.6b01150. M. Govoni and G. Galli, “Gw100: Comparisonof methods and accuracy of results obtained withthe west code,” Journal of Chemical Theory andComputation , 1895–1909 (2018), pMID: 29397712,https://doi.org/10.1021/acs.jctc.7b00952. M. Shishkin and G. Kresse, “Self-consistent gw calcula-tions for semiconductors and insulators,” Phys. Rev. B ,235102 (2007). F. Fuchs, J. Furthm¨uller, F. Bechstedt, M. Shishkin, andG. Kresse, “Quasiparticle band structure based on a gen-eralized kohn-sham scheme,” Phys. Rev. B , 115109(2007). J. c. v. Klimeˇs, M. Kaltak, and G. Kresse, “Predictive gw calculations using plane waves and pseudopotentials,”Phys. Rev. B , 075125 (2014). M. Grumet, P. Liu, M. Kaltak, J. c. v. Klimeˇs, andG. Kresse, “Beyond the quasiparticle approximation: Fullyself-consistent gw calculations,” Phys. Rev. B , 155143(2018). S. K¨orbel, P. Boulanger, I. Duchemin, X. Blase,M. A. L. Marques, and S. Botti, “Benchmark many-body gw and bethesalpeter calculations for small transi-tion metal molecules,” Journal of Chemical Theory andComputation , 3934–3943 (2014), pMID: 26588537,https://doi.org/10.1021/ct5003658. F. Kaplan, M. E. Harding, C. Seiler, F. Weigend, F. Ev-ers, and M. J. van Setten, “Quasi-particle self-consistentgw for molecules,” Journal of Chemical Theory andComputation , 2528–2541 (2016), pMID: 27168352,https://doi.org/10.1021/acs.jctc.5b01238. C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen, “Fullyself-consistent gw calculations for molecules,” Phys. Rev.B , 085103 (2010). H. Wu and L.-S. Wang, “Photoelectron spectroscopy andelectronic structure of scon- (n = 14) and yon- (n = 15):strong electron correlation effects in sco-and yo-,” TheJournal of Physical Chemistry A , 9129–9135 (1998),https://doi.org/10.1021/jp982588q. H. Wu and L. Wang, “Electronic structure of titaniumoxide clusters: Tioy (y=13) and (tio2)n (n=14),” The Journal of Chemical Physics , 8221–8228 (1997),https://doi.org/10.1063/1.475026. H. Wu, S. R. Desai, and L.-S. Wang, “Chemical bond-ing between cu and oxygencopper oxides vs o2 complexes:a study of cuox (x = 06) species by anion photoelectronspectroscopy,” The Journal of Physical Chemistry A ,2103–2111 (1997), https://doi.org/10.1021/jp9631442. V. D. Moravec, S. A. Klopcic, B. Chatterjee, and C. C.Jarrold, “The electronic structure of zno and znf deter-mined by anion photoelectron spectroscopy,” ChemicalPhysics Letters , 313 – 318 (2001). F. Bruneval, N. Vast, L. Reining, M. Izquierdo, F. Sirotti,and N. Barrett, “Exchange and correlation effects in elec-tronic excitations of cu O,” Phys. Rev. Lett. , 267601(2006). L. Y. Isseroff and E. A. Carter, “Importance of referencehamiltonians containing exact exchange for accurate one-shot gw calculations of cu o,” Phys. Rev. B , 235142(2012). M. L. Tiago and J. R. Chelikowsky, “Optical excitationsin organic molecules, clusters, and defects studied by first-principles green’s function methods,” Phys. Rev. B ,205334 (2006). A. Seidl, A. G¨orling, P. Vogl, J. A. Majewski, and M. Levy,“Generalized kohn-sham schemes and the band-gap prob-lem,” Phys. Rev. B , 3764–3774 (1996). T. Sander, E. Maggio, and G. Kresse, “Beyond the tamm-dancoff approximation for extended systems using exactdiagonalization,” Phys. Rev. B , 045209 (2015). Y.-M. Byun and C. A. Ullrich, “Excitons in solids fromtime-dependent density-functional theory: Assessing thetamm-dancoff approximation,” Computation , 9 (2017). D. Golze, J. Wilhelm, M. J. van Setten, andP. Rinke, “Core-level binding energies from gw:An efficient full-frequency approach within a local-ized basis,” Journal of Chemical Theory and Com-putation , 4856–4869 (2018), pMID: 30092140,https://doi.org/10.1021/acs.jctc.8b00458. A. J. Cohen, D. J. Tozer, and N. C. Handy,“Evaluation of s2 in density functional theory,” TheJournal of Chemical Physics , 214104 (2007),https://doi.org/10.1063/1.2737773. A. S. Menon and L. Radom, “Consequences of spincontamination in unrestricted calculations on open-shellspecies: Effect of hartreefock and mllerplesset contri-butions in hybrid and double-hybrid density functionaltheory approaches,” The Journal of Physical Chem-istry A , 13225–13230 (2008), pMID: 18759419,https://doi.org/10.1021/jp803064k. J. G. Hill and J. A. Platts, “Auxiliary basis sets for densityfittingmp2 calculations: Nonrelativistic triple- all-electroncorrelation consistent basis sets for the 3d elements sczn,”The Journal of Chemical Physics , 044104 (2008),https://doi.org/10.1063/1.2826348. J. G. Hill and K. A. Peterson, “Explicitly correlated cou-pled cluster calculations for molecules containing group 11(cu, ag, au) and 12 (zn, cd, hg) elements: Optimized com-plementary auxiliary basis sets for valence and corevalencebasis sets,” Journal of Chemical Theory and Computation , 518–526 (2012), https://doi.org/10.1021/ct200856f. M. S. Hybertsen and S. G. Louie, “Electron correlation insemiconductors and insulators: Band gaps and quasiparti-cle energies,” Phys. Rev. B , 5390–5413 (1986). M. V´eril, P. Romaniello, J. A. Berger, and P.-F. Loos, “Unphysical discontinuities in gw meth-ods,” Journal of Chemical Theory and Compu-tation , 5220–5228 (2018), pMID: 30212627,https://doi.org/10.1021/acs.jctc.8b00745. X. Blase, P. Boulanger, F. Bruneval, M. Fernandez-Serra,and I. Duchemin, “Gw and bethe-salpeter study of smallwater clusters,” The Journal of Chemical Physics ,034109 (2016), https://doi.org/10.1063/1.4940139. T. Kotani, M. van Schilfgaarde, and S. V. Faleev,“Quasiparticle self-consistent gw method: A basis for theindependent-particle approximation,” Phys. Rev. B ,165106 (2007). C. L. Collins, K. G. Dyall, and H. F. Schaefer,“Relativistic and correlation effects in cuh, agh, andauh: Comparison of various relativistic methods,” TheJournal of Chemical Physics , 2024–2031 (1995),https://doi.org/10.1063/1.468724. K. A. Peterson and C. Puzzarini, “Systematically conver-gent basis sets for transition metals. ii. pseudopotential-based correlation consistent basis sets for the group 11 (cu,ag, au) and 12 (zn, cd, hg) elements,” Theoretical Chem-istry Accounts , 283–296 (2005). P. Pyykk¨o and J. P. Desclaux, “Relativity and the periodicsystem of elements,” Accounts of Chemical Research ,276–281 (1979), https://doi.org/10.1021/ar50140a002. P. Pyykk¨o, “Relativistic effects in structural chem-istry,” Chemical Reviews , 563–594 (1988),https://doi.org/10.1021/cr00085a006. P. Pyykk¨o, “Relativistic effects in chemistry: Morecommon than you thought,” Annual Review of Phys-ical Chemistry , 45–64 (2012), pMID: 22404585,https://doi.org/10.1146/annurev-physchem-032511-143755. P. S. Bagus, Y. S. Lee, and K. S. Pitzer, “Effects of rela-tivity and of the lanthanide contraction on the atoms fromhafnium to bismuth,” Chemical Physics Letters , 408 –411 (1975). H. Tatewaki, S. Yamamoto, and Y. Hatano,“Relativistic effects in the electronic structureof atoms,” ACS Omega , 6072–6080 (2017),https://doi.org/10.1021/acsomega.7b00802. C. H¨attig, W. Klopper, A. K¨ohn, and D. P.Tew, “Explicitly correlated electrons in molecules,”Chemical Reviews , 4–74 (2012), pMID: 22206503,https://doi.org/10.1021/cr200168z. F. H¨user, T. Olsen, and K. S. Thygesen, “Quasiparticlegw calculations for solids, molecules, and two-dimensionalmaterials,” Phys. Rev. B , 235132 (2013). T. Rangel, S. M. Hamed, F. Bruneval, and J. B. Neaton,“Evaluating the gw approximation with ccsd(t) for chargedexcitations across the oligoacenes,” Journal of ChemicalTheory and Computation , 2834–2842 (2016), pMID:27123935, https://doi.org/10.1021/acs.jctc.6b00163. L. Hung, F. Bruneval, K. Baishya, and S. ¨O˘g¨ut,“Benchmarking the gw approximation and bethe-salpeter equation for groups ib and iib atoms andmonoxides,” Journal of Chemical Theory and Com-putation , 2135–2146 (2017), pMID: 28387124,https://doi.org/10.1021/acs.jctc.7b00123. B. Shi, S. Weissman, F. Bruneval, L. Kronik, and S. ¨O˘g¨ut,“Photoelectron spectra of copper oxide cluster anions fromfirst principles methods,” The Journal of Chemical Physics , 064306 (2018), https://doi.org/10.1063/1.5038744. F. Bruneval, “Ionization energy of atoms obtained fromgw self-energy or from random phase approximation totalenergies,” The Journal of Chemical Physics , 194107(2012), https://doi.org/10.1063/1.4718428. L. Hung, F. Bruneval, K. Baishya, and S. ¨O˘g¨ut,“Correction to benchmarking the gw approximationand bethe-salpeter equation for groups ib and iibmonoxides,” Journal of Chemical Theory and Com-putation , 5820–5821 (2017), pMID: 29099593,https://doi.org/10.1021/acs.jctc.7b01054. L. Hedin, “On correlation effects in electron spectroscopiesand the gw approximation,” Journal of Physics: Con-densed Matter , R489–R528 (1999). J. S. Zhou, J. J. Kas, L. Sponza, I. Reshetnyak, M. Guzzo,C. Giorgetti, M. Gatti, F. Sottile, J. J. Rehr, andL. Reining, “Dynamical effects in electron spectroscopy,”The Journal of Chemical Physics , 184109 (2015),https://doi.org/10.1063/1.4934965. P. Koval, D. Foerster, and D. S´anchez-Portal, “Fullyself-consistent gw and quasiparticle self-consistent gw formolecules,” Phys. Rev. B , 155417 (2014). P. Liao and E. A. Carter, “Testing variations of the gwapproximation on strongly correlated transition metal ox-ides: hematite ( α -fe2o3) as a benchmark,” Phys. Chem.Chem. Phys. , 15189–15199 (2011). K. Momma and F. Izumi, “Vesta: a three-dimensional vi-sualization system for electronic and structural analysis,”Journal of Applied Crystallography , 653–658 (2008),https://onlinelibrary.wiley.com/doi/pdf/10.1107/S0021889808012016. B. Dai, K. Deng, J. Yang, and Q. Zhu, “Excitedstates of the 3d transition metal monoxides,” TheJournal of Chemical Physics , 9608–9613 (2003),https://doi.org/10.1063/1.1570811. A. J. Bridgeman and J. Rothery, “Periodic trends in thediatomic monoxides and monosulfides of the 3d transitionmetals,” J. Chem. Soc., Dalton Trans. , 211–218 (2000). E. Miliordos and A. Mavridis, “Electronic structure andbonding of the early 3d-transition metal diatomic oxides and their ions: Sco0,, tio0,, cro0,, and mno0,,” The Journalof Physical Chemistry A , 8536–8572 (2010), pMID:20113002, https://doi.org/10.1021/jp910218u. C. W. Bauschlicher and H. Partridge, “A comparison ofzno and zno,” The Journal of Chemical Physics , 8430–8434 (1998), https://doi.org/10.1063/1.477506. M. B. Walsh, R. A. King, and H. F. Schaefer, “The struc-tures, electron affinities, and energetic stabilities of tionand tion (n=13),” The Journal of Chemical Physics ,5224–5230 (1999), https://doi.org/10.1063/1.478418. H. Xian, Z. Cao, X. Xu, X. Lu, and Q. Zhang, “Theoret-ical study of low-lying electronic states of cuo and cuo,”Chemical Physics Letters , 485 – 493 (2000). A. Daoudi, A. T. Benjelloun, J. Flament, and G. Berthier,“Potential energy curves and electronic structure of coppernitrides cun and cun+versus cuo and cuo+,” Journal ofMolecular Spectroscopy , 8 – 16 (1999). P. S. Bagus, C. J. Nelin, and C. W. Bauschlicher, “On thelowlying states of cuo,” The Journal of Chemical Physics , 2975–2981 (1983), https://doi.org/10.1063/1.446126. G. L. Gutsev, L. Andrews, and C. W. Bauschlicher, Jr.,“Similarities and differences in the structure of 3d-metalmonocarbides and monoxides,” Theoretical Chemistry Ac-counts , 298–308 (2003). C. A. Fancher, H. L. de Clercq, O. C. Thomas, D. W.Robinson, and K. H. Bowen, “Zinc oxide and its an-ion: A negative ion photoelectron spectroscopic study,”The Journal of Chemical Physics , 8426–8429 (1998),https://doi.org/10.1063/1.477505. M. Strange, C. Rostgaard, H. H¨akkinen, and K. S. Thyge-sen, “Self-consistent gw calculations of electronic transportin thiol- and amine-linked molecular junctions,” Phys. Rev.B , 115108 (2011). N. Marom, F. Caruso, X. Ren, O. T. Hofmann,T. K¨orzd¨orfer, J. R. Chelikowsky, A. Rubio, M. Schef-fler, and P. Rinke, “Benchmark of gw methods for az-abenzenes,” Phys. Rev. B86