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 3d shells (shallow and deep 3d 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 3d transition metals. We find that the G-only eigenvalue self-consistent GW scheme with W fixed to the PBE level (GnW0@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 3d-transition-metal oxide molecules. The success of GnW0@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 GnW0@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.

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 3d 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 3d-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 GW100 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 3d 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 (G0W0) scheme, the eigenvalue self-consistent GW (evGW) scheme (with two types GnW0 and GnWn, which update eigenvalues only in G and in both G and W, respectively), the QP self-consistent GW (QSGW) scheme using a static and Hermitian approximation to the GW self-energy, and the fully self-consistent GW (SCGW) 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, SCGW and QSGW systematically overestimate the bandgaps of solids,20 displaying worse performance than evGW, 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 G0W0 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 3d transition metals. To this end, we calculate the QP spectra of closed- and open-shell molecular anions with partially and completely filled 3d 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 Cu2O and ZnO, respectively, which are challenging systems for the GW method,28,29 and (iii) shallow and deep 3d 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 3d-electron binding energies, and conclude that the GnW0@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.

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 paper6 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 Π.

In Hedin’s GW approximation, the nonlocal, dynamical, and non-Hermitian self-energy Σσ at frequency ω is given by

Σσ(r,r,ω)=i2πdωeiηωGσ(r,r,ω+ω)W(r,r,ω),
(1)

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 ϵmσ and corresponding wavefunctions φmσ(r) obtained from the Hartree or mean-field approximation,

Gσ(r,r,ω)=iφiσ(r)φiσ(r)ωϵiσiη+aφaσ(r)φaσ(r)ωϵaσ+iη,
(2)

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σ, 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 χ=χ0[1vχ0]1 within the RPA, the screened Coulomb interaction, and the self-energy,

χ0=iσGσGσ,
(3)
W=v+vχ0W=v+vχ0v+vχ0vχ0v+=v+vχv,
(4)
Σσ=iGσW=iGσ(v+vχv)=Σxσ+Σcσ,
(5)

where v denotes the bare (unscreened) Coulomb interaction v(r, r′) = 1/|rr′|, Σ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 Σcσ(ω) are dynamic, whereas v and Σxσ 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

ϵmG0W0,σ=ϵmσ+φmσ|ReΣσ(ϵmG0W0,σ)vxcσ|φmσ,
(6)

where ϵmG0W0,σ are the perturbative one-shot GW QP energies and vxcσ is the xc potential. Experimentally, ϵmG0W0,σ correspond to vertical IEs and EAs in PES and IPES, respectively. Theoretically, ϵmG0W0,σ 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σ,

Ammσ(r,r,ω)=1π|ImGmmσ(r,r,ω)|,
(7)

where Ammσ are the diagonal elements of the spectral function, Gmmσ are the diagonal elements of Green’s function, and Im represents the imaginary part. Note that Ammσ 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 ϵmσ, where m = i or a, by ϵmσ+φmσ|Σσ(ω)vxcσ|φmσ, one can find that Ammσ have Lorentzian peaks at

ω=ϵmσ+φmσ|ReΣσ(ω)vxcσ|φmσ,
(8)

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 ϵmG0W0,σ, 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.

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) ϕμ,

φmσ(r)=μCμmσϕμ(r),
(9)

where Cμmσ 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 scheme31 for (semi-)local functionals, hybrid functionals, and HF],

HσCσ=SCσϵσ,
(10)

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

Sμν=drϕμ(r)ϕν(r),
(11)

and Hσ is the Hamiltonian matrix with elements

Hμνσ=Tμν+Vext,μν+JμναKμνσ+(1α)Vx,μνPBE,σ+Vc,μνPBE,σ,
(12)

where T, Vext, J, and Kσ are the kinetic energy, external potential energy, Hartree, and Fock exchange terms, respectively, Vxσ and Vcσ 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

Jμν=λτ(μν|λτ)σDλτσ,
(13)

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

(μν|λτ)=drdrϕμ(r)ϕν(r)1|rr|ϕλ(r)ϕτ(r),
(14)

and Dσ is the density matrix with elements

Dμνσ=mfmσCμmσCνmσ,
(15)

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

Kμνσ=λτDλτσ(μλ|τν).
(16)

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

ρσ(r)=μνDμνσϕμ(r)ϕν(r).
(17)

The gKS equation in Eq. (10) (the restricted Roothaan-Hall or unrestricted Pople-Nesbet equations) should be solved using the SCF method because J, Kσ, Vxσ, and Vcσ in Eq. (12) depend on the density matrix in Eq. (15), as shown in Eqs. (13), (16), and (17).

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

A BB AXsYs=ΩsXsYs,
(18)

where A and B are the resonant and coupling matrices, respectively, and Ωs and (Xs, Ys) are the eigenvalues (the neutral two-particle excitation energies) and corresponding eigenvectors, respectively. The matrix elements in A and B are given by

Aiaσjbσ=(ϵaσϵiσ)δijδabδσσ+(iaσ|jbσ)+fxc,iaσjbσ,
(19)
Biaσjbσ=(iaσ|bjσ)+fxc,iaσbjσ,
(20)

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

(iaσ|jbσ)=drdrφiσ(r)φaσ(r)1|rr|φjσ(r)φbσ(r).
(21)

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

(iaσ|jbσ)=μνλτCμiσCνaσCλjσCτbσ(μν|λτ),
(22)

which scales as O(N5). 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 (Xs, Ys). 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 (Xs, Ys), 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 Σxσ and Σcσ(ω), respectively, whose diagonal matrix elements are given by

φmσ|Σxσ|φmσ=i(miσ|imσ),
(23)
φmσ|Σcσ(ω)|φmσ=iswmiσswmiσsωϵiσ+Ωsiη+aswmaσswmaσsωϵaσΩs+iη,
(24)

where i runs over occupied states, a runs over empty states, s runs over all excitations, and wmnσs are given by

wmnσs=iaσ(mnσ|iaσ)(Xiaσs+Yiaσs).
(25)

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.

In unrestricted HF and KS calculations for open-shell systems, the expectation value of the total angular momentum ⟨S2⟩ is given by

S2=S(S+1)+NiNjN|φi|φj|2,
(26)

where N and N are the numbers of ↑- and ↓-spin electrons, respectively, and S is (NN)/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.

In Gaussian-based GW, the 4-center integrals (μν|λτ) in Eq. (14), which scale as O(N4), 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

(μν|λτ)PQ(μν|P)(P|Q)1(Q|λτ),
(27)

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, Σxσ, and Σcσ(ω) 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).

TABLE I.

Number of occupied and empty states for the ↑-spin channel (Nocc and Nemp, respectively) used in GW calculations. FC means the frozen-core approximation. AE and ECP mean all electron and effective core potential, respectively. CN means the cardinal number.

Nemp
AnionPotentialFCNoccCN = 2CN = 3CN = 4CN = 5
ScO AE Yes 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
AnionPotentialFCNoccCN = 2CN = 3CN = 4CN = 5
ScO AE Yes 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.

In this work, we used three methods to obtain ϵmG0W0, as is practically impossible to obtain unique, correct, and accurate ϵmG0W0 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),

ϵmG0W0ϵm+Zmlinearφm|Σ(ϵm)vxc|φmϵmG0W0,linear,
(28)

where ϵmG0W0,linear is the perturbative one-shot QP energy obtained from the linearization, and Zmlinear is the QP renormalization factor for the linearization,

Zmlinear(ϵm)=11ωφm|Σ(ω)|φm|ω=ϵm,
(29)

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 G0W0 method, different frequency points give similar results for ϵmG0W0,linear because PPA makes ⟨φmc(ω)|φm⟩ in Eq. (24) smooth around ϵm by drastically reducing the number of self-energy poles.39 However, in the full-frequency G0W0 method, different frequency points can give very different results for ϵmG0W0,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,

Zm(ϵmG0W0)=11ωφm|Σ(ω)|φm|ω=ϵmG0W0.
(30)

The derivative of the self-energy is evaluated at ω = ϵm and ω=ϵmG0W0 in Zmlinear and Zm, respectively. Generally, Zm is smaller than Zmlinear because the self-energy has a steeper slope at ω=ϵmG0W0 than at ω = ϵm. Zm 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 ϵmG0W0 obtained from the graphical solution as ϵmG0W0,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 Amm in Eq. (7). In this work, we refer to ϵmG0W0 obtained from the spectral function as ϵmG0W0,spect and define Zmspect by replacing Zmlinear and ϵmG0W0,linear in Eq. (28) by Zmspect and ϵmG0W0,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 Zm in Eq. (30) because Zm represents the spectral weight, as explained above.

As introduced in Sec. I, there are various levels of self-consistency in the GW approximation (from the lowest to the highest): G0W0, GnW0, GnWn, QSGW, and SCGW. In this work, we used GnW0 and GnWn for simplicity, efficiency, and stability. GnW0 updates only gKS eigenvalues in Gσ [ϵiσ and ϵaσ in Eq. (24)], whereas GnWn updates gKS eigenvalues in Gσ and A [ϵiσ and ϵaσ in Eq. (19)] as well as Casida eigenvalues in Ws in Eq. (24)]. Therefore, GnWn is computationally more expensive than GnW0 by the time to build and completely diagonalize the RPA Casida matrix in Eq. (18). Note that GnWn can be viewed as a diagonal approximation to QSGW.

In this work, we obtained GnW0 and GnWn QP energies (ϵmGnW0 and ϵmGnWn, respectively) by iterating the recurrence relations (n ≥ 3),

ϵmevGW,1=ϵm+ZevGWφm|Σ(ϵm)vxc|φm,
(31)
ϵmevGW,2=ϵmevGW,1+ZevGWφm|Σ(ϵmevGW,1)Σ(ϵm)|φm,
(32)
ϵmevGW,n=ϵmevGW,n1+ZevGWφm|Σ  ϵmevGW,n1  Σ  ϵmevGW,n2  |φm,
(33)

where ϵmevGW,n is ϵmGnW0 or ϵmGnWn, and ZevGW = 1. Summing up Eqs. (31)–(33), we get

ϵmevGW,n=ϵm+φm|Σ(ϵmevGW,n1)vxc|φm,
(34)

which we refer to as the evGW QP equation in this work. When the evGW convergence is reached (ϵmevGW,n=ϵmevGW,n1=ϵmevGW,, where ϵmevGW, are converged evGW QP energies), the evGW QP equation in Eq. (34) becomes similar to the G0W0 QP equation in Eq. (6).

Whereas most GW codes use 0 < ZevGW < 1,7,17,40 MOLGW uses ZevGW = 1. Even though we implemented evGW with 0 < Z < 1 into MOLGW, we adopted evGW with Z = 1 in this work for a few reasons. First, Eq. (33) shows that converged evGW QP energies (ϵmevGW,) 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, evGW with 0 < Z < 1 is suited for a simplified evGW 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, evGW with Z = 1 has no variant and does not need a QP equation solver, but evGW with 0 < Z < 1 has multiple variants, depending on the choice of QP equation solvers. For example, two evGW 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).

ZevGW = 1 and the efficiency comparison of evGW and G0W0 are discussed in the supplementary material.

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-pVnZ (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,

ExcPBEα,σ=αExHF,σ+(1α)ExPBE,σ+EcPBE,σ,
(35)

where ExHF,σ, ExPBE,σ, and EcPBE,σ 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.

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 3d-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 3d-electron binding energies, where scalar relativistic effects are considerable (∼10% of the experimental 3d-electron binding energy). Note that Ref. 13 used only AE even though the GW100 benchmark set contains Ag2, Cu2, 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 3d 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 ⟨S2⟩ 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.

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,

Em=a+bNBF,
(36)
Em=a+bCN3,
(37)

where Em are correlation or mth QP energies, a and b are fitting parameters, NBF 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.

FIG. 1.

Effect of the fitting function type on the GW complete basis set (CBS) limit using G0W0@PBEα(α = 1.00) QP HOMO energies of CuO. a + b/NBF is used on the left, while a + b/CN2 and a + b/CN3 are used on the right. NBF and CN represent the number of basis functions and the cardinal number, respectively. G0W0@PBEα(α = 1.00) QP HOMO energies are obtained from gKS-PBEα(α = 1.00) HOMO-1, which corresponds to gKS-PBE HOMO (see text).

FIG. 1.

Effect of the fitting function type on the GW complete basis set (CBS) limit using G0W0@PBEα(α = 1.00) QP HOMO energies of CuO. a + b/NBF is used on the left, while a + b/CN2 and a + b/CN3 are used on the right. NBF and CN represent the number of basis functions and the cardinal number, respectively. G0W0@PBEα(α = 1.00) QP HOMO energies are obtained from gKS-PBEα(α = 1.00) HOMO-1, which corresponds to gKS-PBE HOMO (see text).

Close modal
FIG. 2.

Effect of the EXX amount on the GW complete basis set (CBS) limit using G0W0@PBE (left) and G0W0@PBEα(α = 1.00) (right) QP HOMO energies of ScO. NBF represents the number of basis functions.

FIG. 2.

Effect of the EXX amount on the GW complete basis set (CBS) limit using G0W0@PBE (left) and G0W0@PBEα(α = 1.00) (right) QP HOMO energies of ScO. NBF represents the number of basis functions.

Close modal

In this work, we obtained gKS and GW CBS results using the fitting function in Eq. (36) with CN = 2,3,4,5. We chose Eq. (36) not because it is superior to Eq. (37) but because it can be used by both Gaussian- and PW-based GW implementations.

2. Number of empty states

By enabling MOLGW 1.F to support the largest 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 Σcσ(ω), 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. G0W0 quasiparticle energy

A full-frequency G0W0 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 ϵmG0W0 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 ϵmG0W0 shows that the graphical-solution method using η = 0.001 hartree and the linearization method randomly give incorrect solutions. Therefore, we obtained ϵmG0W0 from the graphical-solution or spectral-function method using η = 0.002 or 0.005 hartree. When ϵmG0W0,graph and ϵmG0W0,spect are different, we manually selected a correct solution by analyzing Σc(ω) and A(ω).

A large distance between ϵm and ϵmG0W0 and multiple self-energy poles between them make it difficult to obtain correct and accurate ϵmG0W0 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 G0W0@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.

FIG. 3.

Comparison of three QP equation solvers: linearization, graphical-solution, and spectral-function methods. (Top left) All the three methods give correct solutions. (Top right) The graphical-solution method can give an incorrect solution. (Bottom left) The linearization method can give an incorrect solution. (Bottom right) The spectral-function method can give an incorrect solution. EKS represents the gKS-PBE eigenvalue. ElinQP, EgraQP, and EspeQP represent G0W0@PBE QP energies obtained from linearization, graphical-solution, and spectral-function methods, respectively. The green straight line is a tangent to the red curve at ω = EKS. Except for A(ω) at the bottom right, all results are obtained from CN = 2 no-RI AE.

FIG. 3.

Comparison of three QP equation solvers: linearization, graphical-solution, and spectral-function methods. (Top left) All the three methods give correct solutions. (Top right) The graphical-solution method can give an incorrect solution. (Bottom left) The linearization method can give an incorrect solution. (Bottom right) The spectral-function method can give an incorrect solution. EKS represents the gKS-PBE eigenvalue. ElinQP, EgraQP, and EspeQP represent G0W0@PBE QP energies obtained from linearization, graphical-solution, and spectral-function methods, respectively. The green straight line is a tangent to the red curve at ω = EKS. Except for A(ω) at the bottom right, all results are obtained from CN = 2 no-RI AE.

Close modal

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 ϵmG0W0 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 ϵmG0W0 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 ϵmG0W0 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 ωϵmG0W0. 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 ϵmG0W0 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 [G0W0@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 ϵmG0W0,linear=2.9 eV causes a large overestimation error of 0.8 eV in G0W0@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 ϵmG0W0 for deep states not only automatically but also manually.

We conclude this section by summarizing several guidelines to obtain a reliable and reproducible G0W0@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 |ϵmG0W0-ϵm|, which generally decreases with the amount of EXX and increases with the depth of the mth state. For example, when calculating G0W0@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 G0W0@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 G0W0 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 G0W0@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 G0W0@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 G0W0@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 G0W0@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 Zm 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 Zm, 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 G0W0 starting point and evGW 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 Zm (∼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 φmσ|Σcσ(ω)|φmσ in Eq. (24) diverge at ω=ϵiσΩs and ω=ϵaσ+Ω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. GnW0 and GnWn quasiparticle energy

In this work, we used only η = 0.001 hartree for evGW because it is small enough to obtain the convergence of evGW QP energies with respect to η within ∼0.01 eV. The convergence test results for evGW QP energies with respect to the iteration number are shown in the supplementary material. QSGW and our evGW are quasiparticle-only GW schemes with no spectral weight transfer (Z = 1), and GnWn is a diagonal approximation to QSGW. Therefore, we compared the convergence behaviors of QSGW and our evGW and found a couple of similarities and differences between them.

First, the evGW 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 evGW. Unlike evGW, QSGW 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 evGW. For example, we observed that evGW QP energies for HOMO of CuO depend more strongly on the EXX amount in the evGW starting point than those for HOMO of ScO. We attribute this to different amounts of 3d character in HOMOs of ScO and CuO (6% and 23%, respectively, as shown in Table II). In other words, as the 3d character in MO increases, the starting-point dependency of evGW increases. We also observed that GnWn has a weaker (stronger) starting-point dependency for HOMO of ScO (CuO) than GnW0. Our observations for ScO are consistent with Ref. 41, which studied the evGW starting-point dependency using small water clusters and concluded that as the evGW self-consistency level increases from GnW0 and GnWn, the evGW starting-point dependency decreases. However, our observations for CuO are not consistent with this conclusion. This is likely because CuO has strong 3d character in HOMO, whereas ScO and small water clusters do not. This orbital-character-dependent starting-point dependency of evGW may be related to conflicting results for QSGW in the literature: Ref. 58 showed the starting-point independency of QSGW using a small sp−bonded molecule (CH4), while Ref. 59 showed the strong starting-point dependency of QSGW using a d solid (α-Fe2O3).

TABLE II.

TM 3d character in molecular orbitals of TMO anions, obtained from NWChem gKS-PBE and gKS-PBEα(α = 1.00) results with CN = 2 no-RI AE using the Mulliken population analysis. The gKS-PBE orbital order is used for gKS-PBEα(α = 1.00) molecular orbitals (see text). Bold numbers are used to highlight entire TM 3d character.

TiOZnO
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   
TiOZnO
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   

In this section, we compare our GW calculations to anion PES experiments,24–27 focusing especially on the first IE, the lowest 3d-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, G0W0@PBEα calculations as α is varied in steps of 0.25 from 0 to 1), and then, we discuss eigenvalue self-consistent GW (GnW0 and GnWn) 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 3d shells, while CuO and ZnO have completely filled 3d shells. Third, TiO has a shallow 3d state, but CuO and ZnO have deep 3d states. Fourth, 3d-electron photodetachment transitions are observed in TiO and CuO, but not in ScO and ZnO. Fifth, CuO has strong 3d character in HOMO, but ScO, TiO, and ZnO have weak 3d 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 3d 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 3d 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.

FIG. 4.

Contour plots of molecular orbitals of TMO anions, obtained from MOLGW gKS-PBE results with CN = 2 no-RI AE using VESTA.60 Red boxes are used to highlight entire TM 3d character.

FIG. 4.

Contour plots of molecular orbitals of TMO anions, obtained from MOLGW gKS-PBE results with CN = 2 no-RI AE using VESTA.60 Red boxes are used to highlight entire TM 3d character.

Close modal
TABLE III.

Experimental (PES) and calculated (AE GW) electron binding energies of TMO anions (in eV). MAE represents the mean absolute error. Bold numbers are used to highlight 3d-electron binding energies.

G0W0@PBEαGnW0@GnWn@
StateExp.α = 0.00α = 0.25α = 0.50α = 0.75α = 1.00PBEPBEOthers
ScO (1Σ+          
HOMO 1Σ+ 1.35a 0.51 1.15 1.45 1.59 1.63 1.21 1.35 1.28b, 1.26c, 1.19d 
HOMO-1 2Δ 3.10a 3.30 4.77 5.45 5.72 5.82 5.34 6.08 2.41b, 2.78e, 3.31f 
HOMO-2 2Π 3.40a 3.42 4.81 5.40 5.61 5.63 5.39 6.17 3.34b, 3.24e, 3.44f 
MAE   0.84g 0.18g 0.10g 0.24g 0.28g 0.14g 0.00g  
TiO (2Δ)           
↑-HOMO 1Σ+ 2.00h 0.26 2.24 3.70 4.79 5.63 1.83 2.74 2.39b, 2.37e, 2.34f 
↑-HOMO-1 1Δ 1.73h 0.53 1.28 1.65 1.83 1.95 1.38 1.61 1.88b, 1.72f 
↓-HOMO 3Δ 1.30h 0.31 1.00 1.29 1.43 1.55 1.06 1.21 1.19b, 1.18c, 1.14i 
MAE   1.31 0.33 0.60 1.01 1.37 0.25 0.32  
CuO (1Σ+          
HOMO 2Π 1.78j 0.40 1.40 1.58 1.40 0.97 2.19 3.17 1.55b, 1.52c, 0.46k 
HOMO-1 2Σ+ 2.75j 1.39 2.17 2.23 1.98 1.56 2.66 3.58 2.96b, 2.86e, 1.60k, 2.78l, 2.47m, 2.81n 
HOMO-2  4.50o 2.18 4.05 4.60 4.56 4.25 4.88 6.38 4.07b, 4.01e, 4.58l, 4.50m 
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.09p 0.91 1.91 2.20 2.23 2.00 2.11 2.57 2.19b, 2.33c, 2.29q, 2.10r, 1.06k 
↑-HOMO-1 1Π 2.71p 1.36 2.43 3.03 3.72 4.79 2.77 3.73 2.62b, 1.43k 
↑-HOMO-2   3.24 4.04 4.89 5.64 6.90 4.77 5.66 3.50k 
↓-HOMO 3Π 2.40p 1.17 2.17 2.64 3.15 4.02 2.65 3.48 2.41b, 1.20k 
↓-HOMO-1 3Σ+ 3.96p 2.71 3.35 3.47 3.31 3.00 4.11 4.51 4.15b, 2.89k 
MAE   1.25 0.33 0.29 0.64 1.19 0.12 0.78  
G0W0@PBEαGnW0@GnWn@
StateExp.α = 0.00α = 0.25α = 0.50α = 0.75α = 1.00PBEPBEOthers
ScO (1Σ+          
HOMO 1Σ+ 1.35a 0.51 1.15 1.45 1.59 1.63 1.21 1.35 1.28b, 1.26c, 1.19d 
HOMO-1 2Δ 3.10a 3.30 4.77 5.45 5.72 5.82 5.34 6.08 2.41b, 2.78e, 3.31f 
HOMO-2 2Π 3.40a 3.42 4.81 5.40 5.61 5.63 5.39 6.17 3.34b, 3.24e, 3.44f 
MAE   0.84g 0.18g 0.10g 0.24g 0.28g 0.14g 0.00g  
TiO (2Δ)           
↑-HOMO 1Σ+ 2.00h 0.26 2.24 3.70 4.79 5.63 1.83 2.74 2.39b, 2.37e, 2.34f 
↑-HOMO-1 1Δ 1.73h 0.53 1.28 1.65 1.83 1.95 1.38 1.61 1.88b, 1.72f 
↓-HOMO 3Δ 1.30h 0.31 1.00 1.29 1.43 1.55 1.06 1.21 1.19b, 1.18c, 1.14i 
MAE   1.31 0.33 0.60 1.01 1.37 0.25 0.32  
CuO (1Σ+          
HOMO 2Π 1.78j 0.40 1.40 1.58 1.40 0.97 2.19 3.17 1.55b, 1.52c, 0.46k 
HOMO-1 2Σ+ 2.75j 1.39 2.17 2.23 1.98 1.56 2.66 3.58 2.96b, 2.86e, 1.60k, 2.78l, 2.47m, 2.81n 
HOMO-2  4.50o 2.18 4.05 4.60 4.56 4.25 4.88 6.38 4.07b, 4.01e, 4.58l, 4.50m 
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.09p 0.91 1.91 2.20 2.23 2.00 2.11 2.57 2.19b, 2.33c, 2.29q, 2.10r, 1.06k 
↑-HOMO-1 1Π 2.71p 1.36 2.43 3.03 3.72 4.79 2.77 3.73 2.62b, 1.43k 
↑-HOMO-2   3.24 4.04 4.89 5.64 6.90 4.77 5.66 3.50k 
↓-HOMO 3Π 2.40p 1.17 2.17 2.64 3.15 4.02 2.65 3.48 2.41b, 1.20k 
↓-HOMO-1 3Σ+ 3.96p 2.71 3.35 3.47 3.31 3.00 4.11 4.51 4.15b, 2.89k 
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 G0W0@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 G0W0@PBE0 method.

Throughout this work, we use only gKS-PBE TM 3d character except when we discuss the subtle competition between direct and indirect relativistic effects because the EXX amount has a small effect on TM 3d 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 G0W0@PBEα(α = 1.00) HOMO of CuO, as shown in the supplementary material and Fig. 6. In this work, gKS-PBE and G0W0@PBE were found to have the same orbital order.

Figures 5 and 6 show PES and G0W0@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 G0W0@PBEα(0.00 ≤ α ≤ 1.00) results are in perfect agreement with experiment, (ii) G0W0@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) G0W0@PBEα(0.25 ≤ α ≤ 0.50) reduces it to ∼0.1 eV. In the following, we analyze each TMO anion individually.

FIG. 5.

Effect of the EXX amount in the G0W0 starting point on the electronic structure of ScO (left) and TiO (right). A Gaussian distribution function with a smearing width of 0.1 eV is used to broaden the spectra.

FIG. 5.

Effect of the EXX amount in the G0W0 starting point on the electronic structure of ScO (left) and TiO (right). A Gaussian distribution function with a smearing width of 0.1 eV is used to broaden the spectra.

Close modal
FIG. 6.

Effect of the EXX amount in the G0W0 starting point on the electronic structure of CuO (left) and ZnO (right). A Gaussian distribution function with a smearing width of 0.1 eV is used to broaden the spectra.

FIG. 6.

Effect of the EXX amount in the G0W0 starting point on the electronic structure of CuO (left) and ZnO (right). A Gaussian distribution function with a smearing width of 0.1 eV is used to broaden the spectra.

Close modal

1. ScO

Scandium is the first transition metal and has only one 3d electron. DFT and CCSD(T) calculations in Refs. 10 and 11 confirmed the ground state of ScO as 1Σ+ (8σ23π49σ2), correcting the wrongly assumed state 3Δ (8σ23π49σ11δ) in Ref. 24. There is no 3d peak or band in the PES spectrum of ScO, and the top three valence molecular orbitals have weak Sc 3d 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, G0W0@PBE binding energies slightly overestimate PES ones (by 0.20 and 0.02 eV for HOMO-1 and HOMO-2, respectively), whereas G0W0@PBEα(0.25 ≤ α ≤ 1.00) binding energies significantly overestimate PES ones by ∼2 eV [e.g., for HOMO-1, G0W0@PBEα(α = 0.25) and G0W0@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, G0W0@PBE performs better than G0W0@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σ23π49σ2 (1Σ+ ScO) to 8σ23π410σ (B2Σ+ ScO) and to 8σ23π41δ (A2Δ ScO) states, respectively, which GW calculations for quasiparticle excitations cannot account for. In other words, the seemingly excellent agreement between G0W0@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, G0W0@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 G0W0@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 3d 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 3d electron in the δ shell. The transition of the 3d 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 G0W0@PBEα QP spectrum of TiO, ↑-HOMO is of entirely Ti 3d character, as shown in Table II.

The right panel of Fig. 5 clearly shows that G0W0@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 G0W0@PBEα binding energy to a change in α causes a couple of problems. First, G0W0@PBE underestimates the IE and the 3d-electron binding energy of TiO nonuniformly (by 0.99 and 1.74 eV, respectively), leading to the wrong orbital order. In other words, G0W0 does not correct the wrong orbital order produced by PBE. Second, the G0W0 starting-point approach does not give accurate results for both the IE and the 3d-electron binding energy of TiO at the same time. For example, G0W0@PBEα(α = 0.50) gives a better result for the IE of TiO by 0.29 eV, but a worse result for the 3d-electron binding energy of TiO by 1.46 eV, than G0W0@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 G0W0@PBEα IE of both ScO and TiO: For both ScO and TiO, G0W0@PBEα(0.25 ≤ α ≤ 1.00) reduces the underestimation of IE by G0W0@PBE from ∼1 eV to ∼0.1 eV [e.g., G0W0@PBEα(α = 0.25) reduces the difference in IE between PES and G0W0@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 G0W0@PBEα 3d-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 3d electrons. DFT calculations in Ref. 10 confirmed the ground state of CuO as 1Σ+ with the electron configuration of 3d102224. 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 3d electrons (343432) from 1Σ+ CuO 3d102224 to Z CuO 3d92224 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 G0W0@PBEα QP spectrum of CuO, HOMO-2 is of entirely Cu 3d character, as shown in Table II.

In the left panel of Fig. 6, we see that G0W0@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 3d character than HOMO-2, as shown in Table II, and G0W0@PBEα(0.50 ≤ α ≤ 0.75) gives good results for the IE and the 3d-electron binding energy (corresponding to HOMO and HOMO-2, respectively) of CuO at the same time. Scalar relativistic effects in ECP reduce G0W0@PBEα(α = 0.50) and G0W0@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 G0W0@PBEα(0.50 ≤ α ≤ 0.75) gives good results for the 3d-electron binding energy of CuO. Like TiO, the orbital-character-dependent sensitivity of G0W0@PBEα binding energy to a change in α causes G0W0@PBE to underestimate the IE and the 3d-electron binding energy of CuO nonuniformly (by 1.38 and 2.32 eV, respectively). G0W0@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 3d 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 3d 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 3d electrons generally do not participate in bonding.70 DFT calculations in Ref. 10 confirmed the ground-state electron configuration of ZnO as 2Σ+ 10σ19σ24π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 3d 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 G0W0@PBEα QP spectrum of ZnO have weak Zn 3d character, as shown in Table II.

The right panel of Fig. 6 shows that unlike CuO, G0W0@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 3d character and thus their G0W0@PBEα binding energies have similar sensitivity to a change in α. G0W0@PBEα(α = 0.50) gives good results for the IE and the orbital order of ZnO at the same time. Unlike CuO, G0W0@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 CuO2 molecule).

Figure 7 shows PES and evGW QP spectra of ScO, TiO, CuO, and ZnO. In the following, we analyze GnW0@PBE and GnWn@PBE results individually.

FIG. 7.

Effect of the evGW self-consistency level on the electronic structure of ScO (top left), TiO (top right), CuO (bottom left), and ZnO (bottom right). A Gaussian distribution function with a smearing width of 0.1 eV is used to broaden the spectra.

FIG. 7.

Effect of the evGW self-consistency level on the electronic structure of ScO (top left), TiO (top right), CuO (bottom left), and ZnO (bottom right). A Gaussian distribution function with a smearing width of 0.1 eV is used to broaden the spectra.

Close modal

In Fig. 7, we see that as the evGW self-consistency level increases from G0W0 to GnWn, GW binding energies always increase, but this occurs at different rates: As the evGW self-consistency level increases from G0W0 to GnW0, 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 GnW0 to GnWn, 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. G0W0@PBE always underestimates electron binding energies, whereas GnWn@PBE generally overestimates them. GnW0@PBE binding energies are always in between G0W0@PBE and GnWn@PBE ones and generally close to experiment. In other words, G0W0@PBE and GnWn@PBE act as lower and upper bounds for GnW0@PBE, generally producing overscreenings and underscreenings, respectively. This trend of the evGW self-consistency approach in electronic structure of molecules is also observed in the band structure of solids.17 

We also see that the evGW self-consistency has a strong effect on GW binding energies for molecular orbitals with strong 3d character (e.g., ↑-HOMO of TiO and HOMO-2 of CuO). For example, GnW0@PBE reduces the underestimation errors of G0W0@PBE in the IE and the 3d-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, GnW0@PBE corrects the wrong G0W0@PBE orbital order in TiO. Another example is that GnW0@PBE gives small (∼0.1 eV) errors in electron binding energies for all valence molecular orbitals of ZnO, which are uniformly underestimated by G0W0@PBE by ∼1 eV due to similarly weak Zn 3d character. For ZnO, G0W0@PBE and GnW0@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 GnW0@PBE and GnWn@PBE binding energies compared to other TMO anions. This trend is not associated with scalar relativistic effects in ECP, which reduce GnW0@PBE and GnWn@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 3d character in molecular orbitals of CuO. For example, CuO has a larger difference between GnW0@PBE and GnWn@PBE IEs than ScO (0.98 and 0.14 eV, respectively) possibly because CuO has stronger 3d character in HOMO than ScO (23% and 6%, respectively, as shown in Table II).

From our results presented so far, it appears that both G0W0 starting-point and evGW self-consistent approaches can, in principle, be good GW methods for finite systems: both G0W0@PBEα(0.25 ≤ α ≤ 0.50) and GnW0@PBE can reduce the large and orbital-character-dependent nonuniform errors for electron binding energies of TMO anions produced by G0W0@PBE with respect to experiment from ∼1–2 eV to ∼0.1–0.5 eV. Reference 18 obtained similar results for extended systems: both G0W0@PBEα(α = 0.25) and GnW0@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 G0W0@PBEα(α = 0.25) or GnW0@PBE because they give similar results, but (ii) for efficiency, one may want to choose G0W0@PBEα(α = 0.25) over GnW0@PBE because the former is computationally cheaper than the latter. However, in the case of molecular systems, we argue that GnW0@PBE has several practical advantages over G0W0@PBEα.

First, GnW0@PBE does not contain system-dependent adjustable parameters. Unlike extended systems, there is no unique amount of EXX for the G0W0 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, GnW0@PBE is transferable between finite and extended systems. GnW0@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, GnW0@PBE may be applicable to solid-molecule hybrid systems such as molecular junctions and molecules adsorbed on solid surfaces.71 Also, GnW0@PBE may be used for the study of quantum size effects in clusters because it is independent of the cluster size. Third, GnW0@PBE is easy to use and reliable. Unlike PBEα(0.00 < α ≤ 1.00), PBE is safe from the SCF convergence issue, and unlike G0W0, evGW with Z = 1 is immune to the GW multisolution issue. Therefore, GnW0@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, GnW0@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 GnW0@PBE cheaper. Two desirable properties come from the GW part. GnW0 (as well as GnWn) gives faster and more stable GW convergence than QSGW and depends more weakly on the choice of η (e.g., we used a single value of η for evGW in this work), as discussed in Sec. III C 4. Also, GnW0 is cheaper than GnWn, as pointed out in Ref. 17 and discussed in Sec. II G. In fact, GnW0 is the cheapest self-consistent GW scheme.

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

Some of our results for the performance of G0W0 starting-point and evGW 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 G0W0 starting-point approach. Table IV summarizes a few selected results for the optimal amount of EXX in the G0W0 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 G0W0 IE (∼0.06 eV). There are a couple of other factors that have a stronger effect on G0W0 results than implementation differences. One factor is the choice of system and property. As shown in Sec. IV A, G0W0@PBEα(0.25 ≤ α ≤ 1.00) IEs of sp systems are slightly different (by ∼0.1 eV). Most existing G0W0 studies used the IE of sp-bonded systems to determine the optimal amount of EXX in the G0W0 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 G0W0 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 G0W0 starting point and thus is likely to produce the wide range of amounts that exist in the literature.

TABLE IV.

Optimal amount of EXX in the G0W0 starting point for gas-phase small molecules (highlighted in bold).

ReferenceKö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 spa 29 closed-shell sp 34 closed-shell spa 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 ΔSCFc QSGWd Experiment Experiment 
ω integration Contour Fully analytic Fully analytic Fully analytic Fully analytic 
QP equation deformatione Linearization Spectral function Linearizationf and Graphical solution and 
 Linearization   Spectral functiong 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 PAWh AE 
RI Used Not used Used Not applicablei Not used 
ReferenceKö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 spa 29 closed-shell sp 34 closed-shell spa 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 ΔSCFc QSGWd Experiment Experiment 
ω integration Contour Fully analytic Fully analytic Fully analytic Fully analytic 
QP equation deformatione Linearization Spectral function Linearizationf and Graphical solution and 
 Linearization   Spectral functiong 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 PAWh AE 
RI Used Not used Used Not applicablei 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 QSGW 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 evGW self-consistency approach and discuss the origin of apparently conflicting evGW results for IE and starting-point dependency. First, Ref. 7 reported that the GnWn 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 GnW0@PBE gives satisfactory results for the electronic structure (including the IE) of small 3d molecules. A comparison of evGW 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 evGW binding energy to a change in evGW self-consistency level. As shown in Sec. IV B, GnW0@PBE and GnWn@PBE binding energies are slightly different for delocalized HOMO with weak 3d character by ∼0.1 eV but significantly different for localized HOMO with strong 3d 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 evGW IE, as shown in Sec. III C. Accordingly, they are most likely not the reason for the large (0.98 eV) difference between GnW0@PBE and GnWn@PBE IEs of CuO. Second, Ref. 41 reported that in small water clusters, as the evGW self-consistency level increases, the evGW starting-point dependency decreases, whereas we found in Sec. III C 4 that in TMO anions, GnWn sometimes depends more strongly on the starting point than GnW0. As mentioned in Sec. III C 4, we believe that the orbital character influences the evGW starting-point dependency: for molecular orbitals with strong (weak) 3d character, GnWn depends more strongly (weakly) on the starting point than GnW0. Overall, without molecular orbitals with strong 3d character (e.g., HOMO of CuO), our evGW 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 3d character; (ii) we applied 15 different starting-point–self-consistency hybrid GW schemes (G0W0, GnW0, and GnWn; 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 3d-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 (G0W0@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 3d-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 (G0W0@PBE0 and GnW0@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).

FIG. 8.

Effect of the GW starting point and the evGW self-consistency level (G0W0, GnW0, and GnWn on the left, middle, and right, respectively) on the electronic structure of ScO (top) and TiO (bottom). Ebindexp and Ebindcal represent experimental and calculated electron binding energies, respectively. Dashed and solid lines track sp- and d-electron binding energies, respectively. αopt represents an optimal fraction of EXX in the GW starting point.

FIG. 8.

Effect of the GW starting point and the evGW self-consistency level (G0W0, GnW0, and GnWn on the left, middle, and right, respectively) on the electronic structure of ScO (top) and TiO (bottom). Ebindexp and Ebindcal represent experimental and calculated electron binding energies, respectively. Dashed and solid lines track sp- and d-electron binding energies, respectively. αopt represents an optimal fraction of EXX in the GW starting point.

Close modal

In summary, we calculated the electronic structure of closed- and open-shell molecular anions with partially and completely filled 3d shells (shallow and deep 3d 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 3d transition metals.

We found that the perturbative one-shot G0W0@PBE scheme, which is the most widely used GW scheme for extended systems, has a couple of problems for finite systems. Fundamentally, G0W0@PBE underestimates the IE and the 3d-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, G0W0@PBE sometimes gives the incorrect orbital order. Practically, G0W0@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 G0W0 starting-point approach, G0W0@PBEα, can improve G0W0@PBE at the expense of introducing a couple of problems. The G0W0 starting-point approach can give good results for the IE and the 3d-electron binding energy at the same time and, thus, correct the wrong orbital order produced by PBE. Also, the G0W0 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 G0W0 starting point depends strongly on the amount of 3d character in molecular orbitals, leading to the strong sensitivity of 3d-electron binding energy to a change in the EXX amount. Thus, the optimal amount of EXX is strongly system- and property-dependent. More importantly, G0W0@PBEα suffers from the SCF convergence issue in open-shell systems, which is absent in G0W0@PBE.

We found that the eigenvalue self-consistency approaches, GnW0@PBE and GnWn@PBE, can improve G0W0@PBE, too. Especially, GnW0@PBE gives as good results for the IE and the 3d-electron binding energy as G0W0@PBEα without suffering from GW multisolution and SCF convergence issues.

We recommend GnW0@PBE because of its practical advantages: (i) GnW0@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) GnW0@PBE is predictive because it does not need any system- and property-dependent parameters; and (iii) GnW0 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 GnW0@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 G0W0 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 GnW0 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) G0W0@PBEα(0.25 ≤ α ≤ 0.50) and GnW0@PBE give good QP energies for molecular orbitals with both weak and strong 3d character and (ii) the evGW 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 G0W0@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 G0W0, and the insufficient number of eigenvalues to update in evGW.

GnW0@PBE is not a conserving and starting-point-independent GW scheme. It is not the most accurate or efficient GW scheme, either. However, GnW0@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.

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

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.

1.
F.
Bruneval
and
M.
Gatti
, “
Quasiparticle self-consistent GW method for the spectral properties of complex materials
,” in
First Principles Approaches to Spectroscopic Properties of Complex Materials
, edited by
C.
Di Valentin
,
S.
Botti
, and
M.
Cococcioni
(
Springer Berlin Heidelberg
,
Berlin, Heidelberg
,
2014
), pp.
99
135
.
2.
L.
Reining
, “
The GW approximation: Content, successes and limitations
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1344
(
2018
).
3.
X.
Ren
,
P.
Rinke
,
V.
Blum
,
J.
Wieferink
,
A.
Tkatchenko
,
A.
Sanfilippo
,
K.
Reuter
, and
M.
Scheffler
, “
Resolution-of-identity approach to Hartree–Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions
,”
New J. Phys.
14
,
053020
(
2012
).
4.
F.
Bruneval
and
M. A. L.
Marques
, “
Benchmarking the starting points of the GW approximation for molecules
,”
J. Chem. Theory Comput.
9
,
324
329
(
2013
).
5.
E.
Pavarini
,
E.
Koch
,
J.
van den Brink
, and
G.
Sawatzky
,
Quantum Materials: Experiments and Theory, Modeling and Simulation
(
Forschungszentrum Jülich
,
Jülich
,
2016
), Vol. 6, p.
420
.
6.
F.
Bruneval
,
T.
Rangel
,
S. M.
Hamed
,
M.
Shao
,
C.
Yang
, and
J. B.
Neaton
, “
MOLGW 1: Many-body perturbation theory software for atoms, molecules, and clusters
,”
Comput. Phys. Commun.
208
,
149
161
(
2016
).
7.
X.
Blase
,
C.
Attaccalite
, and
V.
Olevano
, “
First-principles GW calculations for fullerenes, porphyrins, phtalocyanine, and other molecules of interest for organic photovoltaic applications
,”
Phys. Rev. B
83
,
115103
(
2011
).
8.
M. J.
van Setten
,
F.
Weigend
, and
F.
Evers
, “
The GW-method for quantum chemistry applications: Theory and implementation
,”
J. Chem. Theory Comput.
9
,
232
246
(
2013
).
9.
J.
Wilhelm
,
M.
Del Ben
, and
J.
Hutter
, “
GW in the Gaussian and plane waves scheme with application to linear acenes
,”
J. Chem. Theory Comput.
12
,
3623
3635
(
2016
).
10.
G. L.
Gutsev
,
B. K.
Rao
, and
P.
Jena
, “
Electronic structure of the 3d metal monoxide anions
,”
J. Phys. Chem. A
104
,
5374
5379
(
2000
).
11.
J. M.
Gonzales
,
R. A.
King
, and
H. F.
Schaefer
, “
Analyses of the ScO and ScO2 photoelectron spectra
,”
J. Chem. Phys.
113
,
567
572
(
2000
).
12.
K. N.
Kudin
and
G. E.
Scuseria
, “
Converging self-consistent field equations in quantum chemistry – recent achievements and remaining challenges
,”
ESAIM: Math. Modell. Numer. Anal.
41
,
281
296
(
2007
).
13.
M. J.
van Setten
,
F.
Caruso
,
S.
Sharifzadeh
,
X.
Ren
,
M.
Scheffler
,
F.
Liu
,
J.
Lischner
,
L.
Lin
,
J. R.
Deslippe
,
S. G.
Louie
,
C.
Yang
,
F.
Weigend
,
J. B.
Neaton
,
F.
Evers
, and
P.
Rinke
, “
GW100: Benchmarking G0W0 for molecular systems
,”
J. Chem. Theory Comput.
11
,
5665
5687
(
2015
).
14.
F.
Caruso
,
M.
Dauth
,
M. J.
van Setten
, and
P.
Rinke
, “
Benchmark of GW approaches for the GW100 test set
,”
J. Chem. Theory Comput.
12
,
5076
5087
(
2016
).
15.
E.
Maggio
,
P.
Liu
,
M. J.
van Setten
, and
G.
Kresse
, “
GW100: A plane wave perspective for small molecules
,”
J. Chem. Theory Comput.
13
,
635
648
(
2017
).
16.
M.
Govoni
and
G.
Galli
, “
GW100: Comparison of methods and accuracy of results obtained with the WEST code
,”
J. Chem. Theory Comput.
14
,
1895
1909
(
2018
).
17.
M.
Shishkin
and
G.
Kresse
, “
Self-consistent GW calculations for semiconductors and insulators
,”
Phys. Rev. B
75
,
235102
(
2007
).
18.
F.
Fuchs
,
J.
Furthmüller
,
F.
Bechstedt
,
M.
Shishkin
, and
G.
Kresse
, “
Quasiparticle band structure based on a generalized Kohn-Sham scheme
,”
Phys. Rev. B
76
,
115109
(
2007
).
19.
J.
Klimeš
,
M.
Kaltak
, and
G.
Kresse
, “
Predictive GW calculations using plane waves and pseudopotentials
,”
Phys. Rev. B
90
,
075125
(
2014
).
20.
M.
Grumet
,
P.
Liu
,
M.
Kaltak
,
J.
Klimeš
, and
G.
Kresse
, “
Beyond the quasiparticle approximation: Fully self-consistent GW calculations
,”
Phys. Rev. B
98
,
155143
(
2018
).
21.
S.
Körbel
,
P.
Boulanger
,
I.
Duchemin
,
X.
Blase
,
M. A. L.
Marques
, and
S.
Botti
, “
Benchmark many-body GW and Bethe–Salpeter calculations for small transition metal molecules
,”
J. Chem. Theory Comput.
10
,
3934
3943
(
2014
).
22.
F.
Kaplan
,
M. E.
Harding
,
C.
Seiler
,
F.
Weigend
,
F.
Evers
, and
M. J.
van Setten
, “
Quasi-particle self-consistent GW for molecules
,”
J. Chem. Theory Comput.
12
,
2528
2541
(
2016
).
23.
C.
Rostgaard
,
K. W.
Jacobsen
, and
K. S.
Thygesen
, “
Fully self-consistent GW calculations for molecules
,”
Phys. Rev. B
81
,
085103
(
2010
).
24.
H.
Wu
and
L.-S.
Wang
, “
Photoelectron spectroscopy and electronic structure of ScOn(n = 1–4) and YOn(n = 1–5): Strong electron correlation effects in ScOand YO
,”
J. Phys. Chem. A
102
,
9129
9135
(
1998
).
25.
H.
Wu
and
L.
Wang
, “
Electronic structure of titanium oxide clusters: TiOy (y = 1–3) and (TiO2)n (n = 1–4)
,”
J. Chem. Phys.
107
,
8221
8228
(
1997
).
26.
H.
Wu
,
S. R.
Desai
, and
L.-S.
Wang
, “
Chemical bonding between Cu and oxygen–copper oxides vs O2 complexes: A study of CuOx (x = 0–6) species by anion photoelectron spectroscopy
,”
J. Phys. Chem. A
101
,
2103
2111
(
1997
).
27.
V. D.
Moravec
,
S. A.
Klopcic
,
B.
Chatterjee
, and
C. C.
Jarrold
, “
The electronic structure of ZnO and ZnF determined by anion photoelectron spectroscopy
,”
Chem. Phys. Lett.
341
,
313
318
(
2001
).
28.
F.
Bruneval
,
N.
Vast
,
L.
Reining
,
M.
Izquierdo
,
F.
Sirotti
, and
N.
Barrett
, “
Exchange and correlation effects in electronic excitations of Cu2O
,”
Phys. Rev. Lett.
97
,
267601
(
2006
).
29.
L. Y.
Isseroff
and
E. A.
Carter
, “
Importance of reference Hamiltonians containing exact exchange for accurate one-shot GW calculations of Cu2O
,”
Phys. Rev. B
85
,
235142
(
2012
).
30.
M. L.
Tiago
and
J. R.
Chelikowsky
, “
Optical excitations in organic molecules, clusters, and defects studied by first-principles Green’s function methods
,”
Phys. Rev. B
73
,
205334
(
2006
).
31.
A.
Seidl
,
A.
Görling
,
P.
Vogl
,
J. A.
Majewski
, and
M.
Levy
, “
Generalized Kohn-Sham schemes and the band-gap problem
,”
Phys. Rev. B
53
,
3764
3774
(
1996
).
32.
T.
Sander
,
E.
Maggio
, and
G.
Kresse
, “
Beyond the Tamm-Dancoff approximation for extended systems using exact diagonalization
,”
Phys. Rev. B
92
,
045209
(
2015
).
33.
Y.-M.
Byun
and
C. A.
Ullrich
, “
Excitons in solids from time-dependent density-functional theory: Assessing the Tamm-Dancoff approximation
,”
Computation
5
,
9
(
2017
).
34.
D.
Golze
,
J.
Wilhelm
,
M. J.
van Setten
, and
P.
Rinke
, “
Core-level binding energies from GW: An efficient full-frequency approach within a localized basis
,”
J. Chem. Theory Comput.
14
,
4856
4869
(
2018
).
35.
A. J.
Cohen
,
D. J.
Tozer
, and
N. C.
Handy
, “
Evaluation of s2 in density functional theory
,”
J. Chem. Phys.
126
,
214104
(
2007
).
36.
A. S.
Menon
and
L.
Radom
, “
Consequences of spin contamination in unrestricted calculations on open-shell species: Effect of Hartree–Fock and Møller–Plesset contributions in hybrid and double-hybrid density functional theory approaches
,”
J. Phys. Chem. A
112
,
13225
13230
(
2008
).
37.
J. G.
Hill
and
J. A.
Platts
, “
Auxiliary basis sets for density fitting–MP2 calculations: Nonrelativistic triple-ζ all-electron correlation consistent basis sets for the 3d elements Sc–Zn
,”
J. Chem. Phys.
128
,
044104
(
2008
).
38.
J. G.
Hill
and
K. A.
Peterson
, “
Explicitly correlated coupled cluster calculations for molecules containing group 11 (Cu, Ag, Au) and 12 (Zn, Cd, Hg) elements: Optimized complementary auxiliary basis sets for valence and core–valence basis sets
,”
J. Chem. Theory Comput.
8
,
518
526
(
2012
).
39.
M. S.
Hybertsen
and
S. G.
Louie
, “
Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies
,”
Phys. Rev. B
34
,
5390
5413
(
1986
).
40.
M.
Véril
,
P.
Romaniello
,
J. A.
Berger
, and
P.-F.
Loos
, “
Unphysical discontinuities in GW methods
,”
J. Chem. Theory Comput.
14
,
5220
5228
(
2018
.
41.
X.
Blase
,
P.
Boulanger
,
F.
Bruneval
,
M.
Fernandez-Serra
, and
I.
Duchemin
, “
GW and Bethe-Salpeter study of small water clusters
,”
J. Chem. Phys.
144
,
034109
(
2016
).
42.
C. L.
Collins
,
K. G.
Dyall
, and
H. F.
Schaefer
, “
Relativistic and correlation effects in CuH, AgH, and AuH: Comparison of various relativistic methods
,”
J. Chem. Phys.
102
,
2024
2031
(
1995
).
43.
K. A.
Peterson
and
C.
Puzzarini
, “
Systematically convergent basis sets for transition metals. II. Pseudopotential-based correlation consistent basis sets for the group 11 (Cu, Ag, Au) and 12 (Zn, Cd, Hg) elements
,”
Theor. Chem. Acc.
114
,
283
296
(
2005
).
44.
P.
Pyykkö
and
J. P.
Desclaux
, “
Relativity and the periodic system of elements
,”
Acc. Chem. Res.
12
,
276
281
(
1979
).
45.
P.
Pyykkö
, “
Relativistic effects in structural chemistry
,”
Chem. Rev.
88
,
563
594
(
1988
).
46.
P.
Pyykkö
, “
Relativistic effects in chemistry: More common than you thought
,”
Annu. Rev. Phys. Chem.
63
,
45
64
(
2012
).
47.
P. S.
Bagus
,
Y. S.
Lee
, and
K. S.
Pitzer
, “
Effects of relativity and of the lanthanide contraction on the atoms from hafnium to bismuth
,”
Chem. Phys. Lett.
33
,
408
411
(
1975
).
48.
H.
Tatewaki
,
S.
Yamamoto
, and
Y.
Hatano
, “
Relativistic effects in the electronic structure of atoms
,”
ACS Omega
2
,
6072
6080
(
2017
.
49.
C.
Hättig
,
W.
Klopper
,
A.
Köhn
, and
D. P.
Tew
, “
Explicitly correlated electrons in molecules
,”
Chem. Rev.
112
,
4
74
(
2012
).
50.
F.
Hüser
,
T.
Olsen
, and
K. S.
Thygesen
, “
Quasiparticle GW calculations for solids, molecules, and two-dimensional materials
,”
Phys. Rev. B
87
,
235132
(
2013
).
51.
T.
Rangel
,
S. M.
Hamed
,
F.
Bruneval
, and
J. B.
Neaton
, “
Evaluating the GW approximation with CCSD(T) for charged excitations across the oligoacenes
,”
J. Chem. Theory Comput.
12
,
2834
2842
(
2016
).
52.
L.
Hung
,
F.
Bruneval
,
K.
Baishya
, and
S.
Öğüt
, “
Benchmarking the GW approximation and Bethe-Salpeter equation for groups Ib and IIb atoms and monoxides
,”
J. Chem. Theory Comput.
13
,
2135
2146
(
2017
).
53.
B.
Shi
,
S.
Weissman
,
F.
Bruneval
,
L.
Kronik
, and
S.
Öğüt
, “
Photoelectron spectra of copper oxide cluster anions from first principles methods
,”
J. Chem. Phys.
149
,
064306
(
2018
).
54.
L.
Hung
,
F.
Bruneval
,
K.
Baishya
, and
S.
Öğüt
, “
Correction to benchmarking the GW approximation and Bethe–Salpeter equation for groups Ib and IIb monoxides
,”
J. Chem. Theory Comput.
13
,
5820
5821
(
2017
).
55.
L.
Hedin
, “
On correlation effects in electron spectroscopies and the GW approximation
,”
J. Phys.: Condens. Matter
11
,
R489
R528
(
1999
).
56.
J. S.
Zhou
,
J. J.
Kas
,
L.
Sponza
,
I.
Reshetnyak
,
M.
Guzzo
,
C.
Giorgetti
,
M.
Gatti
,
F.
Sottile
,
J. J.
Rehr
, and
L.
Reining
, “
Dynamical effects in electron spectroscopy
,”
J. Chem. Phys.
143
,
184109
(
2015
).
57.
F.
Bruneval
, “
Ionization energy of atoms obtained from GW self-energy or from random phase approximation total energies
,”
J. Chem. Phys.
136
,
194107
(
2012
).
58.
P.
Koval
,
D.
Foerster
, and
D.
Sánchez-Portal
, “
Fully self-consistent GW and quasiparticle self-consistent GW for molecules
,”
Phys. Rev. B
89
,
155417
(
2014
).
59.
P.
Liao
and
E. A.
Carter
, “
Testing variations of the GW approximation on strongly correlated transition metal oxides: Hematite (α-Fe2O3) as a benchmark
,”
Phys. Chem. Chem. Phys.
13
,
15189
15199
(
2011
).
60.
K.
Momma
and
F.
Izumi
, “
VESTA: A three-dimensional visualization system for electronic and structural analysis
,”
J. Appl. Crystallogr.
41
,
653
658
(
2008
).
61.
B.
Dai
,
K.
Deng
,
J.
Yang
, and
Q.
Zhu
, “
Excited states of the 3d transition metal monoxides
,”
J. Chem. Phys.
118
,
9608
9613
(
2003
).
62.
A. J.
Bridgeman
and
J.
Rothery
, “
Periodic trends in the diatomic monoxides and monosulfides of the 3d transition metals
,”
J. Chem. Soc., Dalton Trans.
29
,
211
218
(
2000
).
63.
E.
Miliordos
and
A.
Mavridis
, “
Electronic structure and bonding of the early 3d-transition metal diatomic oxides and their ions: ScO0,±, TiO0,±, CrO0,±, and MnO0,±
,”
J. Phys. Chem. A
114
,
8536
8572
(
2010
).
64.
C. W.
Bauschlicher
and
H.
Partridge
, “
A comparison of ZnO and ZnO
,”
J. Chem. Phys.
109
,
8430
8434
(
1998
).
65.
M. B.
Walsh
,
R. A.
King
, and
H. F.
Schaefer
, “
The structures, electron affinities, and energetic stabilities of TiOn and TiOn (n = 1 − 3)
,”
J. Chem. Phys.
110
,
5224
5230
(
1999
).
66.
H.
Xian
,
Z.
Cao
,
X.
Xu
,
X.
Lu
, and
Q.
Zhang
, “
Theoretical study of low-lying electronic states of CuO and CuO
,”
Chem. Phys. Lett.
326
,
485
493
(
2000
).
67.
A.
Daoudi
,
A. T.
Benjelloun
,
J.
Flament
, and
G.
Berthier
, “
Potential energy curves and electronic structure of copper nitrides CuN and CuN+versus CuO and CuO+
,”
J. Mol. Spectrosc.
194
,
8
16
(
1999
).
68.
P. S.
Bagus
,
C. J.
Nelin
, and
C. W.
Bauschlicher
, “
On the low-lying states of CuO
,”
J. Chem. Phys.
79
,
2975
2981
(
1983
).
69.
G. L.
Gutsev
,
L.
Andrews
, and
C. W.
Bauschlicher
, Jr.
, “
Similarities and differences in the structure of 3d-metal monocarbides and monoxides
,”
Theor. Chem. Acc.
109
,
298
308
(
2003
).
70.
C. A.
Fancher
,
H. L.
de Clercq
,
O. C.
Thomas
,
D. W.
Robinson
, and
K. H.
Bowen
, “
Zinc oxide and its anion: A negative ion photoelectron spectroscopic study
,”
J. Chem. Phys.
109
,
8426
8429
(
1998
).
71.
M.
Strange
,
C.
Rostgaard
,
H.
Häkkinen
, and
K. S.
Thygesen
, “
Self-consistent GW calculations of electronic transport in thiol- and amine-linked molecular junctions
,”
Phys. Rev. B
83
,
115108
(
2011
).
72.
N.
Marom
,
F.
Caruso
,
X.
Ren
,
O. T.
Hofmann
,
T.
Körzdörfer
,
J. R.
Chelikowsky
,
A.
Rubio
,
M.
Scheffler
, and
P.
Rinke
, “
Benchmark of GW methods for azabenzenes
,”
Phys. Rev. B
86
,
245127
(
2012
).

Supplementary Material