The *GW* approximation to many-body perturbation theory is a reliable tool for describing charged electronic excitations, and it has been successfully applied to a wide range of extended systems for several decades using a plane-wave basis. However, the *GW* approximation has been used to test limited 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 the quasiparticle spectra of closed- and open-shell molecular anions with partially and completely filled 3*d* shells (shallow and deep 3*d* states, respectively), ScO^{−}, TiO^{−}, CuO^{−}, and ZnO^{−}, using various levels 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 find that the *G*-only eigenvalue self-consistent *GW* scheme with *W* fixed to the PBE level (*G*_{n}*W*_{0}@PBE), which gives the best compromise between accuracy and efficiency for solids, also gives good results for both localized (*d*) and delocalized (*sp*) states of 3*d*-transition-metal oxide molecules. The success of *G*_{n}*W*_{0}@PBE in predicting electronic excitations in these systems reasonably well is likely due to the fortuitous cancellation effect between the overscreening of the Coulomb interaction by PBE and the underscreening by the neglect of vertex corrections. Together with the absence of the self-consistent field convergence error (e.g., spin contamination in open-shell systems) and the *GW* multisolution issue, the *G*_{n}*W*_{0}@PBE scheme gives the possibility to predict the electronic structure of 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 the electronic structure of an interacting many-electron system. In experiment, electron removal and addition energies of both extended and finite systems are measured by direct and inverse photoelectron spectroscopy (PES and IPES, respectively). In theory, it is well known that the *GW* approximation to many-body perturbation theory (MBPT) describes bandgaps and band structures of solids more accurately than local and semilocal approximations to density-functional theory (DFT).^{1,2} However, less is known about the performance of the *GW* approximation on the electronic structure of atoms, molecules, and clusters. Especially, *GW* calculations for the quasiparticle (QP) spectra of open-shell molecules containing 3*d* transition metals are scarce. There are a few reasons for it.

First, it is easier to test only frontier orbital energies such as the ionization energy (IE) and the electron affinity (EA) than the full QP spectrum (all orbital energies). There are mainly two ways to calculate IE and EA of molecules. On the one hand, IE (EA) can be obtained from DFT, HF (Hartree-Fock), MP2 (second-order Møller-Plesset perturbation theory), RPA (random-phase approximation), or CCSD(T) (coupled-cluster singles and doubles plus perturbative triples) total energy differences between a neutral and a cation (anion) within the so-called ΔSCF (self-consistent field) method.^{3} Generally, the ΔSCF method gives good results for frontier orbital energies of finite systems but cannot be applied to extended systems. Also, it is not straightforward for the Δ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 molecular orbital) and the LUMO (lowest occupied molecular orbital), respectively. Due to the simplicity of the ΔSCF method, many studies have evaluated the performance of the *GW* approximation on molecules by comparing *GW* IE and EA to ΔSCF ones,^{4} but that approach does not utilize the full power of the *GW* approximation, which is the ability to provide the QP spectrum for both finite and extended systems. For example, Bethe-Salpeter equation (BSE) calculations for optical excitations require more orbital energies than IE and EA as input.^{5}

Second, it is easier to test closed-shell systems than open-shell ones. Most of the quantum chemistry-based *GW* implementations for finite systems, such as MOLGW,^{6} FIESTA,^{7} TURBOMOLE,^{8} FHI-AIMS,^{3} and CP2K,^{9} use local-orbital basis sets such as Gaussian basis sets. *GW* calculations require mean-field self-consistent calculations, such as restricted and unrestricted Hartree-Fock or Kohn-Sham (RHF or RKS and UHF or UKS, respectively) calculations for closed- and open-shell systems, respectively. The problem is that unlike RHF and RKS self-consistent calculations, UHF and UKS ones are not guaranteed to converge, as their convergence strongly depends on the initial guess wavefunctions. This is especially the case for spin-unrestricted calculations performed with hybrid exchange-correlation (xc) functionals, which include a fraction of exact exchange (EXX), and HF on 3*d*-transition-metal-containing molecules due to the near-degeneracy of energy levels.^{10–12} Partially due to this SCF convergence issue, most existing studies have used only closed-shell systems to assess the performance of the *GW* approximation on finite systems. For example, Refs. 13–16 used the so-called *GW*100 benchmark set, 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 accurately predict the electronic structure of *d* systems (especially 3*d* systems) than *sp* ones because of the strong localization, and thus the strong correlation, of *d* electrons. For example, it is challenging for *GW* to accurately reproduce the experimental bandgap and *d*-band position of bulk ZnO at the same time.^{17–19} Practically, it is computationally 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, which increases 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 the human effort by making it necessary to manually explore many minima with similar energies using many initial guess wavefunctions.^{10–12}

The *GW* approximation is unique, but due to its high computational costs, there are various *GW* schemes and variants. Generally, there are two approaches. One approach is to vary the self-consistency level in the *GW* approximation. The *GW* self-consistent levels from the lowest to the highest include the perturbative non-self-consistent (one-shot) *GW* (*G*_{0}*W*_{0}) scheme, the eigenvalue self-consistent *GW* (ev*GW*) scheme (with two types *G*_{n}*W*_{0} and *G*_{n}*W*_{n}, which update eigenvalues only in *G* and in both *G* and *W*, respectively), the QP self-consistent *GW* (QS*GW*) scheme using a static and Hermitian approximation to the *GW* self-energy, and the fully self-consistent *GW* (SC*GW*) scheme.^{1} Generally, as the *GW* self-consistency level increases, the *GW* approximation depends less on the mean-field starting point and becomes more conserving with respect to the particle number, momentum, and energy. However, the higher *GW* self-consistency level does not necessarily give more accurate QP energies because vertex corrections are missing in the *GW* approximation. For example, SC*GW* and QS*GW* systematically overestimate the bandgaps of solids,^{20} displaying worse performance than ev*GW*, which currently provides the best balance between accuracy and efficiency for solids.^{17}

The other approach is to vary the amount of EXX in the *GW* starting point to reduce the self-interaction error by (semi-)local xc functionals. Typically, the *G*_{0}*W*_{0} scheme chooses this approach to obtain good results at low computational costs. However, the predictive power of this approach is questionable since the optimal amount of EXX in the *GW* starting point is strongly system-dependent. For example, for extended systems, the reported values for the optimal amount of EXX are narrowly spread between 0% and 25%,^{18} while for finite systems, they are widely spread between 25% and 100%.^{4,21–23}

The purpose of this work is to evaluate the performance of the *GW* approximation on the electronic structure of small oxide molecules containing 3*d* transition metals. To this end, we calculate the QP spectra of closed- and open-shell molecular anions with partially and completely filled 3*d* shells, ScO^{−}, TiO^{−}, CuO^{−}, and ZnO^{−}, using various levels of *GW* theory. There are a few reasons why we chose these molecular systems: (i) their anion PES data are available,^{24–27} (ii) CuO^{−} and ZnO^{−} are molecular analogs to bulk Cu_{2}O and ZnO, respectively, which are challenging systems for the *GW* method,^{28,29} and (iii) shallow and deep 3*d* states are measured in TiO^{−} and CuO^{−}, respectively.

This article is organized as follows: First, we give a brief introduction to the *GW* approximation and its implementation in the framework of quantum chemistry. Second, we present various convergence test results and show that care should be taken to obtain reliable and reproducible QP energies of finite systems from Gaussian-based *GW* implementations. Third, we assess various *GW* schemes, focusing on ionization energies and 3*d*-electron binding energies, and conclude that the *G*_{n}*W*_{0}@PBE scheme gives the best performance among *GW* schemes considered in this work in terms of accuracy and efficiency. Last, we discuss the origin of seemingly conflicting *GW* results for finite systems in the literature.

## II. THEORETICAL BACKGROUND

In this section, we briefly review the *GW* approximation and its implementation using local-orbital basis sets. This section contains only a minimal number of equations, which will be needed later. More details can be found in Refs. 6, 8, 13, and 30. Generally, we follow the notation in the MOLGW implementation paper^{6} for consistency: (i) Hartree atomic units are used in all equations; (ii) the complex conjugate notation is not used for wavefunctions 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 referred to 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 nonlocal, dynamical, and non-Hermitian self-energy Σ^{σ} at frequency *ω* is given by

where *σ* is the spin channel (↑ or ↓), *G*^{σ} is the time-ordered one-particle Green’s function, *W* is the dynamically screened Coulomb interaction, and *η* is a positive infinitesimal.

The self-energy in Eq. (1) can be calculated from first principles by solving the coupled Hedin’s equations in order. One starts by constructing the one-particle Green’s function using the one-electron eigenvalues $\u03f5m\sigma $ and corresponding wavefunctions $\phi m\sigma (r)$ obtained from the Hartree or mean-field approximation,

where *i* runs over occupied states and *a* runs over empty states. Note that *G*^{σ} in Eq. (2) is not the interacting (dressed) Green’s function but the noninteracting (bare) one, which are conventionally denoted by *G*^{σ} and $G0\sigma $, respectively. In this work, we use the subscript 0 to distinguish the non-self-consistent *GW* scheme from the self-consistent one.

Using the one-particle Green’s function in Eq. (2), one can successively obtain the noninteracting (irreducible) polarizability *χ*_{0} and the interacting (reducible) polarizability $\chi =\chi 0[1\u2212v\chi 0]\u22121$ within the RPA, the screened Coulomb interaction, and the self-energy,

where $v$ denotes the bare (unscreened) Coulomb interaction $v$(**r**, **r**′) = 1/|**r** − **r**′|, Σ_{x} is the exchange part of the self-energy, and Σ_{c} is the correlation part of the self-energy. Note that in Eqs. (3)–(5), space and frequency variables (**r**, **r**′, *ω*) are omitted for simplicity, and *χ*_{0}(*ω*), *χ*(*ω*), Σ^{σ}(*ω*), and $\Sigma c\sigma (\omega )$ are dynamic, whereas $v$ and $\Sigma x\sigma $ are static. Note also that *W* is obtained without using the dielectric matrix.

Using the real part (Re) of the self-energy in Eq. (5) and the first-order perturbation theory, one can obtain the (diagonal) QP equation

where $\u03f5mG0W0,\sigma $ are the perturbative one-shot *GW* QP energies and $vxc\sigma $ is the xc potential. Experimentally, $\u03f5mG0W0,\sigma $ correspond to vertical IEs and EAs in PES and IPES, respectively. Theoretically, $\u03f5mG0W0,\sigma $ correspond to the positions of poles of Green’s function in the spectral (Lehmann) representation and thereby to the positions of QP peaks and plasmon satellites in the corresponding spectral function *A*^{σ},

where $Amm\sigma $ are the diagonal elements of the spectral function, $Gmm\sigma $ are the diagonal elements of Green’s function, and Im represents the imaginary part. Note that $Amm\sigma $ give the local density of states.

*G*^{σ} in Eq. (7) is the interacting Green’s function, whereas *G*^{σ} in Eq. (2) is the noninteracting one. In other words, by plugging *G*^{σ} in Eq. (2) into Eq. (7) after replacing $\u03f5m\sigma $, where *m* = *i* or *a*, by $\u03f5m\sigma +\u27e8\phi m\sigma |\Sigma \sigma (\omega )\u2212vxc\sigma |\phi m\sigma \u27e9$, one can find that $Amm\sigma $ have Lorentzian peaks at

which shows that solving the QP equation in Eq. (6) and locating the peak positions in the spectral function in Eq. (7) are equivalent ways of obtaining the QP energies.

The QP equation in Eq. (6) is nonlinear because Σ^{σ} depends on $\u03f5mG0W0,\sigma $, 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 nonlinear QP equation and to iteratively solve the coupled Hedin equations will be discussed later.

### B. Self-consistent field method

In order to obtain the ingredients for the one-particle Green’s function in Eq. (2) using local-orbital basis sets, molecular orbitals (MOs) and corresponding MO energies are used as one-electron wavefunctions and corresponding eigenvalues. MOs are expanded as a linear combination of atomic orbitals (AOs) *ϕ*_{μ},

where $C\mu m\sigma $ are MO expansion coefficients. In MOLGW, atom-centered (contracted) Gaussian orbitals are used as AOs.

The MO coefficients in Eq. (9) and MO energies can be obtained by solving the generalized Kohn-Sham (gKS) equation [the Hartree-Fock–Kohn-Sham scheme^{31} for (semi-)local functionals, hybrid functionals, and HF],

where **C**^{σ} is a matrix of MO coefficients, *ϵ*^{σ} is a diagonal matrix of MO energies, **S** is the AO overlap matrix with elements,

and **H**^{σ} is the Hamiltonian matrix with elements

where *T*, $V$_{ext}, *J*, and *K*^{σ} are the kinetic energy, external potential energy, Hartree, and Fock exchange terms, respectively, $Vx\sigma $ and $Vc\sigma $ are the exchange and correlation potentials, respectively, and *α* is the fraction of EXX in hybrid functionals that will be introduced later.

We briefly explain only a few terms in the Hamiltonian matrix in Eq. (12), which will be needed later. The matrix elements of the Hartree term in Eq. (12) are given by

where (*μν*|*λτ*) are the 4-center two-electron Coulomb repulsion integrals,

and **D**^{σ} is the density matrix with elements

where *f*^{σ} is the occupation number (0 or 1). The matrix elements of the Fock exchange term in Eq. (12) are given by

The exchange and correlation potentials in Eq. (12) depend on the density *ρ*^{σ} (and the density gradient ∇*ρ*^{σ}),

### C. *GW* self-energy

In order to obtain the ingredients for the interacting polarizability in Eq. (4), one should solve the Casida equation in matrix form,

where **A** and **B** are the resonant and coupling matrices, respectively, and Ω_{s} and (*X*^{s}, *Y*^{s}) are the eigenvalues (the neutral two-particle excitation energies) and corresponding eigenvectors, respectively. The matrix elements in **A** and **B** are given by

where *i* and *j* are for occupied states, *a* and *b* are for empty states, *f*_{xc} is the time-dependent density-functional theory (TDDFT) xc kernel, and (*iaσ*| *jbσ*′) are the 4-orbital two-electron Coulomb repulsion integrals,

In this work, we used the RPA by setting *f*_{xc} = 0. Note that Ref. 8 showed that TDDFT and RPA polarizabilities make a difference of ∼0.1 eV in *G*_{0}*W*_{0}@PBE QP HOMO energy. Note also that the dimension of the Casida matrix in Eq. (18) scales as *O*(*N*^{2}) with *N* being the system size, so building and completely diagonalizing the Casida matrix scale as *O*(*N*^{4}) and *O*(*N*^{6}), respectively. The MO integrals in Eq. (21) are transformed from the AO integrals in Eq. (14) through the AO-MO integral transformation,

which scales as *O*(*N*^{5}). Note that this integral transformation is a bottleneck in Gaussian-based *GW* and MP2 calculations.

Diagonalizing the Casida matrix in Eq. (18) yields eigenvalues Ω_{s} and eigenvectors (*X*^{s}, *Y*^{s}). In MOLGW, the diagonalization is performed without using the Tamm-Dancoff approximation (TDA), which sets **B** to zero, but efficiently using the so-called beyond-TDA method.^{6,32,33} Using Ω_{s} and (*X*^{s}, *Y*^{s}), one can construct the spectral representation of the interacting polarizability *χ*(*ω*).^{8,13,30} From *χ*(*ω*) and Eq. (4), one can obtain the spectral representation of the screened Coulomb interaction *W*(*ω*).^{6,8,13,30} Using *W*(*ω*) and Eq. (5), and analytically performing the convolution of *G*^{σ}(*ω*) and *W*(*ω*) in the frequency domain, one can obtain the exchange and correlation parts of the *GW* self-energy $\Sigma x\sigma $ and $\Sigma c\sigma (\omega )$, respectively, whose diagonal matrix elements are given by

where *i* runs over occupied states, *a* runs over empty states, *s* runs over all excitations, and $wmn\sigma s$ are given by

Note that unlike the plasmon-pole approximation (PPA), the analytic continuation method, and the contour deformation technique,^{13,34} the fully analytic method employed in RGWBS,^{30} TURBOMOLE, and MOLGW gives the exact *GW* self-energy at all frequency points because it does not rely on any approximation and numerical parameter.

### D. Spin contamination

In unrestricted HF and KS calculations for open-shell systems, the expectation value of the total angular momentum ⟨*S*^{2}⟩ is given by

where *N*_{↑} and *N*_{↓} are the numbers of ↑- and ↓-spin electrons, respectively, and *S* is (*N*_{↑} − *N*_{↓})/2 with *N*_{↑} > *N*_{↓}. The last two terms on the right side of Eq. (26) are called the spin contamination, which is non-negative.^{35,36} The spin contamination becomes large when a ground state is mixed with (contaminated by) excited states.

In restricted calculations for closed-shell systems, the SCF cycle always converges to a global minimum and the spin contamination is zero for all (semi-)local and hybrid functionals as well as HF. In unrestricted calculations for open-shell systems, the SCF convergence and the spin contamination depend on EXX amount and basis size. For (semi-)local functionals, the SCF cycle almost always converges to a global minimum and the spin contamination is small [generally smaller than ∼10% of *S*(*S* + 1)]. For hybrid functionals and HF, there is a chance (which increases with EXX amount and basis size) that the SCF cycle fails, does not converge, or converges to local minima or the spin contamination is large.

There are a few points to note about the spin contamination. First, the spin contamination is just an indicator for the SCF convergence error; therefore, a small spin contamination does not guarantee the correct SCF convergence. Second, the spin contamination generally raises but sometimes lowers the gKS total energy, so the lowest gKS total energy does not guarantee the correct SCF convergence, either. Last, the spin contamination and the SCF cycle are independent of each other. For example, the SCF cycle can converge quickly with large spin contamination or slowly with small spin contamination.

### E. Auxiliary basis sets and multithread parallelization

In Gaussian-based *GW*, the 4-center integrals (*μν*|*λτ*) in Eq. (14), which scale as *O*(*N*^{4}), are a common bottleneck in gKS and *GW* parts in terms of compute time and memory storage. One way to reduce the bottleneck is the resolution-of-identity (RI) approximation (the density-fitting approximation), which expands the product of basis functions *ϕ*_{μ}(**r**)*ϕ*_{ν}(**r**) as a linear combination of auxiliary basis functions *ϕ*_{P}(**r**).^{3,37,38} There are two types of RI approximation: RI-V using a Coulomb metric and RI-SVS using an overlap metric. For example, FIESTA uses both RI-V and RI-SVS, whereas MOLGW uses only RI-V, which is known to be superior to RI-SVS. Within RI-V, the 4-center integrals (*μν*|*λτ*) in Eq. (14) approximate to

where *P* and *Q* run over auxiliary basis functions, (*μν*|*P*) and (*Q*|*λτ*) are the 3-center integrals, and (*P*|*Q*) are the 2-center integrals.

RI can be applied to both gKS [*J* and *K*^{σ} in Eqs. (13) and (16)] and *GW* [**A**, **B**, $\Sigma x\sigma $, and $\Sigma c\sigma (\omega )$ in Eqs. (19), (20), and (22)–(25)] parts. In this work, we refer 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 reduces both compute time and memory storage by about the number of basis functions (by ∼100 times as shown in Table I).

. | . | . | . | $Nemp\u2191$ . | |||
---|---|---|---|---|---|---|---|

Anion . | Potential . | FC . | $Nocc\u2191$ . | CN = 2 . | CN = 3 . | CN = 4 . | CN = 5 . |

ScO^{−} | AE | Yes | 9 | 67 | 124 | 205 | 314 |

ScO^{−} | AE | No | 15 | 67 | 124 | 205 | 314 |

TiO^{−} | AE | Yes | 10 | 66 | 123 | 204 | 313 |

TiO^{−} | AE | No | 16 | 66 | 123 | 204 | 313 |

CuO^{−} | AE | Yes | 13 | 63 | 120 | 201 | 310 |

CuO^{−} | AE | No | 19 | 63 | 120 | 201 | 310 |

CuO^{−} | ECP | Yes | 13 | 63 | 120 | 201 | 310 |

CuO^{−} | ECP | No | 14 | 63 | 120 | 201 | 310 |

ZnO^{−} | AE | Yes | 14 | 62 | 119 | 200 | 309 |

ZnO^{−} | AE | No | 20 | 62 | 119 | 200 | 309 |

ZnO^{−} | ECP | Yes | 14 | 62 | 119 | 200 | 309 |

ZnO^{−} | ECP | No | 15 | 62 | 119 | 200 | 309 |

. | . | . | . | $Nemp\u2191$ . | |||
---|---|---|---|---|---|---|---|

Anion . | Potential . | FC . | $Nocc\u2191$ . | CN = 2 . | CN = 3 . | CN = 4 . | CN = 5 . |

ScO^{−} | AE | Yes | 9 | 67 | 124 | 205 | 314 |

ScO^{−} | AE | No | 15 | 67 | 124 | 205 | 314 |

TiO^{−} | AE | Yes | 10 | 66 | 123 | 204 | 313 |

TiO^{−} | AE | No | 16 | 66 | 123 | 204 | 313 |

CuO^{−} | AE | Yes | 13 | 63 | 120 | 201 | 310 |

CuO^{−} | AE | No | 19 | 63 | 120 | 201 | 310 |

CuO^{−} | ECP | Yes | 13 | 63 | 120 | 201 | 310 |

CuO^{−} | ECP | No | 14 | 63 | 120 | 201 | 310 |

ZnO^{−} | AE | Yes | 14 | 62 | 119 | 200 | 309 |

ZnO^{−} | AE | No | 20 | 62 | 119 | 200 | 309 |

ZnO^{−} | ECP | Yes | 14 | 62 | 119 | 200 | 309 |

ZnO^{−} | ECP | No | 15 | 62 | 119 | 200 | 309 |

RI is an approximation, so it causes an error. There are mixed results for the RI error in the literature, ranging from ∼1 meV to ∼0.1 eV, because different molecular systems, molecular orbitals, levels of theory (DFT vs *GW*), xc functionals (PBE vs HF), and basis sets are used to evaluate the quality of RI.^{13,21}

Another way to reduce the bottleneck without causing an error is the parallelization. We parallelized the 4-center integrals in Eqs. (13), (16), (19), (20), and (22)–(25) [as well as other bottlenecks, such as the integral transformation in Eq. (22) and the correlation part of the self-energy in Eq. (24)] using Open Multi-Processing (OpenMP), which consumes much less memory than the message passing interface by using shared-memory threads. The performance gain by our OpenMP parallelization is shown in the supplementary material. We also optimized our OpenMP implementation to reduce Nonuniform Memory Access (NUMA) effects in modern multicore processors by enhancing the memory bandwidth and reducing the memory latency. Our OpenMP implementation in MOLGW 1.F has recently been merged into MOLGW 2.A.

### F. *G*_{0}*W*_{0} quasiparticle energy

In this work, we used three methods to obtain $\u03f5mG0W0$, as is practically impossible to obtain unique, correct, and accurate $\u03f5mG0W0$ for all energy levels of all molecular systems using a single method, which will be discussed in detail 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 nonlinear QP equation in Eq. (6),

where $\u03f5mG0W0,linear$ is the perturbative one-shot QP energy obtained from the linearization, and $Zmlinear$ is the QP renormalization factor for the linearization,

where the derivative of the self-energy is obtained from the finite difference method using two frequency points at *ϵ*_{m} ± Δ*ω* with Δ*ω* being the frequency grid spacing, which is set to 0.001 hartree in this work.

There are a few points to note about the linearization method. First, one can choose different frequency points for the finite difference method (e.g., *ϵ*_{m} ± 0.1 eV and *ϵ*_{m} ± 0.5 eV in Refs. 15 and 39, respectively). In the PPA *G*_{0}*W*_{0} method, different frequency points give similar results for $\u03f5mG0W0,linear$ because PPA makes ⟨*φ*_{m}|Σ_{c}(*ω*)|*φ*_{m}⟩ in Eq. (24) smooth around *ϵ*_{m} by drastically reducing the number of self-energy poles.^{39} However, in the full-frequency *G*_{0}*W*_{0} method, different frequency points can give very different results for $\u03f5mG0W0,linear$ when the finite difference method fails due to weak self-energy poles around *ϵ*_{m}, which will be discussed later. Second, the derivative of the self-energy can be evaluated analytically,^{8} but it is not used in this work. Last, $Zmlinear$ in Eq. (29) is slightly different from that in Ref. 13,

The derivative of the self-energy is evaluated at *ω* = *ϵ*_{m} and $\omega =\u03f5mG0W0$ in $Zmlinear$ and *Z*_{m}, respectively. Generally, *Z*_{m} is smaller than $Zmlinear$ because the self-energy has a steeper slope at $\omega =\u03f5mG0W0$ than at *ω* = *ϵ*_{m}. *Z*_{m} represents the QP weight (the pole residue of Green’s function), which equals the area under the Lorentzian QP peak and depends on the spectral weight transfer from the QP peak to plasmon satellites and the incoherent background.

The second method is to graphically solve the nonlinear QP equation in Eq. (6) using the secant (quasi-Newton) method. In this work, we refer to $\u03f5mG0W0$ obtained from the graphical solution as $\u03f5mG0W0,graph$. Note that the above linearization method corresponds to the first step of the secant method. Note also that when the nonlinear QP equation in Eq. (6) has multiple solutions, the secant method gives only one of them, which depends strongly on the choice of *η*.

The last method is to use the position of the QP peak with the highest spectral weight in the spectral function *A*_{mm} in Eq. (7). In this work, we refer to $\u03f5mG0W0$ obtained from the spectral function as $\u03f5mG0W0,spect$ and define $Zmspect$ by replacing $Zmlinear$ and $\u03f5mG0W0,linear$ in Eq. (28) by $Zmspect$ and $\u03f5mG0W0,spect$, respectively. Note that we searched for the QP peak at $0<Zmspect<1$ using the peak height instead of the spectral weight (the area under the peak) due to the practical difficulty of determining the peak range. Note also that the highest spectral weight gives the largest *Z*_{m} in Eq. (30) because *Z*_{m} represents the spectral weight, as explained above.

### G. *G*_{n}*W*_{0} and *G*_{n}*W*_{n} quasiparticle energy

As introduced in Sec. I, there are various levels of self-consistency in the *GW* approximation (from the lowest to the highest): *G*_{0}*W*_{0}, *G*_{n}*W*_{0}, *G*_{n}*W*_{n}, QS*GW*, and SC*GW*. In this work, we used *G*_{n}*W*_{0} and *G*_{n}*W*_{n} for simplicity, efficiency, and stability. *G*_{n}*W*_{0} updates only gKS eigenvalues in *G*^{σ} [$\u03f5i\sigma $ and $\u03f5a\sigma $ in Eq. (24)], whereas *G*_{n}*W*_{n} updates gKS eigenvalues in *G*^{σ} and **A** [$\u03f5i\sigma $ and $\u03f5a\sigma $ in Eq. (19)] as well as Casida eigenvalues in *W* [Ω_{s} in Eq. (24)]. Therefore, *G*_{n}*W*_{n} is computationally more expensive than *G*_{n}*W*_{0} by the time to build and completely diagonalize the RPA Casida matrix in Eq. (18). Note that *G*_{n}*W*_{n} can be viewed as a diagonal approximation to QS*GW*.

In this work, we obtained *G*_{n}*W*_{0} and *G*_{n}*W*_{n} QP energies ($\u03f5mGnW0$ and $\u03f5mGnWn$, respectively) by iterating the recurrence relations (*n* ≥ 3),

where $\u03f5mevGW,n$ is $\u03f5mGnW0$ or $\u03f5mGnWn$, and *Z*^{evGW} = 1. Summing up Eqs. (31)–(33), we get

which we refer to as the ev*GW* QP equation in this work. When the ev*GW* convergence is reached ($\u03f5mevGW,n=\u03f5mevGW,n\u22121=\u03f5mevGW,\u221e$, where $\u03f5mevGW,\u221e$ are converged ev*GW* QP energies), the ev*GW* QP equation in Eq. (34) becomes similar to the *G*_{0}*W*_{0} QP equation in Eq. (6).

Whereas most *GW* codes use 0 < *Z*^{evGW} < 1,^{7,17,40} MOLGW uses *Z*^{evGW} = 1. Even though we implemented ev*GW* with 0 < *Z* < 1 into MOLGW, we adopted ev*GW* with *Z* = 1 in this work for a few reasons. First, Eq. (33) shows that converged ev*GW* QP energies ($\u03f5mevGW,\u221e$) are independent of whether 0 < *Z* < 1 or *Z* = 1. Second, *Z* = 1 gives a unique solution that satisfies the QP equation in Eq. (34), which allows us to avoid the *GW* multisolution issue from the graphical-solution and spectral-function methods and the ∼0.1–1 eV error from the linearization method (to be discussed in detail later). Third, ev*GW* with 0 < *Z* < 1 is suited for a simplified ev*GW* variant that updates only a few states near HOMO and LUMO and rigidly shifts all the other states for efficiency,^{17,41} but we updated all eigenvalues in this work for accuracy. Last, ev*GW* with *Z* = 1 has no variant and does not need a QP equation solver, but ev*GW* with 0 < *Z* < 1 has multiple variants, depending on the choice of QP equation solvers. For example, two ev*GW* variants with 0 < *Z* < 1 using the linearization and graphical-solution methods in Refs. 7 and 40, respectively, may give different QP energies because the two 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*_{0}*W*_{0} are discussed in the supplementary material.

## III. COMPUTATIONAL DETAILS AND TEST RESULTS

### A. Computational details

Our gKS calculations were carried out using both MOLGW and NWChem in order to cross-check the results and to ascertain the correct SCF convergence. For *GW* calculations, we used only MOLGW. MOs were expanded using augmented Dunning correlation-consistent Gaussian basis sets, aug-cc-pV*n*Z (*n* = D, T, Q, and 5), which are designed to smoothly converge with basis size. Augmentation using diffuse functions is essential in ground-state calculations for anions and in excited-state calculations for both neutrals and anions. Without augmentation, gKS and *GW* eigenvalues for empty states converge very slowly with basis size.^{9,21} In the following, the cardinal number (CN = 2, 3, 4, and 5) is used to represent the approximate size of diverse basis sets employed in the literature and this work. For example, CN = 4 means def2-QZVP in Ref. 13 and aug-cc-pVQZ in this work. Table I summarizes the exact size of CN = 2,3,4,5 basis sets used in this work. To determine the optimized bond lengths of TMO anions, we used NWChem with PBE and CN = 3. We obtained bond lengths of 1.695, 1.642, 1.697, and 1.765 Å for ScO^{−}, TiO^{−}, CuO^{−}, and ZnO^{−}, respectively.

In order to study the starting-point dependency of the *GW* approximation, we used global hybrid functionals,

where $ExHF,\sigma $, $ExPBE,\sigma $, and $EcPBE,\sigma $ are Fock exact exchange, PBE exchange, and PBE correlation energies, respectively. We refer to the hybrid functionals in Eq. (35) as PBE*α* functionals in this work. While we tested other functionals such as B3LYP, HSE06, BHLYP, and HF, we discuss only PBE*α*(0.00 ≤ *α* ≤ 1.00) results because the EXX amount in the starting point has a stronger effect on *GW* results than other factors such as range separation (to screen the Coulomb interaction) and correlation type. As shown in the supplementary material, HSE06, PBE0, and BHLYP*α*(*α* = 0.25) [where PBE is replaced by LYP in Eq. (35)] give similar *GW* results. Note that the type of the correlation functional is not important (e.g., PBE vs LYP), but the existence of it is. As shown in the supplementary material, PBE*α*(*α* = 1.00) and HF can make a large difference (∼1 eV) in *GW* results for some states.

### B. gKS test results

#### 1. Effective Core Potentials

Unlike Sc and Ti, Cu and Zn have two choices of basis sets: AE (All Electron) and ECP (Effective Core Potential). ECP allows us to remove core electrons and include relativistic effects. We first tested scalar relativistic effects 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 relatively light elements.^{42,43} The test results are presented in the supplementary material. We found that (i) ECP and AE *GW* IEs differ by 0.01–0.15 eV, depending on the subtle competition between direct and indirect relativistic effects (*s* and *p* orbital contraction and stabilization and *d* and *f* orbital expansion and destabilization, respectively),^{44–46} which is consistent with Ref. 47, and (ii) ECP *GW* 3*d*-electron binding energies are smaller than AE ones by 0.16–0.66 eV due to indirect effects, which is consistent with Ref. 48. We next tested the efficiency of ECP with respect to AE. We found that ECP is more efficient than AE because the absence of core states not only makes the basis size smaller, which benefits both gKS and *GW* parts, but also makes the SCF cycle faster and more stable, which benefits only the gKS part. In this work, we present mainly AE results not because AE is superior to ECP but because scalar relativistic effects make smaller changes than large (∼1–2 eV) errors that we encounter. However, we discuss both AE and ECP results for 3*d*-electron binding energies, where scalar relativistic effects are considerable (∼10% of the experimental 3*d*-electron binding energy). Note that Ref. 13 used only AE even though the *GW*100 benchmark set contains Ag_{2}, Cu_{2}, and NCCu molecules.

#### 2. RI for gKS

We did not use RI in this work because our goal is to assess the range of applicability of the *GW* approximation as accurately as possible using small molecules, but RI is unavoidable for the practical *GW* study of large molecules. Thus, we evaluated the quality of RI for both AE and ECP by comparing RI and no-RI gKS eigenvalues and total spins. The evaluation results are presented in the supplementary material. We found that CN = 5 RI ECP causes a large random error in gKS results (e.g., ∼0.2 and ∼0.8 eV for CuO^{−} and ZnO^{−}, respectively, in gKS-PBE IEs), which *decreases* with the EXX amount. It is important to note that unlike the SCF convergence error, which occurs only in open-shell systems with nonzero EXX amounts, this gKS RI error occurs in both closed- and open-shell systems with all EXX amounts. It is difficult to detect the gKS RI error because all SCF cycles with different convergence parameters smoothly converge to the same local minimum with no or small spin contamination. Therefore, we conclude that RI should be used only after the potential gKS RI error is thoroughly tested.

Note that because we did not use RI and the 4-center integrals are computed at each SCF step, a *single* gKS calculation is as expensive as a *single GW* calculation in this work, which is consistent with Ref. 9. Note also that the effect of RI on *GW* results is discussed in the supplementary material.

#### 3. SCF convergence tests

It is not straightforward to obtain the correct mean-field input for *GW* calculations because successful SCF convergence could come from both correct convergence to a global minimum and wrong convergence to some local minima. This is a particularly critical issue in gKS calculations on open-shell systems involving nonzero EXX and large basis. Many minima with similar total energies and total spins due to nearly degenerate energy levels in 3*d* transition metals make it more difficult to obtain correct SCF convergence.^{10–12} For closed-shell systems or open-shell systems with (semi-)local xc functionals, the SCF cycle is generally guaranteed to converge to a global minimum. However, when EXX is used for open-shell systems, wrong SCF convergence occurs frequently and randomly, which makes manual, time-consuming, and error-prone SCF convergence tests mandatory.

In order to obtain the correct mean-field input, we performed three-step SCF convergence tests. First, we used 12 and 96 sets of SCF convergence parameters for MOLGW and NWChem, respectively. Second, we manually searched for correctly converged SCF results using multiple indicators: gKS total energy, total spin ⟨*S*^{2}⟩ in Eq. (26), the number of total SCF cycles, a trend over basis size (CN = 2, 3, 4, 5), and a trend over EXX amount (by manually choosing gKS total energies and total spins that vary smoothly with basis size and EXX amount). Last, we cross-checked all MOLGW and NWChem gKS results. Our SCF convergence test results are presented in the 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. Accordingly, one should extrapolate *GW* QP energies obtained from different basis sizes to the complete basis set (CBS) limit to avoid the incomplete basis set error of ∼0.1 eV.^{3} We, therefore, tested the effect of the fitting function type, EXX amount, and basis size on the *GW* CBS limit.

Two fitting functions are most widely used for the CBS limit,^{49} which we refer to as standard fitting functions in this work,

where *E*_{m} are correlation or *m*th QP energies, *a* and *b* are fitting 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 CBS limit. Note that there are various nonstandard fitting functions used in the literature.^{6,50–53}

Figure 1 compares CBS results obtained from two standard fitting functions in Eqs. (36) and (37) (as well as one nonstandard one used in Refs. 6 and 53). We see that different fitting functions always give different *GW* CBS limits, deviating from each other by up to ∼0.1 eV depending on molecular systems and molecular orbitals. Figure 2 shows the effect of the EXX amount on the *GW* CBS limit. We observe that the incomplete basis set error increases with the EXX amount. CN = 2 occasionally and randomly causes a significant error (∼0.1 eV) in the *GW* CBS limit, which is commonly observed in the literature.^{8,13} Based on these test results, we conclude that it is important to check whether extrapolation is used or not, whether CN = 2 is used or not for extrapolation, and which fitting function is used when analyzing and comparing *GW* results. For example, Ref. 13 reported that IEs obtained from Gaussian- and plane-wave (PW)-based *GW* implementations with and without extrapolation, respectively, differ by ∼0.2 eV, but Refs. 15 and 16 showed that the use of PW *GW* IEs with extrapolation reduces the difference to ∼0.06 eV.

#### 2. Number of empty states

By enabling MOLGW 1.F to support the largest available basis set (CN = 5), we also tested the effect of CN = 5 on the *GW* CBS limit. The test results are presented in the supplementary material. Here, we briefly mention a 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 some cases, CN = 5 has an appreciable (∼0.1 eV) effect on the *GW* CBS limit by reducing the effect of the large random 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 error. Moreover, CN = 5 calculations are expensive (due to the large number of empty states and the slow SCF convergence speed) and error-prone (due to the high chance of SCF convergence and gKS RI errors). Therefore, we conclude that it is more beneficial to obtain the *GW* CBS limit from CN = 3,4 than from CN = 2,3,4,5. Using CN = 4 (∼100 empty states per atom, as shown in Table I) instead of CN = 5 as the largest basis set for the *GW* CBS limit tremendously reduces the computational costs. This conclusion is consistent with Ref. 13, which used only CN = 3,4 for extrapolation, and gives a rough estimate for two important and interdependent convergence parameters for $\Sigma c\sigma (\omega )$, the dimension of the dielectric matrix, and the number of empty states, in sum-over-states PW *GW* implementations.^{13}

The above conclusion holds only for occupied states. The effect of CN = 5 on the *GW* CBS limit for empty states is discussed in the supplementary material. The effect of the number of occupied states on *GW* results using the frozen-core (FC) approximation, which reduces the number of occupied states used in the construction of *G* and *W* and thus speeds up *GW* calculations,^{4} is also discussed in the supplementary material.

#### 3. *G*_{0}*W*_{0} quasiparticle energy

A full-frequency *G*_{0}*W*_{0} method used in this work produces complicated self-energy pole (and spectral-function peak) structures at nonfrontier orbitals, so it is not straightforward to automatically obtain correct and accurate $\u03f5mG0W0$ by using a single value of *η* and a single QP equation solver. Thus, we used three values of *η* (0.001, 0.002, and 0.005 hartree with Δ*ω* = 0.001 hartree) and three QP equation solvers mentioned in Sec. II F (linearization, graphical-solution, and spectral-function methods).

Our analysis of a total of nine solutions for $\u03f5mG0W0$ shows that the graphical-solution method using *η* = 0.001 hartree and the linearization method randomly give incorrect solutions. Therefore, we obtained $\u03f5mG0W0$ from the graphical-solution or spectral-function method using *η* = 0.002 or 0.005 hartree. When $\u03f5mG0W0,graph$ and $\u03f5mG0W0,spect$ are different, we manually selected a correct solution by analyzing Σ_{c}(*ω*) and *A*(*ω*).

A large distance between *ϵ*_{m} and $\u03f5mG0W0$ and multiple self-energy poles between them make it difficult to obtain correct and accurate $\u03f5mG0W0$ for nonfrontier orbitals, which is especially the case for (semi-)local xc functionals. Figure 3 shows various examples of successes and failures of three QP equation solvers for *G*_{0}*W*_{0}@PBE. In the following, we analyze each example individually. Note that we used *η* = 0.0001, 0.0002, and 0.0005 hartree with Δ*ω* = 0.0001 hartree in a few of the following examples to demonstrate the danger of small *η* values for the spectral-function method. Such small *η* values are not recommended, as they significantly increase disk storage requirements while barely improving accuracy.

First, the top two left panels of Fig. 3 show a general example, in which all three methods succeed. We see a few general trends. Typically, all three methods give correct solutions at *m* = HOMO and LUMO, which have a simple pole structure in Σ_{c}(*ω*). Graphical-solution and spectral-function methods generally give *multiple* solutions, whereas the linearization method always gives a *unique* solution, at which two straight lines intersect. Generally, graphical-solution and spectral-function methods give identical (correct and accurate) solutions (in this case, at *ω* = −0.48 eV), whereas the linearization method gives a different (correct but inaccurate) solution (in this case, at *ω* = −0.69 eV) due to an intrinsic error of ∼0.1–1 eV. A very small *η* (0.0002 hartree) sharpens a weak self-energy pole at *ω* = −1.0 eV in Σ_{c}(*ω*), but the sharpened pole *does not* cause an error in the graphical-solution method because it is not between $\u03f5mG0W0$ and *ϵ*_{m}. The very small *η* heightens the weak peak C in *A*(*ω*), but the heightened peak *cannot* cause an error in the spectral-function method because it is still lower than 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}(*ϵ*_{m} ± Δ*ω*) and (ii) *ω* = −1.0 eV is too distant from *ϵ*_{m} to affect the finite difference method.

Second, the top two right panels of Fig. 3 show a special example, in which the graphical-solution method can give an incorrect solution. We see a few special trends. Generally, some of the three methods give incorrect solutions at *m* = HOMO-*n* and LUMO+*n* (*n* = 1, 2, 3, …), which have a complicated pole structure in Σ_{c}(*ω*). A very small *η* (0.0002 hartree) sharpens a weak pole at *ω* = −2.8 eV in Σ_{c}(*ω*), and this sharpened pole causes a large error of 0.4 eV in the graphical-solution method because it is between $\u03f5mG0W0$ and *ϵ*_{m} and sharpened enough to make the secant method fail by causing it to find an incorrect intersection point. The very small *η* heightens a weak peak B in *A*(*ω*), but the heightened peak *cannot* cause an error in the spectral-function method because it is still lower than the other peak A. The very small *η* has little effect on the linearization method because *ω* = −2.8 eV is distant from *ϵ*_{m}.

Third, the bottom two left panels of Fig. 3 show a special example, in which the linearization method can give an incorrect solution. We see a few special trends. A very small *η* (0.0002 hartree) sharpens a weak pole at *ω* ≈ *ϵ*_{m} in Σ_{c}(*ω*), but the sharpened pole *does not* cause an error in the graphical-solution method even though it is between $\u03f5mG0W0$ and *ϵ*_{m} because the pole is not sharpened enough (*η* = 0.0001 hartree, on the other hand, causes a large error of 1.3 eV). The very small *η* heightens a weak peak at *ω* ≈ *ϵ*_{m} in *A*(*ω*), but the heightened peak *cannot* cause an error in the spectral-function method because it is still lower than the other peak at $\omega \u2248\u03f5mG0W0$. The very small *η* sharpens a weak pole at *ω* ≈ *ϵ*_{m} in Σ_{c}(*ω*), and the sharpened pole causes a large error of 2.1 eV in the linearization method because it is too close to *ϵ*_{m}, making the finite difference method fail by causing a large error in the slope of the tangent line at *ω* = *ϵ*_{m}.

Last, the bottom two right panels of Fig. 3 show a special example, in which the spectral-function method can give an incorrect solution. We see a few special trends. A very small *η* (0.0002 hartree) sharpens a weak pole at *ω* = −2.3 eV in Σ_{c}(*ω*), but the sharpened pole *does not* cause an error in the graphical-solution method because it is not between $\u03f5mG0W0$ and *ϵ*_{m}. Two peaks A and B (at *ω* = −3.9 and −2.1 eV, respectively) in *A*(*ω*) have similar spectral weights (and peak heights), so it is not straightforward to unambiguously determine which one is a QP peak or a satellite. Spectral weights (practically, peak heights) of the two peaks depend on the basis size: peak B (A) is higher than peak A (B) for CN = 2, 3, 4 (CN = 5). We *chose* peak B as the QP peak because (i) it is consistent with the solution from the graphical-solution method and (ii) it is consistent with a trend over EXX amount [*G*_{0}*W*_{0}@PBE*α*(0.00 ≤ *α* ≤ 1.00) QP HOMO-2 energies of CuO^{−} using the solution from peak B vary smoothly with *α*, as shown in Fig. 6]. The choice of *η* and CN has little effect on the linearization method, but $\u03f5mG0W0,linear=\u22122.9$ eV causes a large overestimation error of 0.8 eV in *G*_{0}*W*_{0}@PBE binding energy.

There are a few points to note about the above examples: (i) we chose simple examples, in which only one method can give an incorrect solution, for demonstration purposes; multiple methods can give incorrect solutions simultaneously, as shown in the supplementary material, (ii) not only a very small *η* but also a very large *η* (e.g., ∼0.05 hartree in Ref. 54) can cause a large error, (iii) deep states (e.g., HOMO-*n*, where *n* = 5, 6, …) have much more complicated pole [peak] structures in Σ_{c}(*ω*) [*A*(*ω*)] than those in Fig. 3, so it is very difficult to choose correct and accurate $\u03f5mG0W0$ for deep states not only automatically but also manually.

We conclude this section by summarizing several guidelines to obtain a reliable and reproducible *G*_{0}*W*_{0}@PBE QP spectrum. First, one should try multiple *η* (and Δ*ω*) values. There is no single general *η* value that works well for all QP equation solvers, molecular systems, and molecular orbitals. In other words, while *η* is typically viewed as a convergence parameter (the smaller the *η*, the more accurate the *GW* QP energy), it is practically an adjustable parameter, which should be not too small or too large (e.g., 0 and ∼0.05 hartree in Refs. 40 and 54, respectively). The optimal value of *η* depends on $|\u03f5mG0W0$-*ϵ*_{m}|, which generally decreases with the amount of EXX and increases with the depth of the *m*th state. For example, when calculating *G*_{0}*W*_{0}@PBE HOMO and LOMO (lowest occupied molecular orbital) energies, one may try ∼0.1 and ∼1 eV, respectively, for *η*.

Second, we recommend using multiple QP equation solvers. As shown in Fig. 3, the *G*_{0}*W*_{0}@PBE QP spectrum automatically obtained from a single QP equation solver can contain a large (∼1 eV) error at random states.

Third, we recommend using multiple basis sizes. As shown in the bottom right of Fig. 3, different basis sets with different sizes can give very different *G*_{0}*W*_{0} QP energies (by ∼1 eV) at random states. Using multiple basis sizes allows for not only accurate *GW* results without small (∼0.1 eV) systematic errors from the basis set incompleteness but also correct gKS and *GW* results without large (∼1 eV) random errors from SCF convergence and *GW* multisolution issues, respectively.

Fourth, one should be fully aware of the large random errors that the linearization method, which is the most widely used QP equation solver, can cause. Reference 15 suggests the linearization method as a preferable method for a fair comparison of *G*_{0}*W*_{0}@PBE IE (and EA) from different *GW* implementations because it gives a unique solution and thus is free of the *GW* multisolution issue. The idea works well for the IE, but it does not perform as well for the QP spectrum. For HOMO (and LUMO), the linearization method generally succeeds and systemically overestimates the IE only by ∼0.1 eV with respect to the accurate one from the graphical-solution and spectral-function methods, as shown in the top left of Fig. 3, accidentally reducing the ∼0.5 eV underestimation error by *G*_{0}*W*_{0}@PBE with respect to experiment.^{13,15,16} However, for deep states, it randomly succeeds or fails, as shown in the bottom left of Fig. 3, and randomly overestimates or underestimates *G*_{0}*W*_{0}@PBE binding energies by ∼1 eV compared to accurate ones, as shown in the bottom right of Fig. 3 and supplementary material, respectively. This large and unpredictable (with respect to state, magnitude, and direction) error makes the linearization method inadequate for the *G*_{0}*W*_{0}@PBE QP spectrum.

Last, one should be aware that different ways to handle the *GW* multisolution issue are found in the literature. For example, Ref. 13, which suggests that all solutions are physically relevant, manually searched for an actually relevant one by varying *η*, whereas Ref. 40 avoided the issue by automatically selecting the solution with the largest *Z*_{m} in Eq. (30) and using only *η* = 0. In this work, we adopted a combined approach. In other words, we automatically chose the solution with the highest spectral weight, which is identical to the solution with the largest *Z*_{m}, as explained in Sec. II F, when one solution is clearly more relevant than others, but we manually selected one solution that gives smoothly varying *GW* binding energies with a change in the *G*_{0}*W*_{0} starting point and ev*GW* self-consistency level without causing unphysical kinks when multiple solutions are equally relevant [e.g., two solutions at peaks A and B in the bottom two right panels of Fig. 3 give similar *Z*_{m} (∼0.2) for *m* = HOMO-2 of CuO^{−} due to similar slopes of the self-energy (∼ −4.0)]. However, when it comes to *η*, we adopted the approach of Ref. 13 instead of that of Ref. 40 because *η* = 0 frequently causes the secant method in the graphical-solution method to find an incorrect intersection point and makes $\u27e8\phi m\sigma |\Sigma c\sigma (\omega )|\phi m\sigma \u27e9$ in Eq. (24) diverge at $\omega =\u03f5i\sigma \u2212\Omega s$ and $\omega =\u03f5a\sigma +\Omega s$. The cumulant expansion,^{55,56} which describes plasmon satellites better than the *GW* approximation, may allow us to address the *GW* multisolution issue when the QP picture breaks down, but it is beyond the scope of this work.

#### 4. *G*_{n}*W*_{0} and *G*_{n}*W*_{n} quasiparticle energy

In this work, we used only *η* = 0.001 hartree for ev*GW* because it is small enough to obtain the convergence of ev*GW* QP energies with respect to *η* within ∼0.01 eV. The convergence test results for ev*GW* QP energies with respect to the iteration number are shown in the supplementary 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 of QS*GW* and our ev*GW* and found a couple of similarities and differences between them.

First, the ev*GW* convergence is reached after only a few iterations, which is consistent with the literature.^{17,41,51} Due to the fast and stable convergence, a mixing scheme is not used in our ev*GW*. Unlike ev*GW*, QS*GW* generally needs ∼10–20 (up to 60) iterations and a mixing scheme.^{22,57,58}

Second, the orbital character affects the starting-point dependency of ev*GW*. For example, we observed that ev*GW* QP energies for HOMO of CuO^{−} depend more strongly on the EXX amount in the ev*GW* starting point than those for HOMO of ScO^{−}. We attribute this to different amounts of 3*d* character in HOMOs of ScO^{−} and CuO^{−} (6% and 23%, respectively, as shown in Table II). In other words, as the 3*d* character in MO increases, the starting-point dependency of ev*GW* increases. We also observed that *G*_{n}*W*_{n} has a weaker (stronger) starting-point dependency for HOMO of ScO^{−} (CuO^{−}) than *G*_{n}*W*_{0}. Our observations for ScO^{−} are consistent with Ref. 41, which studied the ev*GW* starting-point dependency using small water clusters and concluded that as the ev*GW* self-consistency level increases from *G*_{n}*W*_{0} and *G*_{n}*W*_{n}, the ev*GW* starting-point dependency decreases. However, our observations for CuO^{−} are not consistent with this conclusion. This is likely because CuO^{−} 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 be related to conflicting results for QS*GW* in the literature: Ref. 58 showed the starting-point independency of QS*GW* using a small *sp*−bonded molecule (CH_{4}), while Ref. 59 showed the strong starting-point dependency of QS*GW* using a *d* solid (*α*-Fe_{2}O_{3}).

. | . | TiO^{−}
. | . | ZnO^{−}
. | ||
---|---|---|---|---|---|---|

. | ScO^{−}
. | ↑ . | ↓ . | CuO^{−}
. | ↑ . | ↓ . |

PBE | ||||||

HOMO | 0.06 | 1.00 | 0.00 | 0.23 | 0.00 | 0.00 |

HOMO-1 | 0.20 | 0.09 | 0.17 | 0.32 | 0.00 | 0.00 |

HOMO-2 | 0.18 | 0.21 | 0.20 | 1.00 | 0.05 | |

HOMO-3 | 0.24 | 0.75 | ||||

HOMO-4 | 0.48 | |||||

PBEα(α = 1.00) | ||||||

HOMO | 0.00 | 1.00 | 0.00 | 0.10 | 0.00 | 0.00 |

HOMO-1 | 0.13 | 0.07 | 0.10 | 0.05 | 0.00 | 0.00 |

HOMO-2 | 0.15 | 0.14 | 0.14 | 1.00 | 0.06 | |

HOMO-3 | 0.22 | 0.94 | ||||

HOMO-4 | 0.80 |

. | . | TiO^{−}
. | . | ZnO^{−}
. | ||
---|---|---|---|---|---|---|

. | ScO^{−}
. | ↑ . | ↓ . | CuO^{−}
. | ↑ . | ↓ . |

PBE | ||||||

HOMO | 0.06 | 1.00 | 0.00 | 0.23 | 0.00 | 0.00 |

HOMO-1 | 0.20 | 0.09 | 0.17 | 0.32 | 0.00 | 0.00 |

HOMO-2 | 0.18 | 0.21 | 0.20 | 1.00 | 0.05 | |

HOMO-3 | 0.24 | 0.75 | ||||

HOMO-4 | 0.48 | |||||

PBEα(α = 1.00) | ||||||

HOMO | 0.00 | 1.00 | 0.00 | 0.10 | 0.00 | 0.00 |

HOMO-1 | 0.13 | 0.07 | 0.10 | 0.05 | 0.00 | 0.00 |

HOMO-2 | 0.15 | 0.14 | 0.14 | 1.00 | 0.06 | |

HOMO-3 | 0.22 | 0.94 | ||||

HOMO-4 | 0.80 |

## IV. RESULTS AND DISCUSSION

In this section, we compare our *GW* calculations to anion PES experiments,^{24–27} focusing especially on the first IE, the lowest 3*d*-electron binding energy, and the orbital order. We present our results from two approaches separately: First, we discuss non-self-consistent *GW* with different starting-points (namely, *G*_{0}*W*_{0}@PBE*α* calculations as *α* is varied in steps of 0.25 from 0 to 1), and then, we discuss eigenvalue self-consistent *GW* (*G*_{n}*W*_{0} and *G*_{n}*W*_{n}) with the PBE starting point. We only briefly discuss our *GW* results for the starting-point–self-consistency hybrid approach because (i), fundamentally, we found that the hybrid approach does not give any better results than the two separate approaches and (ii), practically, the hybrid approach inherits disadvantages from both approaches.

ScO^{−}, TiO^{−}, CuO^{−}, and ZnO^{−} are similar but different systems in several aspects. First, ScO^{−} and CuO^{−} are closed-shell systems, whereas TiO^{−} and ZnO^{−} are open-shell systems. Second, ScO^{−} and TiO^{−} have partially filled 3*d* shells, while CuO^{−} and ZnO^{−} have completely 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 in TiO^{−} 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^{−}, but not in TiO^{−}, CuO^{−}, and ZnO^{−}. Due to these similarities and differences, TMO anions are an ideal set of systems for assessment of the performance of *GW* schemes.

Table II shows the amount of TM 3*d* character in all molecular orbitals considered in this work, obtained from CN = 2 gKS-PBE and gKS-PBE*α*(*α* = 1.00) using the Mulliken population analysis. We see that gKS-PBE ↑-HOMO of TiO^{−} and gKS-PBE HOMO-2 of CuO^{−} have entirely TM 3*d* character. Figure 4 shows the contour plots of all molecular orbitals considered in this work, obtained from CN = 2 gKS-PBE. It is clearly seen that ↑-HOMO of TiO^{−} and HOMO-2 of CuO^{−} are strongly localized on Ti and Cu, respectively. Table III summarizes our *GW* calculations with comparison to PES experiments and existing calculations in the literature.

. | . | . | G_{0}W_{0}@PBEα
. | G_{n}W_{0}@
. | G_{n}W_{n}@
. | . | ||||
---|---|---|---|---|---|---|---|---|---|---|

. | State . | Exp. . | α = 0.00
. | α = 0.25
. | α = 0.50
. | α = 0.75
. | α = 1.00
. | PBE . | PBE . | Others . |

ScO^{−} (^{1}Σ^{+}) | ||||||||||

HOMO | ^{1}Σ^{+} | 1.35^{a} | 0.51 | 1.15 | 1.45 | 1.59 | 1.63 | 1.21 | 1.35 | 1.28^{b}, 1.26^{c}, 1.19^{d} |

HOMO-1 | ^{2}Δ | 3.10^{a} | 3.30 | 4.77 | 5.45 | 5.72 | 5.82 | 5.34 | 6.08 | 2.41^{b}, 2.78^{e}, 3.31^{f} |

HOMO-2 | ^{2}Π | 3.40^{a} | 3.42 | 4.81 | 5.40 | 5.61 | 5.63 | 5.39 | 6.17 | 3.34^{b}, 3.24^{e}, 3.44^{f} |

MAE | 0.84^{g} | 0.18^{g} | 0.10^{g} | 0.24^{g} | 0.28^{g} | 0.14^{g} | 0.00^{g} | |||

TiO^{−} (^{2}Δ) | ||||||||||

↑-HOMO | ^{1}Σ^{+} | 2.00^{h} | 0.26 | 2.24 | 3.70 | 4.79 | 5.63 | 1.83 | 2.74 | 2.39^{b}, 2.37^{e}, 2.34^{f} |

↑-HOMO-1 | ^{1}Δ | 1.73^{h} | 0.53 | 1.28 | 1.65 | 1.83 | 1.95 | 1.38 | 1.61 | 1.88^{b}, 1.72^{f} |

↓-HOMO | ^{3}Δ | 1.30^{h} | 0.31 | 1.00 | 1.29 | 1.43 | 1.55 | 1.06 | 1.21 | 1.19^{b}, 1.18^{c}, 1.14^{i} |

MAE | 1.31 | 0.33 | 0.60 | 1.01 | 1.37 | 0.25 | 0.32 | |||

CuO^{−} (^{1}Σ^{+}) | ||||||||||

HOMO | ^{2}Π | 1.78^{j} | 0.40 | 1.40 | 1.58 | 1.40 | 0.97 | 2.19 | 3.17 | 1.55^{b}, 1.52^{c}, 0.46^{k} |

HOMO-1 | ^{2}Σ^{+} | 2.75^{j} | 1.39 | 2.17 | 2.23 | 1.98 | 1.56 | 2.66 | 3.58 | 2.96^{b}, 2.86^{e}, 1.60^{k}, 2.78^{l}, 2.47^{m}, 2.81^{n} |

HOMO-2 | 4.50^{o} | 2.18 | 4.05 | 4.60 | 4.56 | 4.25 | 4.88 | 6.38 | 4.07^{b}, 4.01^{e}, 4.58^{l}, 4.50^{m} | |

HOMO-3 | 2.89 | 4.49 | 5.06 | 5.04 | 4.67 | 4.96 | 6.60 | |||

HOMO-4 | 3.70 | 4.49 | 4.63 | 4.53 | 4.24 | 5.16 | 6.39 | |||

MAE | 1.69 | 0.47 | 0.27 | 0.40 | 0.75 | 0.29 | 1.37 | |||

ZnO^{−} (^{2}Σ^{+}) | ||||||||||

↑-HOMO | ^{1}Σ^{+} | 2.09^{p} | 0.91 | 1.91 | 2.20 | 2.23 | 2.00 | 2.11 | 2.57 | 2.19^{b}, 2.33^{c}, 2.29^{q}, 2.10^{r}, 1.06^{k} |

↑-HOMO-1 | ^{1}Π | 2.71^{p} | 1.36 | 2.43 | 3.03 | 3.72 | 4.79 | 2.77 | 3.73 | 2.62^{b}, 1.43^{k} |

↑-HOMO-2 | 3.24 | 4.04 | 4.89 | 5.64 | 6.90 | 4.77 | 5.66 | 3.50^{k} | ||

↓-HOMO | ^{3}Π | 2.40^{p} | 1.17 | 2.17 | 2.64 | 3.15 | 4.02 | 2.65 | 3.48 | 2.41^{b}, 1.20^{k} |

↓-HOMO-1 | ^{3}Σ^{+} | 3.96^{p} | 2.71 | 3.35 | 3.47 | 3.31 | 3.00 | 4.11 | 4.51 | 4.15^{b}, 2.89^{k} |

MAE | 1.25 | 0.33 | 0.29 | 0.64 | 1.19 | 0.12 | 0.78 |

. | . | . | G_{0}W_{0}@PBEα
. | G_{n}W_{0}@
. | G_{n}W_{n}@
. | . | ||||
---|---|---|---|---|---|---|---|---|---|---|

. | State . | Exp. . | α = 0.00
. | α = 0.25
. | α = 0.50
. | α = 0.75
. | α = 1.00
. | PBE . | PBE . | Others . |

ScO^{−} (^{1}Σ^{+}) | ||||||||||

HOMO | ^{1}Σ^{+} | 1.35^{a} | 0.51 | 1.15 | 1.45 | 1.59 | 1.63 | 1.21 | 1.35 | 1.28^{b}, 1.26^{c}, 1.19^{d} |

HOMO-1 | ^{2}Δ | 3.10^{a} | 3.30 | 4.77 | 5.45 | 5.72 | 5.82 | 5.34 | 6.08 | 2.41^{b}, 2.78^{e}, 3.31^{f} |

HOMO-2 | ^{2}Π | 3.40^{a} | 3.42 | 4.81 | 5.40 | 5.61 | 5.63 | 5.39 | 6.17 | 3.34^{b}, 3.24^{e}, 3.44^{f} |

MAE | 0.84^{g} | 0.18^{g} | 0.10^{g} | 0.24^{g} | 0.28^{g} | 0.14^{g} | 0.00^{g} | |||

TiO^{−} (^{2}Δ) | ||||||||||

↑-HOMO | ^{1}Σ^{+} | 2.00^{h} | 0.26 | 2.24 | 3.70 | 4.79 | 5.63 | 1.83 | 2.74 | 2.39^{b}, 2.37^{e}, 2.34^{f} |

↑-HOMO-1 | ^{1}Δ | 1.73^{h} | 0.53 | 1.28 | 1.65 | 1.83 | 1.95 | 1.38 | 1.61 | 1.88^{b}, 1.72^{f} |

↓-HOMO | ^{3}Δ | 1.30^{h} | 0.31 | 1.00 | 1.29 | 1.43 | 1.55 | 1.06 | 1.21 | 1.19^{b}, 1.18^{c}, 1.14^{i} |

MAE | 1.31 | 0.33 | 0.60 | 1.01 | 1.37 | 0.25 | 0.32 | |||

CuO^{−} (^{1}Σ^{+}) | ||||||||||

HOMO | ^{2}Π | 1.78^{j} | 0.40 | 1.40 | 1.58 | 1.40 | 0.97 | 2.19 | 3.17 | 1.55^{b}, 1.52^{c}, 0.46^{k} |

HOMO-1 | ^{2}Σ^{+} | 2.75^{j} | 1.39 | 2.17 | 2.23 | 1.98 | 1.56 | 2.66 | 3.58 | 2.96^{b}, 2.86^{e}, 1.60^{k}, 2.78^{l}, 2.47^{m}, 2.81^{n} |

HOMO-2 | 4.50^{o} | 2.18 | 4.05 | 4.60 | 4.56 | 4.25 | 4.88 | 6.38 | 4.07^{b}, 4.01^{e}, 4.58^{l}, 4.50^{m} | |

HOMO-3 | 2.89 | 4.49 | 5.06 | 5.04 | 4.67 | 4.96 | 6.60 | |||

HOMO-4 | 3.70 | 4.49 | 4.63 | 4.53 | 4.24 | 5.16 | 6.39 | |||

MAE | 1.69 | 0.47 | 0.27 | 0.40 | 0.75 | 0.29 | 1.37 | |||

ZnO^{−} (^{2}Σ^{+}) | ||||||||||

↑-HOMO | ^{1}Σ^{+} | 2.09^{p} | 0.91 | 1.91 | 2.20 | 2.23 | 2.00 | 2.11 | 2.57 | 2.19^{b}, 2.33^{c}, 2.29^{q}, 2.10^{r}, 1.06^{k} |

↑-HOMO-1 | ^{1}Π | 2.71^{p} | 1.36 | 2.43 | 3.03 | 3.72 | 4.79 | 2.77 | 3.73 | 2.62^{b}, 1.43^{k} |

↑-HOMO-2 | 3.24 | 4.04 | 4.89 | 5.64 | 6.90 | 4.77 | 5.66 | 3.50^{k} | ||

↓-HOMO | ^{3}Π | 2.40^{p} | 1.17 | 2.17 | 2.64 | 3.15 | 4.02 | 2.65 | 3.48 | 2.41^{b}, 1.20^{k} |

↓-HOMO-1 | ^{3}Σ^{+} | 3.96^{p} | 2.71 | 3.35 | 3.47 | 3.31 | 3.00 | 4.11 | 4.51 | 4.15^{b}, 2.89^{k} |

MAE | 1.25 | 0.33 | 0.29 | 0.64 | 1.19 | 0.12 | 0.78 |

^{a}

Reference 24.

^{b}

Reference 61 using 6–3111 + G^{*} basis sets.

^{c}

Reference 10 using the B3LYP functional.

^{d}

Reference 11 using the B3LYP functional.

^{e}

Reference 62.

^{f}

Reference 63 using the multireference configuration interaction (MRCI) method.

^{g}

HOMO-1 and HOMO-2 are not included because our *GW* calculations cannot account for two-electron transitions (see text).

^{h}

Reference 25.

^{i}

Reference 65 using the B3LYP functional.

^{j}

Reference 26.

^{k}

Reference 54 using the *G*_{0}*W*_{0}@PBE method.

^{l}

Reference 66 using the CCSD(T) method.

^{m}

Reference 67.

^{n}

Reference 68 using the single and double excitation configuration interaction (SDCI) method.

^{o}

We chose this value from the Z band in the PES spectrum of CuO^{−} (see text).

^{p}

Reference 27.

^{q}

Reference 64 using the B3LYP functional.

^{r}

Reference 21 using the *G*_{0}*W*_{0}@PBE0 method.

Throughout this work, we use only gKS-PBE TM 3*d* character except when we discuss the subtle competition between direct and indirect relativistic effects because the EXX amount has a small effect on TM 3*d* character. Also, throughout this work, we use only the gKS-PBE orbital order to avoid confusion. The orbital order depends strongly on the amount of EXX (e.g., PBE vs HF) and the level of theory (e.g., DFT vs *GW*).^{8,16} For example, gKS-PBE HOMO of CuO^{−} corresponds to gKS-PBE*α*(*α* = 1.00) HOMO-1 and *G*_{0}*W*_{0}@PBE*α*(*α* = 1.00) HOMO of CuO^{−}, as shown in the supplementary material and Fig. 6. In this work, gKS-PBE and *G*_{0}*W*_{0}@PBE were found to have the same orbital order.

### A. *G*_{0}*W*_{0} starting points

Figures 5 and 6 show PES and *G*_{0}*W*_{0}@PBE*α*(0.00 ≤ *α* ≤ 1.00) QP spectra of ScO^{−}, TiO^{−}, CuO^{−}, and ZnO^{−}. In PES spectra, vertical dashed and solid lines represent experimental *sp*- and *d*-electron binding energies, respectively. In *GW* spectra, oblique dashed and solid lines track calculated *sp*- and *d*-electron binding energies, respectively. In Figs. 5 and 6, we find a few general trends common in all TMO anions considered in this work: (i) no *G*_{0}*W*_{0}@PBE*α*(0.00 ≤ *α* ≤ 1.00) results are in perfect agreement with experiment, (ii) *G*_{0}*W*_{0}@PBE underestimates the IE of TMO anions by ∼1 eV, which is larger than the typical underestimation for *sp* molecules (∼0.5 eV),^{13,15,16} and (iii) *G*_{0}*W*_{0}@PBE*α*(0.25 ≤ *α* ≤ 0.50) reduces it to ∼0.1 eV. In the following, we analyze each TMO anion individually.

#### 1. ScO^{−}

Scandium is the first transition metal and has only one 3*d* electron. DFT and CCSD(T) calculations in Refs. 10 and 11 confirmed the ground state of ScO^{−} as ^{1}Σ^{+} (8*σ*^{2}3*π*^{4}9*σ*^{2}), correcting the wrongly assumed state ^{3}Δ^{−} (8*σ*^{2}3*π*^{4}9*σ*^{1}1*δ*) in Ref. 24. There is no 3*d* peak or band in the PES spectrum of ScO^{−}, and the top three valence 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*_{0}*W*_{0}@PBE binding energies slightly overestimate PES ones (by 0.20 and 0.02 eV for HOMO-1 and HOMO-2, respectively), whereas *G*_{0}*W*_{0}@PBE*α*(0.25 ≤ *α* ≤ 1.00) binding energies significantly overestimate PES ones by ∼2 eV [e.g., for HOMO-1, *G*_{0}*W*_{0}@PBE*α*(*α* = 0.25) and *G*_{0}*W*_{0}@PBE*α*(*α* = 1.00) binding energies are greater than PES ones by 1.67 and 2.72 eV, respectively]. This seems to suggest that for HOMO-1 and HOMO-2 of ScO^{−}, *G*_{0}*W*_{0}@PBE performs better than *G*_{0}*W*_{0}@PBE*α*(0.25 ≤ *α* ≤ 1.00), but this is not the case due to the nature of the corresponding peaks in the PES experiment. Reference 11 suggests that the second and third peaks in the PES spectrum of ScO^{−} are likely due to two-electron transitions from 8*σ*^{2}3*π*^{4}9*σ*^{2} (^{1}Σ^{+} ScO^{−}) to 8*σ*^{2}3*π*^{4}10*σ* (*B*^{2}Σ^{+} ScO) and to 8*σ*^{2}3*π*^{4}1*δ* (*A*′^{2}Δ ScO) states, respectively, which *GW* calculations for quasiparticle excitations cannot account for. In other words, the seemingly excellent agreement between *G*_{0}*W*_{0}@PBE and PES binding energies for HOMO-1 and HOMO-2 of ScO^{−} is accidental. Therefore, we exclude HOMO-1 and HOMO-2 of ScO^{−} from our evaluation of the performance of *GW* schemes in the following.

We also see that as *α* increases, *G*_{0}*W*_{0}@PBE*α* IE always increases, but this happens at different rates: As *α* increases from 0.00 to 0.25, it increases rapidly, whereas as *α* increases from 0.25 to 1.00, it increases slowly. The weak sensitivity of *G*_{0}*W*_{0}@PBE*α*(0.25 ≤ *α* ≤ 1.00) IE to a change in *α* gives a large margin for an optimal amount of EXX: 25%–100%.

#### 2. TiO^{−}

Titanium is the second transition metal and has two 3*d* electrons. Several theoretical studies in Refs. 10, 61, 63, 65, and 69 confirmed 9*σ*^{2}*δ*^{1} (^{2}Δ) as the ground-state electron configuration of TiO^{−}, correcting the wrongly assigned configuration 9*σ*^{1}*δ*^{2} (^{1}Σ^{+}) in Ref. 25. Unlike ScO^{−}, which has an empty *δ* shell, TiO^{−} has one 3*d* electron in the *δ* shell. The transition of the 3*d* electron from 9*σ*^{2}*δ*^{1} (^{2}Δ TiO^{−}) to 9*σ*^{2} (^{1}Σ^{+} TiO) states produces the third peak in the PES spectrum of TiO^{−} at 2.0 eV. In the *G*_{0}*W*_{0}@PBE*α* QP spectrum of TiO^{−}, ↑-HOMO is of entirely Ti 3*d* character, as shown in Table II.

The right panel of Fig. 5 clearly shows that *G*_{0}*W*_{0}@PBE*α*(0.00 ≤ *α* ≤ 1.00) binding energy for ↑-HOMO of TiO^{−} is much more sensitive to a change in *α* than those for other occupied molecular orbitals with mainly *sp* character, as shown in Table II. The orbital-character-dependent sensitivity of *G*_{0}*W*_{0}@PBE*α* binding energy to a change in *α* causes a couple of problems. First, *G*_{0}*W*_{0}@PBE underestimates the IE and the 3*d*-electron binding energy of TiO^{−} nonuniformly (by 0.99 and 1.74 eV, respectively), leading to the wrong orbital order. In other words, *G*_{0}*W*_{0} does not correct the wrong orbital order produced by PBE. Second, the *G*_{0}*W*_{0} starting-point approach does not give accurate results for both the IE and the 3*d*-electron binding energy of TiO^{−} at the same time. For example, *G*_{0}*W*_{0}@PBE*α*(*α* = 0.50) gives a better result for the IE of TiO^{−} by 0.29 eV, but a worse result for the 3*d*-electron binding energy of TiO^{−} by 1.46 eV, than *G*_{0}*W*_{0}@PBE*α*(*α* = 0.25). This type of behavior is not uncommon in *GW* predictions for transition metal oxides; for example, no existing *GW* scheme can accurately reproduce both the bandgap and the *d*-band position in the band structure of bulk ZnO at the same time.^{17–19}

The increase in *α* from 0 to 1 has a similar effect on *G*_{0}*W*_{0}@PBE*α* IE of both ScO^{−} and TiO^{−}: For both ScO^{−} and TiO^{−}, *G*_{0}*W*_{0}@PBE*α*(0.25 ≤ *α* ≤ 1.00) reduces the underestimation of IE by *G*_{0}*W*_{0}@PBE from ∼1 eV to ∼0.1 eV [e.g., *G*_{0}*W*_{0}@PBE*α*(*α* = 0.25) reduces the difference in IE between PES and *G*_{0}*W*_{0}@PBE from 0.84 eV to 0.20 eV and from 0.99 eV to 0.30 eV, respectively]. However, unlike ScO^{−}, the strong sensitivity of *G*_{0}*W*_{0}@PBE*α* 3*d*-electron binding energy of TiO^{−} to a change in *α* gives a small margin for an optimal amount of EXX: ∼25%.

#### 3. CuO^{−}

Copper is the 11th transition metal and has ten 3*d* electrons. DFT calculations in Ref. 10 confirmed the ground state of CuO^{−} as ^{1}Σ^{+} with the electron configuration of 3*d*^{10}2*pσ*^{2}2*pπ*^{4}. There are three bands (named as X, Y, and Z in Ref. 26) in the PES spectrum of CuO^{−}, as shown in the top of the left panel of Fig. 6. Reference 26 suggested that the photodetachment transition of 3*d* electrons (3*dδ*^{4}3*dπ*^{4}3*dσ*^{2}) from ^{1}Σ^{+} CuO^{−} 3*d*^{10}2*pσ*^{2}2*pπ*^{4} to Z CuO 3*d*^{9}2*pσ*^{2}2*pπ*^{4} states produces the broad Z band in the PES spectrum of CuO^{−} at ∼4.5 eV (which we selected from the position of the highest peak in the Z band) and assumed that the Z band is unusually broad likely due to a large geometry change from the anion to the neutral. In the *G*_{0}*W*_{0}@PBE*α* QP spectrum of CuO^{−}, HOMO-2 is of entirely Cu 3*d* character, as shown in Table II.

In the left panel of Fig. 6, we see that *G*_{0}*W*_{0}@PBE*α*(0.00 ≤ *α* ≤ 0.50) binding energy for HOMO-2 of CuO^{−} is more sensitive to a change in *α* than those for other occupied molecular orbitals with weaker Cu 3*d* character than HOMO-2, as shown in Table II, and *G*_{0}*W*_{0}@PBE*α*(0.50 ≤ *α* ≤ 0.75) gives good results for the IE and the 3*d*-electron binding energy (corresponding to HOMO and HOMO-2, respectively) of CuO^{−} at the same time. Scalar relativistic effects in ECP reduce *G*_{0}*W*_{0}@PBE*α*(*α* = 0.50) and *G*_{0}*W*_{0}@PBE*α*(*α* = 0.75) binding energies for HOMO-2 of CuO^{−} by 0.31 and 0.24 eV, as shown in the supplementary material, without changing the conclusion that *G*_{0}*W*_{0}@PBE*α*(0.50 ≤ *α* ≤ 0.75) gives good results for the 3*d*-electron binding energy of CuO^{−}. Like TiO^{−}, the orbital-character-dependent sensitivity of *G*_{0}*W*_{0}@PBE*α* binding energy to a change in *α* causes *G*_{0}*W*_{0}@PBE to underestimate the IE and the 3*d*-electron binding energy of CuO^{−} nonuniformly (by 1.38 and 2.32 eV, respectively). *G*_{0}*W*_{0}@PBE*α*(0.50 ≤ *α* ≤ 1.00) binding energies for all valence molecular orbitals considered in this work are weakly sensitive to a change in *α*. This trend suggests that PBE*α*(0.50 ≤ *α* ≤ 1.00) orbitals with large amounts of EXX are good for localized states of CuO^{−} with strong 3*d* character [i.e., for CuO^{−}, PBE*α*(0.50 ≤ *α* ≤ 1.00) wavefunctions are close to QP ones] and is consistent with the relatively good performance of HF on molecules with weak screening.^{23,50}

#### 4. ZnO^{−}

Zinc is the 12th transition metal and has ten 3*d* electrons. Zinc is rather distinct from other first row transition metals due to its closed-shell electron configuration. In other words, zinc is more similar to alkaline earth metals than other transition metals because Zn 3*d* electrons generally do not participate in bonding.^{70} DFT calculations in Ref. 10 confirmed the ground-state electron configuration of ZnO^{−} as ^{2}Σ^{+} 10*σ*^{1}9*σ*^{2}4*π*^{4}*δ*^{4}. There are four bands 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 to insufficient photon energy of 4.66 eV.^{27} Unlike CuO^{−}, all valence molecular orbitals in the *G*_{0}*W*_{0}@PBE*α* QP spectrum of ZnO^{−} have weak Zn 3*d* character, as shown in Table II.

The right panel of Fig. 6 shows that unlike CuO^{−}, *G*_{0}*W*_{0}@PBE underestimates electron binding energies for 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) possibly because all valence molecular orbitals have similar Zn 3*d* character and thus their *G*_{0}*W*_{0}@PBE*α* binding energies have similar sensitivity to a change in *α*. *G*_{0}*W*_{0}@PBE*α*(*α* = 0.50) gives good results for the IE and the orbital order of ZnO^{−} at the same time. Unlike CuO^{−}, *G*_{0}*W*_{0}@PBE*α*(0.50 ≤ *α* ≤ 1.00) binding energies for all valence molecular orbitals, except for ↑-HOMO and ↓-HOMO-1, are strongly sensitive to a change in *α*. This trend can be explained in terms of the importance of spin-splitting in open-shell molecules (as has been discussed in Ref. 53 for the Cu$O2\u2212$ molecule).

### B. ev*GW* self-consistency levels

Figure 7 shows PES and ev*GW* QP spectra of ScO^{−}, TiO^{−}, CuO^{−}, and ZnO^{−}. In the following, we analyze *G*_{n}*W*_{0}@PBE and *G*_{n}*W*_{n}@PBE results individually.

In Fig. 7, we see that as the ev*GW* self-consistency level increases from *G*_{0}*W*_{0} to *G*_{n}*W*_{n}, *GW* binding energies always increase, but this occurs at different rates: As the ev*GW* self-consistency level increases from *G*_{0}*W*_{0} to *G*_{n}*W*_{0}, *GW* binding energies increase rapidly (e.g., the IE increases by 0.70, 0.75, 1.79, and 1.20 eV for ScO^{−}, TiO^{−}, CuO^{−}, and ZnO^{−}, respectively), while as it increases from *G*_{n}*W*_{0} to *G*_{n}*W*_{n}, they increase slowly (e.g., the IE increases by 0.14, 0.15, and 0.46 eV for ScO^{−}, TiO^{−}, and ZnO^{−}, respectively) except for CuO^{−} (0.98 eV), which will be discussed later. *G*_{0}*W*_{0}@PBE always underestimates electron binding energies, whereas *G*_{n}*W*_{n}@PBE generally overestimates them. *G*_{n}*W*_{0}@PBE binding energies are always in between *G*_{0}*W*_{0}@PBE and *G*_{n}*W*_{n}@PBE ones and generally close to experiment. In other words, *G*_{0}*W*_{0}@PBE and *G*_{n}*W*_{n}@PBE act as lower and upper bounds for *G*_{n}*W*_{0}@PBE, generally producing overscreenings and underscreenings, respectively. This trend of the ev*GW* self-consistency approach in electronic structure of molecules is also observed in the band structure of solids.^{17}

We also see that the ev*GW* self-consistency has a strong effect on *GW* binding energies for molecular orbitals with strong 3*d* character (e.g., ↑-HOMO of TiO^{−} and HOMO-2 of CuO^{−}). For example, *G*_{n}*W*_{0}@PBE reduces the underestimation errors of *G*_{0}*W*_{0}@PBE in the IE and the 3*d*-electron binding energy of TiO^{−} with respect to experiment from 0.99 to 1.74 eV to 0.24 and 0.17 eV, respectively. As a result, *G*_{n}*W*_{0}@PBE corrects the wrong *G*_{0}*W*_{0}@PBE orbital order in TiO^{−}. Another example is that *G*_{n}*W*_{0}@PBE gives small (∼0.1 eV) errors in electron binding energies for all valence molecular orbitals of ZnO^{−}, which are uniformly underestimated by *G*_{0}*W*_{0}@PBE by ∼1 eV due to similarly weak Zn 3*d* character. For ZnO^{−}, *G*_{0}*W*_{0}@PBE and *G*_{n}*W*_{0}@PBE yield mean absolute errors (MAEs) of 1.25 and 0.12 eV, respectively, as shown in Table III.

CuO^{−} exhibits particularly large differences between *G*_{n}*W*_{0}@PBE and *G*_{n}*W*_{n}@PBE binding energies compared to other TMO anions. This trend is not associated with scalar relativistic effects in ECP, which reduce *G*_{n}*W*_{0}@PBE and *G*_{n}*W*_{n}@PBE binding energies by similar amounts (e.g., by 0.57 and 0.66 eV, respectively, for HOMO-2 of CuO^{−}, as shown in the supplementary material). We attribute this trend to strong 3*d* character in molecular orbitals of CuO^{−}. For example, CuO^{−} has a larger difference between *G*_{n}*W*_{0}@PBE and *G*_{n}*W*_{n}@PBE IEs than ScO^{−} (0.98 and 0.14 eV, respectively) possibly because CuO^{−} has stronger 3*d* character in HOMO than ScO^{−} (23% and 6%, respectively, as shown in Table II).

### C. Comparison of *G*_{0}*W*_{0} starting-point and ev*GW* self-consistency approaches

From our results presented so far, it appears that both *G*_{0}*W*_{0} starting-point and ev*GW* self-consistent approaches can, in principle, be good *GW* methods for finite systems: both *G*_{0}*W*_{0}@PBE*α*(0.25 ≤ *α* ≤ 0.50) and *G*_{n}*W*_{0}@PBE can reduce the large and orbital-character-dependent nonuniform errors for electron binding energies of TMO anions produced by *G*_{0}*W*_{0}@PBE with respect to experiment from ∼1–2 eV to ∼0.1–0.5 eV. Reference 18 obtained similar results for extended systems: both *G*_{0}*W*_{0}@PBE*α*(*α* = 0.25) and *G*_{n}*W*_{0}@PBE give satisfactory results for the bandgap and the *d*-electron binding energy of solids and drew the conclusions that (i) for accuracy, one can choose either *G*_{0}*W*_{0}@PBE*α*(*α* = 0.25) or *G*_{n}*W*_{0}@PBE because they give similar results, but (ii) for efficiency, one may want to choose *G*_{0}*W*_{0}@PBE*α*(*α* = 0.25) over *G*_{n}*W*_{0}@PBE because the former is computationally cheaper than the latter. However, in the case of molecular systems, we argue that *G*_{n}*W*_{0}@PBE has several practical advantages over *G*_{0}*W*_{0}@PBE*α*.

First, *G*_{n}*W*_{0}@PBE does not contain system-dependent adjustable parameters. Unlike extended systems, there is no unique amount of EXX for the *G*_{0}*W*_{0} starting point, which works well for all finite systems. For example, we show in Sec. IV A that 25% EXX is optimal for ScO^{−} and TiO^{−}, whereas 50% EXX is optimal for CuO^{−} and ZnO^{−}. Also, it appears that atoms and small molecules require more amount of EXX than clusters and large molecules.^{41} Second, *G*_{n}*W*_{0}@PBE is transferable between finite and extended systems. *G*_{n}*W*_{0}@PBE works well for both molecules and solids (e.g., ZnO^{−} anion and bulk ZnO, respectively).^{17} This greatly extends the range of applicability of the *GW* method. For example, *G*_{n}*W*_{0}@PBE may be applicable to solid-molecule hybrid systems such as molecular junctions and molecules adsorbed on solid surfaces.^{71} Also, *G*_{n}*W*_{0}*@*PBE may be used for the study of quantum size effects in clusters because it is independent of the cluster size. Third, *G*_{n}*W*_{0}*@*PBE is easy to use and reliable. Unlike PBE*α*(0.00 < *α* ≤ 1.00), PBE is safe from the SCF convergence issue, and unlike *G*_{0}*W*_{0}, ev*GW* with *Z* = 1 is immune to the *GW* multisolution issue. Therefore, *G*_{n}*W*_{0}*@*PBE does not need manual, time-consuming, and error-prone tests to address the two issues, which are explained in detail in Sec. III.

Furthermore, *G*_{n}*W*_{0}*@*PBE has a few desirable properties. One of them comes from the PBE part. PBE causes the smallest incomplete basis set error, as shown in Fig. 2, allowing one to use smaller basis sets for the CBS limit, which makes *G*_{n}*W*_{0}*@*PBE cheaper. Two desirable properties come from the *GW* part. *G*_{n}*W*_{0} (as well as *G*_{n}*W*_{n}) gives faster and more stable *GW* convergence than QS*GW* and depends more weakly on the choice of *η* (e.g., we used a single value of *η* for ev*GW* in this work), as discussed in Sec. III C 4. Also, *G*_{n}*W*_{0} is cheaper than *G*_{n}*W*_{n}, as pointed out in Ref. 17 and discussed in Sec. II G. In fact, *G*_{n}*W*_{0} is the cheapest self-consistent *GW* scheme.

One may argue that *G*_{0}*W*_{0}*@*PBE*α* should be a choice of *GW* methods because it is computationally more efficient than *G*_{n}*W*_{0}*@*PBE by the number of self-consistent *G*_{n}*W*_{0} iterations. However, as discussed in Sec. II G and the supplementary material, this is not the case since the compute time difference between *G*_{0}*W*_{0}*@*PBE*α* and *G*_{n}*W*_{0}*@*PBE does not depend only on the number of *G*_{n}*W*_{0} iterations; there are other factors such as the number of eigenvalues to update for $\u03f5mGnW0$, the number of frequency points to use for Σ_{c}(*ω*), the number of Δ*ω* and *η* values to test for $\u03f5mG0W0$, and the number of initial guess wavefunctions to test for gKS calculations. Some factors can cancel each other out; for example, *G*_{n}*W*_{0}*@*PBE requires a few *G*_{n}*W*_{0} iterations, but one typically needs to test a few *η* values for *G*_{0}*W*_{0}*@*PBE*α*. In other words, when all factors are taken into account, the total compute time to obtain reliable and reproducible QP spectra at *G*_{0}*W*_{0}*@*PBE*α* and *G*_{n}*W*_{0}*@*PBE levels of theory can be comparable, as is especially the case for open-shell systems.

### D. Comparison with results in the literature

Some of our results for the performance of *G*_{0}*W*_{0} starting-point and ev*GW* self-consistency approaches in this work may seem to be at odds with some of the results in the literature. In this section, we discuss the origin of the apparent differences between them.

We begin with the *G*_{0}*W*_{0} starting-point approach. Table IV summarizes a few selected results for the optimal amount of EXX in the *G*_{0}*W*_{0} starting point out of numerous results, such as Refs. 14 and 72, in the literature. Interestingly, we see that there is a wide range of EXX amounts from 25% to 100%, and Refs. 4 and 23 obtained different results (50% and 100%, respectively) from the same set of molecules. It seems that 75% and 100% are too large compared to our results: 25%–50%. One may guess that the large difference is due to implementation differences such as basis type (e.g., Gaussian vs PW) and frequency integration type (e.g., analytical vs numerical). However, Refs. 15 and 16 showed that such implementation differences have little effect on *G*_{0}*W*_{0} IE (∼0.06 eV). There are a couple of other factors that have a stronger effect on *G*_{0}*W*_{0} results than implementation differences. One factor is the choice of system and property. As shown in Sec. IV A, *G*_{0}*W*_{0}*@*PBE*α*(0.25 ≤ *α* ≤ 1.00) IEs of *sp* systems are slightly different (by ∼0.1 eV). Most existing *G*_{0}*W*_{0} studies used the IE of *sp*-bonded systems to determine the optimal amount of EXX in the *G*_{0}*W*_{0} starting point. The other factor is that the choice of QP equation solver and CBS extrapolation method. As shown in Secs. III C 1 and III C 3, the linearization method and the CBS extrapolation method (e.g., whether to extrapolate or not and which fitting function and basis set to use for extrapolation) can cause a difference in *G*_{0}*W*_{0} IE on the order of ∼0.1 eV. Overall, the combination of the two factors gives a large margin for the optimal amount of EXX in the *G*_{0}*W*_{0} starting point and thus is likely to produce the wide range of amounts that exist in the literature.

Reference . | Körbel et al.^{21}
. | Bruneval et al.^{4}
. | Kaplan et al.^{22}
. | Rostgaard et al.^{23}
. | This work . |
---|---|---|---|---|---|

Code | FIESTA | MOLGW | TURBOMOLE | GPAW | MOLGW |

Optimal EXX (%) | 25 | 50 | 75 | 100 | 25–50 |

Tested EXX (%) | 25 and 100 | 0, 20, 25 and 50 | 0, 25 and 75 | 0 and 100 | 0, 25, 50, 75 and 100 |

System | 39 closed-shell 3,4, | 34 closed-shell sp^{a} | 29 closed-shell sp | 34 closed-shell sp^{a} | 4 closed- and open- |

5d and 9 closed-shell | shell 3d | ||||

sp | |||||

System size | 2–7 atoms | 2–8 atoms | 2–18 atoms | 2–8 atoms | 2 atoms |

Property | HOMO and LUMO | HOMO | HOMO and HOMO-n | HOMO | HOMO-n (n = 0, 1, …) |

(n = 0, 1, …)^{b} | (focusing on 3d MO) | ||||

Reference data | Experiment | ΔSCF^{c} | QSGW^{d} | Experiment | Experiment |

ω integration | Contour | Fully analytic | Fully analytic | Fully analytic | Fully analytic |

QP equation | deformation^{e} | Linearization | Spectral function | Linearization^{f} and | Graphical solution and |

Linearization | Spectral function^{g} | Spectral function | |||

η | Not available | Not available | 0.001 eV | 0 eV | 0.002 or 0.005 hartree |

CBS 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} | AE |

RI | Used | Not used | Used | Not applicable^{i} | Not used |

Reference . | Körbel et al.^{21}
. | Bruneval et al.^{4}
. | Kaplan et al.^{22}
. | Rostgaard et al.^{23}
. | This work . |
---|---|---|---|---|---|

Code | FIESTA | MOLGW | TURBOMOLE | GPAW | MOLGW |

Optimal EXX (%) | 25 | 50 | 75 | 100 | 25–50 |

Tested EXX (%) | 25 and 100 | 0, 20, 25 and 50 | 0, 25 and 75 | 0 and 100 | 0, 25, 50, 75 and 100 |

System | 39 closed-shell 3,4, | 34 closed-shell sp^{a} | 29 closed-shell sp | 34 closed-shell sp^{a} | 4 closed- and open- |

5d and 9 closed-shell | shell 3d | ||||

sp | |||||

System size | 2–7 atoms | 2–8 atoms | 2–18 atoms | 2–8 atoms | 2 atoms |

Property | HOMO and LUMO | HOMO | HOMO and HOMO-n | HOMO | HOMO-n (n = 0, 1, …) |

(n = 0, 1, …)^{b} | (focusing on 3d MO) | ||||

Reference data | Experiment | ΔSCF^{c} | QSGW^{d} | Experiment | Experiment |

ω integration | Contour | Fully analytic | Fully analytic | Fully analytic | Fully analytic |

QP equation | deformation^{e} | Linearization | Spectral function | Linearization^{f} and | Graphical solution and |

Linearization | Spectral function^{g} | Spectral function | |||

η | Not available | Not available | 0.001 eV | 0 eV | 0.002 or 0.005 hartree |

CBS 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} | AE |

RI | Used | Not used | Used | Not applicable^{i} | Not used |

^{a}

The same set of molecules is used.

^{b}

For naphthalene only.

^{c}

Reference 4 showed that ΔSCF using CCSD(T) with CN = 4 causes an error of ∼0.1 eV in the IE of small *sp* molecules with respect to experiment (the largest being 0.67 eV for NaCl).

^{d}

Reference 57 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 to experiment (the largest being ∼0.4 eV for O).

^{e}

References 13 and 34 showed that the contour deformation technique produces almost the same *GW* self-energy as the fully analytic method for frontier and nonfrontier 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.

Next, we move on to the ev*GW* self-consistency approach and discuss the origin of apparently conflicting ev*GW* results for IE and starting-point dependency. First, Ref. 7 reported that the *G*_{n}*W*_{n} approach with a local-density approximation starting point gives good results for the IE of large *sp* molecules, whereas we found in Sec. IV B that *G*_{n}*W*_{0}*@*PBE gives satisfactory results for the electronic structure (including the IE) of small 3*d* molecules. A comparison of ev*GW* implementations in Ref. 7 and this work is provided in the supplementary material. We believe that the main origin of the different results is the orbital-character-dependent sensitivity of ev*GW* binding energy to a change in ev*GW* self-consistency level. As shown in Sec. IV B, *G*_{n}*W*_{0}*@*PBE and *G*_{n}*W*_{n}*@*PBE binding energies are slightly different for delocalized HOMO with weak 3*d* character by ∼0.1 eV but significantly different for localized HOMO with strong 3*d* character by ∼1 eV. Unlike this work, Ref. 7 used the linearization method, employed pseudopotentials and RI, and did not use the CBS limit, but these cause small (∼0.1 eV) differences in ev*GW* IE, as shown in Sec. III C. Accordingly, they are most likely not the reason for the large (0.98 eV) difference between *G*_{n}*W*_{0}*@*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 found in Sec. III C 4 that in TMO anions, *G*_{n}*W*_{n} sometimes depends more strongly on the starting point than *G*_{n}*W*_{0}. As mentioned in Sec. III C 4, we believe that the orbital character influences the ev*GW* starting-point dependency: for molecular orbitals with strong (weak) 3*d* character, *G*_{n}*W*_{n} depends more strongly (weakly) on the starting point than *G*_{n}*W*_{0}. Overall, without molecular orbitals with strong 3*d* character (e.g., HOMO of CuO^{−}), our ev*GW* results for IE and starting-point dependency in this work are consistent with those in Refs. 7 and 41.

To verify our idea about the origin of the seemingly different results between this work and the literature, we performed a simple test: (i) we chose ScO^{−} and TiO^{−} as our analogs of *sp* molecules in the literature because their valence molecular orbitals have weak transition-metal character, except for *↑*-HOMO of TiO^{−} with entirely Ti 3*d* character; (ii) we applied 15 different starting-point–self-consistency hybrid *GW* schemes (*G*_{0}*W*_{0}, *G*_{n}*W*_{0}, and *G*_{n}*W*_{n}; 0%, 25%, 50%, 75%, and 100% EXX) to them; and (iii) we searched for *GW* schemes that give a reasonably small error of less than 0.5 eV in the IE and the 3*d*-electron binding energy with respect to experiment. Figure 8 shows the results of the test. We see that *GW* IEs of ScO^{−} and TiO^{−} (red dashed lines) depend weakly on the starting point and the self-consistency level, giving a large margin for the choice of *GW* schemes. 14 *GW* schemes out of 15 (*G*_{0}*W*_{0}*@*PBE is an exception as expected) give a small error (less than 0.5 eV), which explains 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* 3*d*-electron binding energy of TiO^{−} (green solid lines) depends strongly on the starting point and the self-consistency level, yielding a small margin for the choice of *GW* schemes. Only two *GW* schemes (*G*_{0}*W*_{0}*@*PBE0 and *G*_{n}*W*_{0}*@*PBE) out of 15 give a small error (less than 0.5 eV), which is why we obtained a small number of good *GW* schemes for the electronic structure of *d* molecules in this work. Overall, we confirm that evaluation results for the performance of *GW* schemes depend strongly on the choice of system and property (e.g., the IE with mainly *sp* character vs the electronic structure containing *d* states).

## V. SUMMARY AND CONCLUSIONS

In summary, we calculated the electronic structure of closed- and open-shell molecular anions with partially and completely filled 3*d* shells (shallow and deep 3*d* states, respectively) using various *GW* schemes and compared calculated *GW* QP spectra to anion PES experiments to evaluate the performance of the *GW* approximation on both localized and delocalized states of small molecules containing 3*d* transition metals.

We found that the perturbative one-shot *G*_{0}*W*_{0}*@*PBE scheme, which is the most widely used *GW* scheme for extended systems, has a couple of problems for finite systems. Fundamentally, *G*_{0}*W*_{0}*@*PBE underestimates the IE and the 3*d*-electron binding energy by ∼1 eV and ∼2 eV, respectively, which are considerably larger than the widely reported underestimation error of ∼0.5 eV. Due to the orbital-character-dependent nonuniform underestimations of *GW* binding energies, *G*_{0}*W*_{0}*@*PBE sometimes gives the incorrect orbital order. Practically, *G*_{0}*W*_{0}*@*PBE suffers from the *GW* multisolution issue due to the large distance between QP and gKS-PBE eigenvalues and the complicated pole (peak) structure in the self-energy (the spectral function).

We found that the *G*_{0}*W*_{0} starting-point approach, *G*_{0}*W*_{0}*@*PBE*α*, can improve *G*_{0}*W*_{0}*@*PBE at the expense of introducing a couple of problems. The *G*_{0}*W*_{0} starting-point approach can give good results for the IE and the 3*d*-electron binding energy at the same time and, thus, correct the wrong orbital order produced by PBE. Also, the *G*_{0}*W*_{0} starting-point approach can mitigate the *GW* multisolution issue by reducing the distance between QP and gKS eigenvalues. However, the optimal amount of EXX in the *G*_{0}*W*_{0} starting point depends strongly on the amount of 3*d* character in molecular orbitals, leading to the strong sensitivity of 3*d*-electron binding energy to a change in the EXX amount. Thus, the optimal amount of EXX is strongly system- and property-dependent. More importantly, *G*_{0}*W*_{0}*@*PBE*α* suffers from the SCF convergence issue in open-shell systems, which is absent in *G*_{0}*W*_{0}*@*PBE.

We found that the eigenvalue self-consistency approaches, *G*_{n}*W*_{0}*@*PBE and *G*_{n}*W*_{n}*@*PBE, can improve *G*_{0}*W*_{0}*@*PBE, too. Especially, *G*_{n}*W*_{0}*@*PBE gives as good results for the IE and the 3*d*-electron binding energy as *G*_{0}*W*_{0}*@*PBE*α* without suffering from *GW* multisolution and SCF convergence issues.

We recommend *G*_{n}*W*_{0}*@*PBE because of its *practical* advantages: (i) *G*_{n}*W*_{0}*@*PBE is transferable because it gives satisfactorily accurate results for both finite and extended systems, for both closed- and open-shell systems, and for both localized and delocalized states; (ii) *G*_{n}*W*_{0}*@*PBE is predictive because it does not need any system- and property-dependent parameters; and (iii) *G*_{n}*W*_{0} is efficient and easy to use because it does not require computational and human efforts to address SCF convergence and *GW* multisolution issues.

We attribute the good performance of *G*_{n}*W*_{0}*@*PBE to the fortuitous cancellation effect: the overscreening of the Coulomb interaction due to the over-delocalized PBE wavefunction is canceled by the underscreening due to the neglect of vertex corrections. In other words, for *G*_{0}*W*_{0} applied to finite systems, PBE is a “bad” starting point in the sense that it causes a large (∼1–2 eV) and orbital-character-dependent underestimation error in electron binding energy, but for *G*_{n}*W*_{0} applied to finite and extended systems, PBE is a “good” starting point in the sense that it accidentally produces the overscreening just as much as vertex corrections do, which is missing in self-consistent *GW* schemes.

Our results in this work—(i) *G*_{0}*W*_{0}*@*PBE*α*(0.25 ≤ *α* ≤ 0.50) and *G*_{n}*W*_{0}*@*PBE give good QP energies for molecular orbitals with both weak and strong 3*d* character and (ii) the ev*GW* starting-point dependency is more related to 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 disagreement is that except for *G*_{0}*W*_{0}*@*PBE*α*(0.00 ≤ *α* ≤ 0.25), varying the self-consistency level and the starting point generally makes a small (∼0.1 eV) change in QP energy for HOMO with mainly *sp* character, which is accidentally comparable to individual or combined errors from multiple sources, such as the incomplete basis set, the linearization method in *G*_{0}*W*_{0}, and the insufficient number of eigenvalues to update in ev*GW*.

*G*_{n}*W*_{0}*@*PBE is not a conserving and starting-point-independent *GW* scheme. It is not the most accurate or efficient *GW* scheme, either. However, *G*_{n}*W*_{0}*@*PBE gives satisfactory and reliable results for a wide range of systems, such as solids with strong screening and molecules with weak screening, at moderate computational and minimal human efforts, and thus is ideal for automated mass *GW* and BSE calculations for high-throughput screening and machine learning. Further studies on the performance of more diverse *GW* schemes on larger and more complex systems containing a broader range of transition metals are needed to extend the range of applicability of the *GW* approximation.

## SUPPLEMENTARY MATERIAL

See supplementary material for more details, results, and discussion.

## ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy Grant No. DE-SC0017824. The computation for this work was performed on the high performance computing infrastructure provided by Research Computing Support Services at the University of Missouri, Columbia. This research also used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We would also like to thank Bin Shi and Meisam Rezaei for useful discussions in the earlier stage of this work.

## REFERENCES

^{−}and $ScO2\u2212$ photoelectron spectra

_{0}W

_{0}for molecular systems

_{n}

^{−}(n = 1–4) and YO

_{n}

^{−}(n = 1–5): Strong electron correlation effects in ScO

^{−}and YO

^{−}

_{y}(y = 1–3) and (TiO

_{2})

_{n}(n = 1–4)

_{2}complexes: A study of CuO

_{x}(x = 0–6) species by anion photoelectron spectroscopy

_{2}O

_{2}O

_{2}O

_{3}) as a benchmark

^{0,±}, TiO

^{0,±}, CrO

^{0,±}, and MnO

^{0,±}

_{n}and TiO

_{n}

^{−}(n = 1 − 3)

^{+}versus CuO and CuO

^{+}

^{d}-metal monocarbides and monoxides