Practical $GW$ scheme for electronic structure of 3$d$-transition-metal monoxide anions: ScO$^{-}$, TiO$^{-}$, CuO$^{-}$, and ZnO$^{-}$

TL;DR

This study evaluates the performance of the $GW$ approximation, especially the $G_{n}W_{0}$@PBE scheme, in predicting electronic spectra of small 3$d$ transition-metal oxide molecules, demonstrating its effectiveness and potential for complex systems.

Contribution

The paper extends the application of the $GW$ approximation to small molecules with 3$d$ transition metals, showing that $G_{n}W_{0}$@PBE provides accurate results for both localized and delocalized states.

Findings

$G_{n}W_{0}$@PBE yields good agreement with experiments for 3$d$ transition-metal oxides.

The scheme benefits from a cancellation effect between overscreening and underscreening.

It offers a reliable, efficient method for electronic structure prediction in complex systems.

Abstract

The $GW$ approximation to many-body perturbation theory is a reliable tool for describing charged electronic excitations, and it has been successfully applied to a wide range of extended systems for several decades using a plane-wave basis. However, the $GW$ approximation has been used to test limited spectral properties of a limited set of finite systems (e.g. frontier orbital energies of closed-shell $sp$ molecules) only for about a decade using a local-orbital basis. Here, we calculate the quasiparticle spectra of closed- and open-shell molecular anions with partially and completely filled 3$d$ shells (shallow and deep 3$d$ states, respectively), ScO$^{-}$, TiO$^{-}$, CuO$^{-}$, and ZnO$^{-}$, using various levels of $GW$ theory, and compare them to experiments to evaluate the performance of the $GW$ approximation on the electronic structure of small molecules containing 3$d$ transition metals. We find that the $G$-only eigenvalue self-consistent $GW$ scheme with $W$ fixed to the PBE level ($G_{n}W_{0}$@PBE), which gives the best compromise between accuracy and efficiency for solids, also gives good results for both localized ($d$) and delocalized ($sp$) states of 3$d$-transition-metal oxide molecules. The success of $G_{n}W_{0}$@PBE in predicting electronic excitations in these systems reasonably well is likely due to the fortuitous cancellation effect between the overscreening of the Coulomb interaction by PBE and the underscreening by the neglect of vertex corrections. Together with the absence of the self-consistent field convergence error (e.g. spin contamination in open-shell systems) and the $GW$ multi-solution issue, the $G_{n}W_{0}$@PBE scheme gives the possibility to predict the electronic structure of complex real systems (e.g. molecule-solid and $sp$-$d$ hybrid systems) accurately and efficiently.

Tables

Table 1. Table 1: Number of occupied and empty states for the ↑ ↑ \uparrow -spin channel ( N occ ↑ superscript subscript 𝑁 occ ↑ N_{\text{occ}}^{\uparrow} and N emp ↑ superscript subscript 𝑁 emp ↑ N_{\text{emp}}^{\uparrow} , respectively) used in G ​ W 𝐺 𝑊 GW calculations. FC means the frozen-core approximation. AE and ECP mean all electron and effective core potential, respectively. CN means the cardinal number.

			$N_{\text{occ}}^{\uparrow }$	$N_{\text{emp}}^{\uparrow }$
Anion	Potential	FC		CN=2	CN=3	CN=4	CN=5
ScO-	AE	Yes	9	67	124	205	314
ScO-	AE	No	15	67	124	205	314
TiO-	AE	Yes	10	66	123	204	313
TiO-	AE	No	16	66	123	204	313
CuO-	AE	Yes	13	63	120	201	310
CuO-	AE	No	19	63	120	201	310
CuO-	ECP	Yes	13	63	120	201	310
CuO-	ECP	No	14	63	120	201	310
ZnO-	AE	Yes	14	62	119	200	309
ZnO-	AE	No	20	62	119	200	309
ZnO-	ECP	Yes	14	62	119	200	309
ZnO-	ECP	No	15	62	119	200	309

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

	ScO-	TiO-		CuO-	ZnO-
		$\uparrow$	$\downarrow$		$\uparrow$	$\downarrow$
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$\alpha$($\alpha$=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

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

			$G_0{W}_0@\text{PBE}\alpha$					$G_n{W}_0@$	$G_n{W}_n@$
	State	Exp.	$\alpha$=0.00	$\alpha$=0.25	$\alpha$=0.50	$\alpha$=0.75	$\alpha$=1.00	PBE	PBE	Others
ScO- (${}_{}{}^1\Sigma _{}^{+}$)
HOMO	${}_{}{}^1\Sigma _{}^{+}$	1.35111Ref. Wu and Wang, 1998	0.51	1.15	1.45	1.59	1.63	1.21	1.35	1.28777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 1.26101010Ref. Gutsev, Rao, and Jena, 2000 using the B3LYP functional, 1.19111111Ref. Gonzales, King, and Schaefer, 2000 using the B3LYP functional
HOMO-1	${}_{}{}^2\Delta$	3.10111Ref. Wu and Wang, 1998	3.30	4.77	5.45	5.72	5.82	5.34	6.08	2.41777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 2.78888Ref. Bridgeman and Rothery, 2000, 3.31999Ref. Miliordos and Mavridis, 2010 using the multi-reference configuration interaction (MRCI) method
HOMO-2	${}_{}{}^2\Pi$	3.40111Ref. Wu and Wang, 1998	3.42	4.81	5.40	5.61	5.63	5.39	6.17	3.34777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 3.24888Ref. Bridgeman and Rothery, 2000, 3.44999Ref. Miliordos and Mavridis, 2010 using the multi-reference configuration interaction (MRCI) method
MAE			0.84555HOMO-1 and HOMO-2 are not included because our $GW$ calculations cannot account for two-electron transitions (see text).	0.18555HOMO-1 and HOMO-2 are not included because our $GW$ calculations cannot account for two-electron transitions (see text).	0.10555HOMO-1 and HOMO-2 are not included because our $GW$ calculations cannot account for two-electron transitions (see text).	0.24555HOMO-1 and HOMO-2 are not included because our $GW$ calculations cannot account for two-electron transitions (see text).	0.28555HOMO-1 and HOMO-2 are not included because our $GW$ calculations cannot account for two-electron transitions (see text).	0.14555HOMO-1 and HOMO-2 are not included because our $GW$ calculations cannot account for two-electron transitions (see text).	0.00555HOMO-1 and HOMO-2 are not included because our $GW$ calculations cannot account for two-electron transitions (see text).
TiO- (${}_{}{}^2\Delta$)
$\uparrow$-HOMO	${}_{}{}^1\Sigma _{}^{+}$	2.00222Ref. Wu and Wang, 1997	0.26	2.24	3.70	4.79	5.63	1.83	2.74	2.39777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 2.37888Ref. Bridgeman and Rothery, 2000, 2.34999Ref. Miliordos and Mavridis, 2010 using the multi-reference configuration interaction (MRCI) method
$\uparrow$-HOMO-1	${}_{}{}^1\Delta$	1.73222Ref. Wu and Wang, 1997	0.53	1.28	1.65	1.83	1.95	1.38	1.61	1.88777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 1.72999Ref. Miliordos and Mavridis, 2010 using the multi-reference configuration interaction (MRCI) method
$\downarrow$-HOMO	${}_{}{}^3\Delta$	1.30222Ref. Wu and Wang, 1997	0.31	1.00	1.29	1.43	1.55	1.06	1.21	1.19777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 1.18101010Ref. Gutsev, Rao, and Jena, 2000 using the B3LYP functional, 1.14141414Ref. Walsh, King, and Schaefer, 1999 using the B3LYP functional
MAE			1.31	0.33	0.60	1.01	1.37	0.25	0.32
CuO- (${}_{}{}^1\Sigma _{}^{+}$)
HOMO	${}_{}{}^2\Pi$	1.78333Ref. Wu, Desai, and Wang, 1997	0.40	1.40	1.58	1.40	0.97	2.19	3.17	1.55777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 1.52101010Ref. Gutsev, Rao, and Jena, 2000 using the B3LYP functional, 0.46151515Ref. Hung et al., 2017b using the $G_0{W}_0@\text{PBE}$ method
HOMO-1	${}_{}{}^2\Sigma _{}^{+}$	2.75333Ref. Wu, Desai, and Wang, 1997	1.39	2.17	2.23	1.98	1.56	2.66	3.58	2.96777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 2.86888Ref. Bridgeman and Rothery, 2000, 1.60151515Ref. Hung et al., 2017b using the $G_0{W}_0@\text{PBE}$ method, 2.78161616Ref. Xian et al., 2000 using the CCSD(T) method, 2.47171717Ref. Daoudi et al., 1999, 2.81181818Ref. Bagus, Nelin, and Bauschlicher, 1983 using the single and double excitation configuration interaction (SDCI) method
HOMO-2		4.50666We chose this value from the Z band in the PES spectrum of CuO- (see text).	2.18	4.05	4.60	4.56	4.25	4.88	6.38	4.07777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 4.01888Ref. Bridgeman and Rothery, 2000, 4.58161616Ref. Xian et al., 2000 using the CCSD(T) method, 4.50171717Ref. Daoudi et al., 1999
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\Sigma _{}^{+}$)
$\uparrow$-HOMO	${}_{}{}^1\Sigma _{}^{+}$	2.09444Ref. Moravec et al., 2001	0.91	1.91	2.20	2.23	2.00	2.11	2.57	2.19777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 2.33101010Ref. Gutsev, Rao, and Jena, 2000 using the B3LYP functional, 2.29121212Ref. Bauschlicher and Partridge, 1998 using the B3LYP functional, 2.10131313Ref. Körbel et al., 2014 using the $G_0{W}_0@\text{PBE}$0 method, 1.06151515Ref. Hung et al., 2017b using the $G_0{W}_0@\text{PBE}$ method
$\uparrow$-HOMO-1	${}_{}{}^1\Pi$	2.71444Ref. Moravec et al., 2001	1.36	2.43	3.03	3.72	4.79	2.77	3.73	2.62777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 1.43151515Ref. Hung et al., 2017b using the $G_0{W}_0@\text{PBE}$ method
$\uparrow$-HOMO-2			3.24	4.04	4.89	5.64	6.90	4.77	5.66	3.50151515Ref. Hung et al., 2017b using the $G_0{W}_0@\text{PBE}$ method
$\downarrow$-HOMO	${}_{}{}^3\Pi$	2.40444Ref. Moravec et al., 2001	1.17	2.17	2.64	3.15	4.02	2.65	3.48	2.41777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 1.20151515Ref. Hung et al., 2017b using the $G_0{W}_0@\text{PBE}$ method
$\downarrow$-HOMO-1	${}_{}{}^3\Sigma _{}^{+}$	3.96444Ref. Moravec et al., 2001	2.71	3.35	3.47	3.31	3.00	4.11	4.51	4.15777Ref. Dai et al., 2003 using 6-3111+G* basis sets, 2.89151515Ref. Hung et al., 2017b using the $G_0{W}_0@\text{PBE}$ method
MAE			1.25	0.33	0.29	0.64	1.19	0.12	0.78

Table 4. Table 4: Optimal amount of EXX in the G 0 ​ W 0 subscript 𝐺 0 subscript 𝑊 0 G_{0}W_{0} starting point for gas-phase small molecules (highlighted in bold).

Reference	Körbel et al. Körbel et al. (2014)	Bruneval et al. Bruneval and Marques (2013)	Kaplan et al. Kaplan et al. (2016)	Rostgaard et al. Rostgaard, Jacobsen, and Thygesen (2010)	This work
Code	FIESTA	MOLGW	TURBOMOLE	GPAW	MOLGW
Optimal EXX	25%	50%	75%	100%	25–50%
Tested EXX	25 & 100%	0, 20, 25 & 50%	0, 25 & 75%	0 & 100%	0, 25, 50, 75 & 100%
System	39 closed-shell 3,4,5$d$	34 closed-shell $sp$111The same set of molecules is used.	29 closed-shell $sp$	34 closed-shell $sp$111The same set of molecules is used.	4 closed- & open-shell 3$d$
	& 9 closed-shell $sp$
System size	2–7 atoms	2–8 atoms	2–18 atoms	2–8 atoms	2 atoms
Property	HOMO & LUMO	HOMO	HOMO &	HOMO	HOMO-$n$ ($n$ = 0, 1, …)
			HOMO-$n$ ($n$ = 0, 1, …)222For naphthalene only		(focusing on 3$d$ MO)
Reference data	Experiment	$\Delta$SCF333Ref. Bruneval and Marques, 2013 showed that $\Delta$SCF using CCSD(T) with CN=4 causes an error of $\sim$0.1 eV in the IE of small $sp$ molecules with respect to experiment (the largest being 0.67 eV for NaCl).	QS$GW$444Ref. Bruneval, 2012 showed that QS$GW$ with CN=5 causes a mean absolute error of 0.18 eV in the IE of the first row atoms with respect to experiment (the largest being $\sim$0.4 eV for O).	Experiment	Experiment
$\omega$ integration	Contour deformation555Refs. van Setten et al., 2015; Golze et al., 2018 showed that the contour deformation technique produces almost the same $GW$ self-energy as the fully analytic method for frontier and non-frontier orbitals, respectively.	Fully analytic	Fully analytic	Fully analytic	Fully analytic
QP equation	Linearization	Linearization	Spectral function	Linearization666For 0% EXX &	Graphical solution &
				Spectral function777For 100% EXX	Spectral function
$\eta$	Not available	Not available	0.001 eV	0 eV	0.002 or 0.005 Ha
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	PAW888Projector-Augmented Wave	AE
RI	Used	Not used	Used	Not applicable999GPAW uses augmented Wannier basis sets, whereas FIESTA, MOLGW, and TURBOMOLE use Gaussian basis sets.	Not used

Equations

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

$$G^{\sigma }(\mathbf{r},{\mathbf{r}}^{\prime },\omega )=\sum _i\frac{{\phi }_i^{\sigma }(\mathbf{r}){\phi }_i^{\sigma }({\mathbf{r}}^{\prime })}{\omega -{ϵ}_i^{\sigma }-i\eta }+\sum _a\frac{{\phi }_a^{\sigma }(\mathbf{r}){\phi }_a^{\sigma }({\mathbf{r}}^{\prime })}{\omega -{ϵ}_a^{\sigma }+i\eta },$$

$${\chi }_0=-i\sum _{\sigma }{G}^{\sigma }{G}^{\sigma },$$

$$W=v+v{\chi }_{0}W=v+v{\chi }_{0}v+v{\chi }_{0}v{\chi }_{0}v+...=v+v\chi v,$$

$${\Sigma }^{\sigma }=i{G}^{\sigma }W=i{G}^{\sigma }(v+v\chi v)={\Sigma }_{\text{x}}^{\sigma }+{\Sigma }_{\text{c}}^{\sigma },$$

$${ϵ}_m^{{\text{G}}_0{\text{W}}_0,\sigma }={ϵ}_m^{\sigma }+⟨{\phi }_m^{\sigma }|\text{Re}{\Sigma }^{\sigma }({ϵ}_m^{{\text{G}}_0{\text{W}}_0,\sigma })-v_{\text{xc}}^{\sigma }|{\phi }_m^{\sigma }⟩,$$

$$A_{mm}^{\sigma }(\mathbf{r},{\mathbf{r}}^{\prime },\omega )=\frac{1}{\pi }|\text{Im}{G}_{mm}^{\sigma }(\mathbf{r},{\mathbf{r}}^{\prime },\omega )|,$$

$$\omega ={ϵ}_m^{\sigma }+⟨{\phi }_m^{\sigma }|\text{Re}{\Sigma }^{\sigma }(\omega )-v_{\text{xc}}^{\sigma }|{\phi }_m^{\sigma }⟩,$$

$${\phi }_m^{\sigma }(\mathbf{r})=\sum _{\mu }{C}_{\mu m}^{\sigma }{\varphi }_{\mu }(\mathbf{r}),$$

$${\mathbf{H}}^{\sigma }{\mathbf{C}}^{\sigma }=\mathbf{S}{\mathbf{C}}^{\sigma }{ϵ}^{\sigma },$$

$$S_{\mu \nu }=\int d\mathbf{r}{\varphi }_{\mu }(\mathbf{r}){\varphi }_{\nu }(\mathbf{r}),$$

$$H_{\mu \nu }^{\sigma }=T_{\mu \nu }+V_{\text{ext},\mu \nu }+J_{\mu \nu }-\alpha K_{\mu \nu }^{\sigma }+(1-\alpha )V_{\text{x},\mu \nu }^{\text{PBE},\sigma }+V_{\text{c},\mu \nu }^{\text{PBE},\sigma },$$

$$J_{\mu \nu }=\sum _{\lambda \tau }(\mu \nu |\lambda \tau )\sum _{\sigma }{D}_{\lambda \tau }^{\sigma },$$

$$(\mu \nu |\lambda \tau )=\iint d\mathbf{r}d{\mathbf{r}}^{\prime }{\varphi }_{\mu }(\mathbf{r}){\varphi }_{\nu }(\mathbf{r})\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\varphi }_{\lambda }({\mathbf{r}}^{\prime }){\varphi }_{\tau }({\mathbf{r}}^{\prime }),$$

$$D_{\mu \nu }^{\sigma }=\sum _m{f}_m^{\sigma }{C}_{\mu m}^{\sigma }{C}_{\nu m}^{\sigma },$$

$$K_{\mu \nu }^{\sigma }=\sum _{\lambda \tau }{D}_{\lambda \tau }^{\sigma }(\mu \lambda |\tau \nu ).$$

$${\rho }^{\sigma }(\mathbf{r})=\sum _{\mu \nu }{D}_{\mu \nu }^{\sigma }{\varphi }_{\mu }(\mathbf{r}){\varphi }_{\nu }(\mathbf{r}).$$

$$\left(\begin{array}{cc}[r]\mathbf{A}& \mathbf{B}\\ -\mathbf{B}& -\mathbf{A}\end{array}\right)\left(\begin{array}{c}{X}^s\\ Y^s\end{array}\right)={\Omega }_s\left(\begin{array}{c}{X}^s\\ Y^s\end{array}\right),$$

$$A_{ia\sigma }^{jb{\sigma }^{\prime }}$$

$$B_{ia\sigma }^{jb{\sigma }^{\prime }}$$

$$(ia\sigma |jb{\sigma }^{\prime })=\iint d\mathbf{r}d{\mathbf{r}}^{\prime }{\phi }_i^{\sigma }(\mathbf{r}){\phi }_a^{\sigma }(\mathbf{r})\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\phi }_j^{{\sigma }^{\prime }}({\mathbf{r}}^{\prime }){\phi }_b^{{\sigma }^{\prime }}({\mathbf{r}}^{\prime }).$$

$$(ia\sigma |jb{\sigma }^{\prime })=\sum _{\mu \nu \lambda \tau }{C}_{\mu i}^{\sigma }{C}_{\nu a}^{\sigma }{C}_{\lambda j}^{{\sigma }^{\prime }}{C}_{\tau b}^{{\sigma }^{\prime }}(\mu \nu |\lambda \tau ),$$

$$⟨{\phi }_m^{\sigma }|{\Sigma }_{\text{x}}^{\sigma }|{\phi }_m^{\sigma }⟩$$

$$⟨{\phi }_m^{\sigma }|{\Sigma }_{\text{c}}^{\sigma }(\omega )|{\phi }_m^{\sigma }⟩$$

$$+\sum _{as}\frac{w_{ma\sigma }^s{w}_{ma\sigma }^s}{\omega -{ϵ}_a^{\sigma }-{\Omega }_s+i\eta },$$

$$w_{mn\sigma }^s=\sum _{ia{\sigma }^{\prime }}(mn\sigma |ia{\sigma }^{\prime })(X_{ia{\sigma }^{\prime }}^s+Y_{ia{\sigma }^{\prime }}^s).$$

$$⟨S^2⟩=S(S+1)+N_{\downarrow }-\sum _i^{N_{\uparrow }}\sum _j^{N_{\downarrow }}|⟨{\phi }_i^{\uparrow }|{\phi }_j^{\downarrow }⟩{|}^2,$$

$$(\mu \nu |\lambda \tau )\approx \sum _{PQ}(\mu \nu |P)(P|Q{)}^{-1}(Q|\lambda \tau ),$$

$${ϵ}_m^{{\text{G}}_0{\text{W}}_0}\approx {ϵ}_m+Z_m^{\text{linear}}⟨{\phi }_m|\Sigma ({ϵ}_m)-v_{\text{xc}}|{\phi }_m⟩\equiv {ϵ}_m^{{\text{G}}_0{\text{W}}_0,\text{linear}},$$

$$Z_m^{\text{linear}}({ϵ}_m)=\frac{1}{1-\frac{\partial }{\partial \omega }⟨{\phi }_m|\Sigma (\omega )|{\phi }_m⟩{|}_{\omega ={ϵ}_m}},$$

$$Z_m({ϵ}_m^{{\text{G}}_0{\text{W}}_0})=\frac{1}{1-\frac{\partial }{\partial \omega }⟨{\phi }_m|\Sigma (\omega )|{\phi }_m⟩{|}_{\omega ={ϵ}_m^{{\text{G}}_0{\text{W}}_0}}}.$$

$${ϵ}_m^{\text{evGW},1}$$

Full text

I Introduction

It is a challenging task to accurately determine the electronic structure of an interacting many-electron system. In experiment, electron removal and addition energies of both extended and finite systems are measured by direct and inverse photoelectron spectroscopy (PES and IPES, respectively). In theory, it is well known that the $GW$ approximation to many-body perturbation theory (MBPT) describes bandgaps and band structures of solids more accurately than local and semi-local approximations to density-functional theory (DFT). Bruneval and Gatti (2014); Reining (2018) However, less is known about the performance of the $GW$ approximation on the electronic structure of atoms, molecules, and clusters. Especially, $GW$ calculations for the quasiparticle (QP) spectra of open-shell molecules containing 3$d$ transition metals are scarce. There are a few reasons for it.

First, it is easier to test only frontier orbital energies such as the ionization energy (IE) and the electron affinity (EA) than the full QP spectrum (all orbital energies). There are mainly two ways to calculate IE and EA of molecules. On 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 $\Delta$SCF (self-consistent field) method. Ren et al. (2012) Generally, the $\Delta$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 $\Delta$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 $\Delta$SCF method, many studies have evaluated the performance of the $GW$ approximation on molecules by comparing $GW$ IE and EA to $\Delta$SCF ones, Bruneval and Marques (2013) 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. Pavarini et al. (2016)

Second, it is easier to test closed-shell systems than open-shell ones. Most of quantum chemistry-based $GW$ implementations for finite systems, such as MOLGW, Bruneval et al. (2016) FIESTA, Blase, Attaccalite, and Olevano (2011) TURBOMOLE, van Setten, Weigend, and Evers (2013) FHI-AIMS, Ren et al. (2012) and CP2K, Wilhelm, Del Ben, and Hutter (2016) use local-orbital basis sets such as Gaussian basis sets. $GW$ calculations require mean-field self-consistent calculations, such as restricted and unresticted Hartree-Fock or Kohn-Sham (RHF or RKS and UHF or UKS, respectively) calculations for closed- and open-shell systems, respectively. The problem is that unlike RHF and RKS self-consistent calculations, UHF and UKS ones are not guaranteed to converge, as their convergence strongly depends on the initial guess wavefunctions. This is especially the case for spin-unrestricted calculations performed with hybrid exchange-correlation (xc) functionals, which include a fraction of exact exchange (EXX), and HF on 3$d$-transition-metal-containing molecules due to the near-degeneracy of energy levels. Gutsev, Rao, and Jena (2000); Gonzales, King, and Schaefer (2000); Kudin, Konstantin N. and Scuseria, Gustavo E. (2007) 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. van Setten et al., 2015; Caruso et al., 2016; Maggio et al., 2017; Govoni and Galli, 2018 used the so-called $GW$100 benchmark set, which is composed of only closed-shell molecules.

Last, it is easier to test $sp$-electron systems than $d$-electron ones. Fundamentally, it is more difficult to accurately predict the electronic structure of $d$ systems (especially, 3$d$ systems) than $sp$ ones because of the strong localization, and thus the strong correlation, of $d$ electrons. For example, it is challenging for $GW$ to accurately reproduce the experimental bandgap and $d$-band position of bulk ZnO at the same time. Shishkin and Kresse (2007); Fuchs et al. (2007); Klimeš, Kaltak, and Kresse (2014) 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. Gutsev, Rao, and Jena (2000); Gonzales, King, and Schaefer (2000); Kudin, Konstantin N. and Scuseria, Gustavo E. (2007)

The $GW$ approximation is unique, but due to its high computational costs, there are various $GW$ schemes and variants. Generally, there are two approaches. One approach is to vary the self-consistency level in the $GW$ approximation. The $GW$ self-consistent levels from the lowest to the highest include the perturbative non-self-consistent (one-shot) $GW$ ($G_{0}W_{0}$) scheme, the eigenvalue self-consistent $GW$ (ev$GW$) scheme (with two types $G_{n}W_{0}$ and $G_{n}W_{n}$, which update eigenvalues only in $G$ and in both $G$ and $W$, respectively), the QP self-consistent $GW$ (QS$GW$) scheme using a static and Hermitian approximation to the $GW$ self-energy, and the fully self-consistent $GW$ (SC$GW$) scheme. Bruneval and Gatti (2014) 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 particle number, momentum, and energy. However, the higher $GW$ self-consistency level does not necessarily give more accurate QP energies because vertex corrections are missing in the $GW$ approximation. For example, SC$GW$ and QS$GW$ systematically overestimate the bandgaps of solids, Grumet et al. (2018) displaying worse performance than ev$GW$, which currently provides the best balance between accuracy and efficiency for solids. Shishkin and Kresse (2007)

The other approach is to vary the amount of EXX in the $GW$ starting point to reduce the self-interaction error by (semi-)local xc functionals. Typically, the $G_{0}W_{0}$ scheme chooses this approach to obtain good results at low computational costs. However, the predictive power of this approach is questionable, since the optimal amount of EXX in the $GW$ starting point is strongly system-dependent. For example, for extended systems, the reported values for the optimal amount of EXX are narrowly spread between 0% and 25%, Fuchs et al. (2007), while for finite systems, they are widely spread between 25% and 100%. Körbel et al. (2014); Bruneval and Marques (2013); Kaplan et al. (2016); Rostgaard, Jacobsen, and Thygesen (2010)

The purpose of this work is to evaluate the performance of the $GW$ approximation on the electronic structure of small oxide molecules containing 3$d$ transition metals. To this end, we calculate the QP spectra of closed- and open-shell molecular anions with partially and completely filled 3$d$ shells, ScO*-, TiO-, CuO-, and ZnO-, using various levels of $GW$ theory. There are a few reasons why we chose these molecular systems: (i) their anion PES data is available, Wu and Wang (1998, 1997); Wu, Desai, and Wang (1997); Moravec et al. (2001) (ii) CuO-* and ZnO*-* are molecular analogs to bulk Cu2O and ZnO, respectively, which are challenging systems for the $GW$ method, Bruneval et al. (2006); Isseroff and Carter (2012) and (iii) shallow and deep 3$d$ states are measured in TiO*-* and CuO*-*, respectively.

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

II Theoretical Background

In this section, we briefly review the $GW$ approximation and its implementation using local-orbital basis sets. This section contains only a minimal number of equations, which will be needed later. More details can be found in Refs. Bruneval et al., 2016; van Setten et al., 2015; van Setten, Weigend, and Evers, 2013; Tiago and Chelikowsky, 2006. Generally, we follow the notation in the MOLGW implementation paper Bruneval et al. (2016) 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 refered to as the polarizability instead of the susceptibility, and (v) $\chi$ is used for the polarizability instead of $P$ and $\Pi$.

II.1 $GW$ Approximation

In Hedin’s $GW$ approximation, the non-local, dynamical, and non-Hermitian self-energy $\Sigma^{\sigma}$ at frequency $\omega$ is given by

[TABLE]

where $\sigma$ is the spin channel ($\uparrow$ or $\downarrow$), $G^{\sigma}$ is the time-ordered one-particle Green’s function, $W$ is the dynamically screened Coulomb interaction, and $\eta$ 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 $\epsilon^{\sigma}_{m}$ and corresponding wavefunctions $\varphi^{\sigma}_{m}(\mathbf{r})$ obtained from the Hartree or mean-field approximation:

[TABLE]

where $i$ runs over occupied states and $a$ runs over empty states. Note that $G^{\sigma}$ in Eq. (2) is not the interacting (dressed) Green’s function, but the non-interacting (bare) one, which are conventionally denoted by $G^{\sigma}$ and $G^{\sigma}_{0}$, 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 non-interacting (irreducible) polarizability $\chi_{0}$ and the interacting (reducible) polarizability $\chi=\chi_{0}[1-v\chi_{0}]^{-1}$ within the RPA, the screened Coulomb interaction, and the self-energy:

[TABLE]

[TABLE]

[TABLE]

where $v$ denotes the bare (unscreened) Coulomb interaction $v(\mathbf{r},\mathbf{r^{\prime}})=1/|\mathbf{r}-\mathbf{r^{\prime}}|$, $\Sigma_{\text{x}}$ is the exchange part of the self-energy, and $\Sigma_{\text{c}}$ is the correlation part of the self-energy. Note that in Eqs. (3), (4), and (5), space and frequency variables ($\mathbf{r}$, $\mathbf{r^{\prime}}$, $\omega$) are omitted for simplicity, and $\chi_{0}(\omega)$, $\chi(\omega)$, $\Sigma^{\sigma}(\omega)$, and $\Sigma^{\sigma}_{\text{c}}(\omega)$ are dynamic, whereas $v$ and $\Sigma^{\sigma}_{\text{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:

[TABLE]

where $\epsilon^{\text{G}_{0}\text{W}_{0},\sigma}_{m}$ are the perturbative one-shot $GW$ QP energies and $v^{\sigma}_{\text{xc}}$ is the xc potential. Experimentally, $\epsilon^{\text{G}_{0}\text{W}_{0},\sigma}_{m}$ correspond to vertical IEs and EAs in PES and IPES, respectively. Theoretically, $\epsilon^{\text{G}_{0}\text{W}_{0},\sigma}_{m}$ correspond to the positions of poles of the 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^{\sigma}$:

[TABLE]

where $A^{\sigma}_{mm}$ are the diagonal elements of the spectral function, $G^{\sigma}_{mm}$ are the diagonal elements of the Green’s function, and Im represents the imaginary part. Note that $A^{\sigma}_{mm}$ give the local density of states.

$G^{\sigma}$ in Eq. (7) is the interacting Green’s function, whereas $G^{\sigma}$ in Eq. (2) is the non-interacting one. In other words, by plugging $G^{\sigma}$ in Eq. (2) into Eq. (7) after replacing $\epsilon^{\sigma}_{m}$, where $m$ = $i$ or $a$, by $\epsilon^{\sigma}_{m}+\langle\varphi^{\sigma}_{m}|\Sigma^{\sigma}(\omega)-v^{\sigma}_{\text{xc}}|\varphi^{\sigma}_{m}\rangle$, one can find that $A^{\sigma}_{mm}$ have Lorentzian peaks at

[TABLE]

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 non-linear, because $\Sigma^{\sigma}$ depends on $\epsilon^{\text{G}_{0}\text{W}_{0},\sigma}_{m}$, so it should be solved numerically. Additionally, Hedin equations are coupled, because $W$ and $\Sigma^{\sigma}$ depend on $G^{\sigma}$, so they should be solved self-consistently. Multiple ways to numerically solve the non-linear QP equation and to iteratively solve the coupled Hedin equations will be discussed later.

II.2 Self-Consistent Field Method

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

[TABLE]

where $C^{\sigma}_{\mu 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 scheme Seidl et al. (1996) for (semi-)local functionals, hybrid functionals, and HF]:

[TABLE]

where $\mathbf{C}^{\sigma}$ is a matrix of MO coefficients, $\mathbf{\epsilon}^{\sigma}$ is a diagonal matrix of MO energies, $\mathbf{S}$ is the AO overlap matrix with elements:

[TABLE]

and $\mathbf{H}^{\sigma}$ is the Hamiltonian matrix with elements:

[TABLE]

where $T$, $V_{\text{ext}}$, $J$, and $K^{\sigma}$ are the kinetic energy, external potential energy, Hartree, and Fock exchange terms, respectively, $V^{\sigma}_{\text{x}}$ and $V^{\sigma}_{\text{c}}$ are the exchange and correlation potentials, respectively, and $\alpha$ 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

[TABLE]

where $(\mu\nu|\lambda\tau)$ are the 4-center two-electron Coulomb repulsion integrals:

[TABLE]

and $\mathbf{D}^{\sigma}$ is the density matrix with elements:

[TABLE]

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

[TABLE]

The exchange and correlation potentials in Eq. (12) depend on the density $\rho^{\sigma}$ (and the density gradient $\nabla\rho^{\sigma}$):

[TABLE]

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^{\sigma}$, $V^{\sigma}_{\text{x}}$, and $V^{\sigma}_{\text{c}}$ in Eq. (12) depend on the density matrix in Eq. (15), as shown in Eqs. (13), (16), and (17).

II.3 $GW$ Self-Energy

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

[TABLE]

where $\mathbf{A}$ and $\mathbf{B}$ are the resonant and coupling matrices, respectively, and $\Omega_{s}$ and $(X^{s},Y^{s})$ are the eivenvalues (the neutral two-particle excitation energies) and corresponding eigenvectors, respectively. The matrix elements in $\mathbf{A}$ and $\mathbf{B}$ are given by

[TABLE]

where $i$ and $j$ are for occupied states, $a$ and $b$ are for empty states, $f_{\text{xc}}$ is the time-dependent density-functional theory (TDDFT) xc kernel, and $(ia\sigma|jb\sigma^{\prime})$ are the 4-orbital two-electron Coulomb repulsion integrals:

[TABLE]

In this work, we used the RPA by setting $f_{\text{xc}}=0$. Note that Ref. van Setten, Weigend, and Evers, 2013 showed that TDDFT and RPA polarizabilities make a difference of $\sim$0.1 eV in $G_{0}W_{0}@\text{PBE}$ QP HOMO energy. Note also that the dimension of the Casida matrix in Eq. (18) scales as $O(N^{2})$ with $N$ being the system size, so building and completely diagonalizing the Casida matrix scale as $O(N^{4})$ and $O(N^{6})$, respectively. The MO integrals in Eq. (21) are transformed from the AO integrals in Eq. (14) through the AO-MO integral transformation:

[TABLE]

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

Diagonalizing the Casida matrix in Eq. (18) yields eigenvalues $\Omega_{s}$ and eigenvectors $(X^{s},Y^{s})$. In MOLGW, the diagonalization is performed without using the Tamm-Dancoff approximation (TDA), which sets $\mathbf{B}$ to zero, but efficiently using the so-called beyond-TDA method. Sander, Maggio, and Kresse (2015); Bruneval et al. (2016); Byun and Ullrich (2017). Using $\Omega_{s}$ and $(X^{s},Y^{s})$, one can construct the spectral representation of the interacting polarizability $\chi({\omega})$. Tiago and Chelikowsky (2006); van Setten et al. (2015); van Setten, Weigend, and Evers (2013) From $\chi({\omega})$ and Eq. (4), one can obtain the spectral representation of the screened Coulomb interaction $W({\omega})$. Tiago and Chelikowsky (2006); van Setten, Weigend, and Evers (2013); van Setten et al. (2015); Bruneval et al. (2016) Using $W({\omega})$ and Eq. (5), and analytically performing the convolution of $G^{\sigma}(\omega)$ and $W(\omega)$ in the frequency domain, one can obtain the exchange and correlation parts of the $GW$ self-energy $\Sigma^{\sigma}_{\text{x}}$ and $\Sigma^{\sigma}_{\text{c}}(\omega)$, respectively, whose diagonal matrix elements are given by

[TABLE]

where $i$ runs over occupied states, $a$ runs over empty states, $s$ runs over all excitations, and $w^{s}_{mn\sigma}$ are given by

[TABLE]

Note that unlike the plasmon-pole approximation (PPA), the analytic continuation method, and the contour deformation technique, van Setten et al. (2015); Golze et al. (2018) the fully analytic method employed in RGWBS, Tiago and Chelikowsky (2006) TURBOMOLE, and MOLGW gives the exact $GW$ self-energy at all frequency points because it does not rely on any approximation and numerical parameter.

II.4 Spin Contamination

In unresticted HF and KS calculations for open-shell systems, the expectation value of the total angular momentum $\langle S^{2}\rangle$ is given by

[TABLE]

where $N_{\uparrow}$ and $N_{\downarrow}$ are the numbers of $\uparrow$- and $\downarrow$-spin electrons, respectively, and $S$ is (${N_{\uparrow}-N_{\downarrow}}$)/2 with $N_{\uparrow}>N_{\downarrow}$. The last two terms on the right side of Eq. (26) are called the spin contamination, which is non-negative. Cohen, Tozer, and Handy (2007); Menon and Radom (2008) 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 $\sim$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, that 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.

II.5 Auxiliary Basis Sets and Multi-thread Parallelization

In Gaussian-based $GW$, the 4-center integrals $(\mu\nu|\lambda\tau)$ in Eq. (14), which scale as $O(N^{4})$, are a common bottleneck in gKS and $GW$ parts in terms of compute time and memory storage. One way to reduce the bottleneck is the resolution-of-identity (RI) approximation (the density-fitting approximation), which expands the product of basis functions $\phi_{\mu}(\mathbf{r})\phi_{\nu}(\mathbf{r})$ as a linear combination of auxiliary basis functions $\phi_{\text{P}}(\mathbf{r})$. Ren et al. (2012); Hill and Platts (2008); Hill and Peterson (2012) There are two types of the 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 $(\mu\nu|\lambda\tau)$ in Eq. (14) approximate to

[TABLE]

where $P$ and $Q$ run over auxiliary basis functions, $(\mu\nu|P)$ and $(Q|\lambda\tau)$ are the 3-center integrals, and $(P|Q)$ are the 2-center integrals.

RI can be applied to both gKS [$J$ and $K^{\sigma}$ in Eqs. (13) and (16)] and $GW$ [$\mathbf{A}$, $\mathbf{B}$, $\Sigma^{\sigma}_{\text{x}}$, and $\Sigma^{\sigma}_{\text{c}}(\omega)$ in Eqs. (19), (20), (22), (23), (24), and (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 $\sim$100 times as shown in Table 1).

RI is an approximation, so it causes an error. There are mixed results for the RI error in the literature, ranging from $\sim$1 meV to $\sim$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. Körbel et al. (2014); van Setten et al. (2015)

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), (22), (23), (24), and (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 Message Passing Interface by using shared-memory threads. The performance gain by our OpenMP parallelization is shown in supplementary material. We also optimized our OpenMP implementation to reduce Non-Uniform Memory Access (NUMA) effects in modern multi-core 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.

II.6 $G_{0}W_{0}$ Quasiparticle Energy

In this work, we used three methods to obtain $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$, as it is practically impossible to obtain unique, correct, and accurate $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ 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 $\sigma$ and the real part Re are omitted for simplicity.

The first method is to linearize the non-linear QP equation in Eq. (6):

[TABLE]

where $\epsilon^{\text{G}_{0}\text{W}_{0},\text{linear}}_{m}$ is the perturbative one-shot QP energy obtained from the linearization, and $Z^{\text{linear}}_{m}$ is the QP renormalization factor for the linearization:

[TABLE]

where the derivative of the self-energy is obtained from the finite difference method using two frequency points at $\epsilon_{m}\pm\Delta\omega$ with $\Delta\omega$ being the frequency grid spacing, which is set to 0.001 Ha 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. $\epsilon_{m}\pm 0.1$ eV and $\epsilon_{m}\pm 0.5$ eV in Refs. Maggio et al., 2017; Hybertsen and Louie, 1986, respectively). In the PPA $G_{0}W_{0}$ method, different frequency points give similar results for $\epsilon^{\text{G}_{0}\text{W}_{0},\text{linear}}_{m}$, because PPA makes $\langle\varphi_{m}|\Sigma_{\text{c}}(\omega)|\varphi_{m}\rangle$ in Eq. (24) smooth around $\epsilon_{m}$ by drastically reducing the number of self-energy poles. Hybertsen and Louie (1986) However, in the full-frequency $G_{0}W_{0}$ method, different frequency points can give very different results for $\epsilon^{\text{G}_{0}\text{W}_{0},\text{linear}}_{m}$ when the finite difference method fails due to weak self-energy poles around $\epsilon_{m}$, which will be discussed later. Second, the derivative of the self-energy can be evaluated analytically, van Setten, Weigend, and Evers (2013) but it is not used in this work. Last, $Z^{\text{linear}}_{m}$ in Eq. (29) is slightly different from that in Ref. van Setten et al., 2015,

[TABLE]

The derivative of the self-energy is evalutated at $\omega=\epsilon_{m}$ and $\omega=\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ in $Z^{\text{linear}}_{m}$ and $Z_{m}$, respectively. Generally, $Z_{m}$ is smaller than $Z^{\text{linear}}_{m}$ because the self-energy has a steeper slope at $\omega=\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ than at $\omega=\epsilon_{m}$. $Z_{m}$ represents the QP weight (the pole residue of the 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 non-linear QP equation in Eq. (6) using the secant (quasi-Newton) method. In this work, we refer to $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ obtained from the graphical solution as $\epsilon^{\text{G}_{0}\text{W}_{0},\text{graph}}_{m}$. Note that the above linearization method corresponds to the first step of the secant method. Note also that when the non-linear QP equation in Eq. (6) has multiple solutions, the secant method gives only one of them, which depends strongly on the choice of $\eta$.

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

II.7 $G_{n}W_{0}$ and $G_{n}W_{n}$ Quasiparticle Energy

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

In this work, we obtained $G_{n}W_{0}$ and $G_{n}W_{n}$ QP energies ($\epsilon^{\text{G}_{n}\text{W}_{0}}_{m}$ and $\epsilon^{\text{G}_{n}\text{W}_{n}}_{m}$, respectively) by iterating the recurrence relations ($n\geq 3$):

[TABLE]

where $\epsilon^{\text{evGW},n}_{m}$ is $\epsilon^{\text{G}_{n}\text{W}_{0}}_{m}$ or $\epsilon^{\text{G}_{n}\text{W}_{n}}_{m}$, and $Z^{\text{evGW}}=1$. Summing up Eqs. (31), (32), and (33), we get

[TABLE]

which we refer to as the ev$GW$ QP equation in this work. When the ev$GW$ convergence is reached ($\epsilon^{\text{evGW},n}_{m}=\epsilon^{\text{evGW},n-1}_{m}=\epsilon^{\text{evGW},\infty}_{m}$, where $\epsilon^{\text{evGW},\infty}_{m}$ are converged ev$GW$ QP energies), the ev$GW$ QP equation in Eq. (34) becomes similar to the $G_{0}W_{0}$ QP equation in Eq. (6).

Whereas most $GW$ codes use $0<Z^{\text{evGW}}<1$, Shishkin and Kresse (2007); Blase, Attaccalite, and Olevano (2011); Véril et al. (2018) MOLGW uses $Z^{\text{evGW}}=1$. Even though we implemented ev$GW$ with $0<Z<1$ into MOLGW, we adopted ev$GW$ with $Z=1$ in this work for a few reasons. First, Eq. (33) shows that converged ev$GW$ QP energies ($\epsilon^{\text{evGW},\infty}_{m}$) is 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$ multi-solution issue from the graphical-solution and spectral-function methods and the $\sim$0.1–1 eV error from the linearization method (to be discussed in detail later). Third, ev$GW$ with $0<Z<1$ is suited for a simplified ev$GW$ variant that updates only a few states near HOMO and LUMO and rigidly shifts all the other states for efficiency, Shishkin and Kresse (2007); Blase et al. (2016) but we updated all eigenvalues in this work for accuracy. Last, ev$GW$ with $Z=1$ has no variant and does not need a QP equation solver, but ev$GW$ with $0<Z<1$ has multiple variants, depending on the choice of QP equation solvers. For example, two ev$GW$ variants with $0<Z<1$ using the linearization and graphical-solution methods in Refs. Blase, Attaccalite, and Olevano, 2011; Véril et al., 2018, respectively, may give different QP energies because the two QP equation solvers give different solutions (especially, for states far away from HOMO and LUMO).

$Z^{\text{evGW}}=1$ and the efficiency comparison of ev$GW$ and $G_{0}W_{0}$ are discusssed in supplementary material.

III Computational Details and Test Results

III.1 Computational Details

Our gKS calculations were carried out using both MOLGW and NWChem in order to cross-check the results and to ascertain the correct SCF convergence. For $GW$ calculations, we used only MOLGW. MOs were expanded using augmented Dunning correlation-consistent Gaussian basis sets, aug-cc-pV$n$Z ($n$ = D, T, Q, and 5), which are designed to smoothly converge with basis size. Augmentation using diffuse functions is essential in ground-state calculations for anions and in excited-state calculations for both neutrals and anions. Without augmentation, gKS and $GW$ eigenvalues for empty states converge very slowly with basis size. Körbel et al. (2014); Wilhelm, Del Ben, and Hutter (2016) 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. van Setten et al., 2015, and aug-cc-pVQZ in this work. Table 1 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:

[TABLE]

where $E^{\text{HF},\sigma}_{\text{x}}$, $E^{\text{PBE},\sigma}_{\text{x}}$, and $E^{\text{PBE},\sigma}_{\text{c}}$ are Fock exact exchange, PBE exchange, and PBE correlation energies, respectively. We refer to the hybrid functionals in Eq. (35) as PBE$\alpha$ functionals in this work. While we tested other functionals such as B3LYP, HSE06, BHLYP, and HF, we discuss only PBE$\alpha$($0.00\leq\alpha\leq 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 supplementary material, HSE06, PBE0, and BHLYP$\alpha$($\alpha$=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 supplementary material, PBE$\alpha$($\alpha$=1.00) and HF can make a large difference ($\sim$1 eV) in $GW$ results for some states.

III.2 gKS Test Results

III.2.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 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. Collins, Dyall, and Schaefer (1995); Peterson and Puzzarini (2005) The test results are presented in 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), Pyykkö and Desclaux (1979); Pyykkö (1988, 2012) which is consistent with Ref. Bagus, Lee, and Pitzer, 1975, and (ii) ECP $GW$ 3$d$-electron binding energies are smaller than AE ones by 0.16–0.66 eV due to indirect effects, which is consistent with Ref. Tatewaki, Yamamoto, and Hatano, 2017. 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 ($\sim$1–2 eV) errors that we encounter. However, we discuss both AE and ECP results for 3$d$-electron binding energies, where scalar relativistic effects are considerable ($\sim$10% of the experimental 3$d$-electron binding energy). Note that Ref. van Setten et al., 2015 used only AE even though the $GW$100 benchmark set contains Ag2, Cu2 and NCCu molecules.

III.2.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 supplementary material. We found that CN=5 RI ECP causes a large random error in gKS results (e.g. $\sim$0.2 and $\sim$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 non-zero 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. Wilhelm, Del Ben, and Hutter, 2016. Note also that we discussed the effect of RI on $GW$ results in supplementary material.

III.2.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 as well as wrong convergence to some local minima. This is a particularly critical issue in gKS calculations on open-shell systems involving non-zero EXX and large basis. Many minima with similar total energies and total spins due to nearly degenerate energy levels in 3$d$ transtion metals make it more difficult to obtain correct SCF convergence. Gutsev, Rao, and Jena (2000); Gonzales, King, and Schaefer (2000); Kudin, Konstantin N. and Scuseria, Gustavo E. (2007) 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 $\langle S^{2}\rangle$ 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 supplementary material.

Note that because of our heavy SCF convergence tests, total gKS calculations are more expensive than total $GW$ calculations in this work.

III.3 $GW$ Test Results

III.3.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 $\sim$0.1 eV. Ren et al. (2012) We, therefore, tested the effect of fitting function type, EXX amount, and basis size on the $GW$ CBS limit.

Two fitting functions are most widely used for the CBS limit, Hättig et al. (2012) which we refer to as standard fitting functions in this work:

[TABLE]

where $E_{m}$ are correlation or $m$th QP energies, $a$ and $b$ are fitting parameters, $N_{\text{BF}}$ is the number of basis functions (see Table 1), 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 non-standard fitting functions used in the literature. Hüser, Olsen, and Thygesen (2013); Rangel et al. (2016); Bruneval et al. (2016); Hung et al. (2017a); Shi et al. (2018)

Fig. 1 compares CBS results obtained from two standard fitting functions in Eqs. (36) and (37) (as well as one non-standard one used in Refs. Bruneval et al., 2016; Shi et al., 2018). We see that different fitting functions always give different $GW$ CBS limits, deviating from each other by up to $\sim$0.1 eV depending on molecular systems and molecular orbitals. Fig. 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 ($\sim$0.1 eV) in the $GW$ CBS limit, which is commonly observed in the literature. van Setten, Weigend, and Evers (2013); van Setten et al. (2015) 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. van Setten et al., 2015 reported that IEs obtained from Gaussian- and planewave (PW)-based $GW$ implementations with and without extrapolation, respectively, differ by $\sim$0.2 eV, but Refs. Maggio et al., 2017; Govoni and Galli, 2018 showed that the use of PW $GW$ IEs with extrapolation reduces the difference to $\sim$0.06 eV.

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.

III.3.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 supplementary material. Here, we briefly mention a couple of trends. In most cases, CN=5 has a small ($\sim$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 ($\sim$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 ($\sim$100 empty states per atom, as shown in Table 1) 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. van Setten et al., 2015, which used only CN=3,4 for extrapolation, and gives a rough estimate for two important and inter-dependent convergence parameters for $\Sigma^{\sigma}_{\text{c}}(\omega)$, the dimension of the dielectrix matrix and the number of empty states, in sum-over-states PW $GW$ implementations. van Setten et al. (2015)

The above conclusion holds only for occupied states. The effect of CN=5 on the $GW$ CBS limit for empty states is discussed in 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, Bruneval and Marques (2013) is also discussed in supplementary material.

III.3.3 $G_{0}W_{0}$ Quasiparticle Energy

A full-frequency $G_{0}W_{0}$ method used in this work produces complicated self-energy pole (and spectral-function peak) structures at non-frontier orbitals, so it is not straightforward to automatically obtain correct and accurate $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ by using a single value of $\eta$ and a single QP equation solver. Thus, we used three values of $\eta$ (0.001, 0.002, and 0.005 Ha with $\Delta\omega$ = 0.001 Ha) and three QP equation solvers mentioned in Section II.6 (linearization, graphical-solution, and spectral-function methods).

Our analysis of a total of nine solutions for $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ shows that the graphical-solution method using $\eta$ = 0.001 Ha and the linearization method randomly give incorrect solutions. Therefore, we obtained $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ from the graphical-solution or spectral-function method using $\eta$ = 0.002 or 0.005 Ha. When $\epsilon^{\text{G}_{0}\text{W}_{0},\text{graph}}_{m}$ and $\epsilon^{\text{G}_{0}\text{W}_{0},\text{spect}}_{m}$ are different, we manually selected a correct solution by analyzing $\Sigma_{\text{c}}(\omega)$ and $A(\omega)$.

A large distance between $\epsilon_{m}$ and $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ and multiple self-energy poles between them makes it difficult to obtain correct and accurate $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ for non-frontier orbitals, which is especially the case for (semi-)local xc functionals. Fig. 3 shows various examples of successes and failures of three QP equation solvers for $G_{0}W_{0}@\text{PBE}$. In the following, we analyze each example individually. Note that we used $\eta$ = 0.0001, 0.0002, and 0.0005 Ha with $\Delta\omega$ = 0.0001 Ha in a few of the following examples to demonstrate the danger of small $\eta$ values for the spectral-function method. Such small $\eta$ values are not recommended, as they significantly increase disk storage requirements while barely improving accuracy.

First, the top two left panels of Fig. 3 show a general example, in which all three methods succeed. We see a few general trends. Typically, all three methods give correct solutions at $m$ = HOMO and LUMO, which have a simple pole structure in $\Sigma_{\text{c}}(\omega)$. 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 $\omega=-0.48$ eV), whereas the linearization method gives a different (correct but inaccurate) solution (in this case, at $\omega=-0.69$ eV) due to an intrinsic error of $\sim$0.1–1 eV. A very small $\eta$ (0.0002 Ha) sharpens a weak self-energy pole at $\omega=-1.0$ eV in $\Sigma_{\text{c}}(\omega)$, but the sharpened pole does not cause an error in the graphical-solution method because it is not between $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ and $\epsilon_{m}$. The very small $\eta$ heightens the weak peak C in $A(\omega)$, 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 $\eta$ than the graphical-solution method. The very small $\eta$ has little effect on the linearization method, because (i) the linearization method in this work depends only on $\Sigma_{\text{c}}(\epsilon_{m}\pm\Delta\omega)$, and (ii) $\omega=-1.0$ eV is too distant from $\epsilon_{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 $\Sigma_{\text{c}}(\omega)$. A very small $\eta$ (0.0002 Ha) sharpens a weak pole at $\omega=-2.8$ eV in $\Sigma_{\text{c}}(\omega)$, and this sharpened pole causes a large error of 0.4 eV in the graphical-solution method because it is between $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ and $\epsilon_{m}$ and sharpened enough to make the secant method fail by causing it to find an incorrect intersection point. The very small $\eta$ heightens a weak peak B in $A(\omega)$, 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 $\eta$ has little effect on the linearization method because $\omega=-2.8$ eV is distant from $\epsilon_{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 $\eta$ (0.0002 Ha) sharpens a weak pole at $\omega\approx\epsilon_{m}$ in $\Sigma_{\text{c}}(\omega)$, but the sharpened pole does not cause an error in the graphical-solution method even though it is between $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ and $\epsilon_{m}$, because the pole is not sharpened enough ($\eta$ = 0.0001 Ha, on the other hand, causes a large error of 1.3 eV). The very small $\eta$ heightens a weak peak at $\omega\approx\epsilon_{m}$ in $A(\omega)$, but the heightened peak cannot cause an error in the spectral-function method because it is still lower than the other peak at $\omega\approx\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$. The very small $\eta$ sharpens a weak pole at $\omega\approx\epsilon_{m}$ in $\Sigma_{\text{c}}(\omega)$, and the sharpened pole causes a large error of 2.1 eV in the linearization method because it is too close to $\epsilon_{m}$, making the finite difference method fail by causing a large error in the slope of the tangent line at $\omega=\epsilon_{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 $\eta$ (0.0002 Ha) sharpens a weak pole at $\omega=-2.3$ eV in $\Sigma_{\text{c}}(\omega)$, but the sharpened pole does not cause an error in the graphical-solution method because it is not between $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ and $\epsilon_{m}$. Two peaks A and B (at $\omega=-3.9$ and $-2.1$ eV, respectively) in $A(\omega)$ 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: the peak B (A) is higher than the peak A (B) for CN=2,3,4 (CN=5). We chose peak B as the QP peak, because (i) it is consistent with the solution from the graphical-solution method, and (ii) it is consistent with a trend over EXX amount [$G_{0}W_{0}@\text{PBE}\alpha(0.00\leq\alpha\leq 1.00)$ QP HOMO-2 energies of CuO*-* using the solution from the peak B vary smoothly with $\alpha$, as shown in Fig. 6]. The choice of $\eta$ and CN has little effect on the linearization method, but $\epsilon^{\text{G}_{0}\text{W}_{0},\text{linear}}_{m}=-2.9$ eV causes a large overestimation error of 0.8 eV in $G_{0}W_{0}@\text{PBE}$ binding energy.

There are several 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 supplementary material, (ii) not only a very small $\eta$, but also a very large $\eta$ (e.g. $\sim$0.05 Ha in Ref. Hung et al., 2017b) can cause a large error, (iii) deep states (e.g. HOMO-$n$, where $n$ = 5, 6, …) have much more complicated pole [peak] structures in $\Sigma_{\text{c}}(\omega)$ [$A(\omega)$] than those in Fig. 3, so it is very difficult to choose correct and accurate $\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ for deep states not only automatically, but also manually.

We conclude this section by summarizing several guidelines to obtain a reliable and reproducible $G_{0}W_{0}@\text{PBE}$ QP spectrum. First, one should try multiple $\eta$ (and $\Delta\omega$) values. There is no single general $\eta$ value that works well for all QP equation solvers, molecular systems, and molecular orbitals. In other words, while $\eta$ is typically viewed as a convergence parameter (the smaller $\eta$, the more accurate $GW$ QP energy), it is practically an adjustable parameter, which should be not too small or too large (e.g. 0 and $\sim$0.05 Ha in Refs. Véril et al., 2018; Hung et al., 2017b, respectively). The optimal value of $\eta$ depends on $|\epsilon^{\text{G}_{0}\text{W}_{0}}_{m}$ - $\epsilon_{m}|$, which generally decreases with the amount of EXX and increases with the depth of the $m$th state. For example, when calculating $G_{0}W_{0}@\text{PBE}$ HOMO and LOMO (lowest occupied molecular orbital) energies, one may try $\sim$0.1 and $\sim$1 eV, respectively, for $\eta$.

Second, we recommend using multiple QP equation solvers. As shown in Fig. 3, the $G_{0}W_{0}@\text{PBE}$ QP spectrum automatically obtained from a single QP equation solver can contain a large ($\sim$1 eV) error at random states.

Third, we recommend using multiple basis sizes. As shown in the bottom right of Fig. 3, different basis sets with different sizes can give very different $G_{0}W_{0}$ QP energies (by $\sim$1 eV) at random states. Using multiple basis sizes allows for not only accurate $GW$ results without small ($\sim$0.1 eV) systematic errors from the basis set incompleteness, but also correct gKS and $GW$ results without large ($\sim$1 eV) random errors from SCF convergence and $GW$ multi-solution 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. Ref. Maggio et al., 2017 suggests the linearization method as a preferable method for a fair comparison of $G_{0}W_{0}@\text{PBE}$ IE (and EA) from different $GW$ implementations, because it gives a unique solution and thus is free of the $GW$ multi-solution 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 $\sim$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 $\sim$0.5 eV underestimation error by $G_{0}W_{0}@\text{PBE}$ with respect to experiment. van Setten et al. (2015); Maggio et al. (2017); Govoni and Galli (2018) However, for deep states, it randomly succeeds or fails, as shown in the bottom left of Fig. 3, and randomly overestimates or underestimates $G_{0}W_{0}@\text{PBE}$ binding energies by $\sim$1 eV compared to accurate ones, as shown in the bottom right of Fig. 3 and supplementary material, respectively. This large and unpredictable (with respect to state, magnitude, and direction) error makes the linearization method inadequate for the $G_{0}W_{0}@\text{PBE}$ QP spectrum.

Last, one should be aware that different ways to handle the $GW$ multi-solution issue are found in the literature. For example, Ref. van Setten et al., 2015, which suggests that all solutions are physically relevant, manually searched for an actually relevant one by varying $\eta$, whereas Ref. Véril et al., 2018 avoided the issue by automatically selecting the solution with the largest $Z_{m}$ in Eq. (30) and using only $\eta=0$. In this work, we adopted a combined approach. In other words, we automatically chose the solution with the highest spectral weight, which is identical to the solution with the largest $Z_{m}$, as explained in Section II.6, 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 $G_{0}W_{0}$ starting point and ev$GW$ self-consistency level without causing unphysical kinks when multiple solutions are equally relevant [e.g. two solutions at peaks A and B in the bottom two right panels of Fig. 3 give similar $Z_{m}$ ($\sim$0.2) for $m=\text{HOMO-2}$ of CuO*-* due to similar slopes of the self-energy (approximately, $-4.0$)]. However, when it comes to $\eta$, we adopted the approach of Ref. van Setten et al., 2015 instead of that of Ref. Véril et al., 2018 because $\eta=0$ frequently causes the secant method in the graphical-solution method to find an incorrect intersection point and makes $\langle\varphi^{\sigma}_{m}|\Sigma^{\sigma}_{\text{c}}(\omega)|\varphi^{\sigma}_{m}\rangle$ in Eq. (24) diverge at $\omega=\epsilon^{\sigma}_{i}-\Omega_{s}$ and $\omega=\epsilon^{\sigma}_{a}+\Omega_{s}$. The cumulant expansion, Hedin (1999); Zhou et al. (2015) which describes plasmon satellites better than the $GW$ approximation, may allow us to address the $GW$ multi-solution issue when the QP picture breaks down, but it is beyond the scope of this work.

III.3.4 $G_{n}W_{0}$ and $G_{n}W_{n}$ Quasiparticle Energy

In this work, we used only $\eta=0.001$ Ha for ev$GW$ because it is small enough to obtain the convergence of ev$GW$ QP energies with respect to $\eta$ within $\sim$0.01 eV. The convergence test results for ev$GW$ QP energies with respect to the iteration number are shown in supplementary material. QS$GW$ and our ev$GW$ are quasiparticle-only $GW$ schemes with no spectral weight transfer ($Z=1$), and $G_{n}W_{n}$ is a diagonal approximation to QS$GW$. Therefore, we compared the convergence behaviors of QS$GW$ and our ev$GW$ and found a couple of similarities and differences between them.

First, the ev$GW$ convergence is reached after only a few iterations, which is consistent with the literature. Shishkin and Kresse (2007); Rangel et al. (2016); Blase et al. (2016) Due to the fast and stable convergence, a mixing scheme is not used in our ev$GW$. Unlike ev$GW$, QS$GW$ generally needs $\sim$10–20 (up to 60) iterations and a mixing scheme. Bruneval (2012); Koval, Foerster, and Sánchez-Portal (2014); Kaplan et al. (2016)

Second, the orbital character affects the starting-point dependency of ev$GW$. For example, we observed that ev$GW$ QP energies for HOMO of CuO*-* depend more strongly on the EXX amount in the ev$GW$ starting point than those for HOMO of ScO*-. We attribute this to different amounts of 3$d$ character in HOMOs of ScO-* and CuO*-* (6% and 23%, respectively, as shown in Table 2). In other words, as the 3$d$ character in MO increases, the starting-point dependency of ev$GW$ increases. We also observed that $G_{n}W_{n}$ has a weaker (stronger) starting-point dependency for HOMO of ScO*-* (CuO*-) than $G_{n}W_{0}$. Our observations for ScO-* are consistent with Ref. Blase et al., 2016, which studied the ev$GW$ starting-point dependency using small water clusters and concluded that as the ev$GW$ self-consistency level increases from $G_{n}W_{0}$ and $G_{n}W_{n}$, the ev$GW$ starting-point dependency decreases. However, our observations for CuO*-* are not consistent with this conclusion. This is likely because CuO*-* has strong 3$d$ character in HOMO, whereas ScO*-* and small water clusters do not. This orbital-character-dependent starting-point dependency of ev$GW$ may be related to conflicting results for QS$GW$ in the literature: Ref. Koval, Foerster, and Sánchez-Portal, 2014 showed the starting-point independency of QS$GW$ using a small $sp-$bonded molecule (CH4), while Ref. Liao and Carter, 2011 showed the strong starting-point dependency of QS$GW$ using a $d$ solid ($\alpha$-Fe2O3).

IV Results and Discussion

In this section, we compare our $GW$ calculations to anion PES experiments, Wu and Wang (1998, 1997); Wu, Desai, and Wang (1997); Moravec et al. (2001) focusing especially on the first IE, the lowest 3$d$-electron binding energy, and the orbital order. We present our results from two approaches seperately: First, we discuss non-self-consistent $GW$ with different starting-points (namely, $G_{0}W_{0}@\text{PBE}\alpha$ calculations as $\alpha$ is varied in steps of 0.25 from 0 to 1), and then, we discuss eigenvalue self-consistent $GW$ ($G_{n}W_{0}$ and $G_{n}W_{n}$) with PBE starting point. We only briefly discuss our $GW$ results for the starting-point–self-consistency hybrid approach, because (i) fundamentally, we found that the hybrid approach does not give any better results than the two separate approaches, and (ii) practically, the hybrid approach inherits disadvantages from both approaches.

ScO*-, TiO-, CuO-, and ZnO-* are similar but different systems in several aspects. First, ScO*-* and CuO*-* are closed-shell systems, whereas TiO*-* and ZnO*-* are open-shell systems. Second, ScO*-* and TiO*-* have partially filled 3$d$ shells, while CuO*-* and ZnO*-* have completely filled 3$d$ shells. Third, TiO*-* has a shallow 3$d$ state, but CuO*-* and ZnO*-* have deep 3$d$ states. Fourth, 3$d$-electron photodetachment transitions are observed in TiO*-* and CuO*-, but not in ScO-* and ZnO*-. Fifth, CuO-* has strong 3$d$ character in HOMO, but ScO*-, TiO-, and ZnO-* have weak 3$d$ character in HOMO. Last, two-electron transitions are observed in ScO*-, but not in TiO-, CuO-, and ZnO-*. Due to these similarities and differences, TMO anions are an ideal set of systems for assessment of the performance of $GW$ schemes.

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

Throughout this work, we use only gKS-PBE TM 3$d$ character except when we discuss the subtle competition between direct and indirect relativistic effects, because the EXX amount has a small effect on TM 3$d$ character. Also, throughout this work, we use only the gKS-PBE orbital order to avoid confusion. The orbital order depends strongly on the amount of EXX (e.g. PBE vs HF) and the level of theory (e.g. DFT vs $GW$). van Setten, Weigend, and Evers (2013); Govoni and Galli (2018) For example, gKS-PBE HOMO of CuO*-* corresponds to gKS-PBE$\alpha$($\alpha$=1.00) HOMO-1 and $G_{0}W_{0}@\text{PBE}\alpha$($\alpha$=1.00) HOMO of CuO*-*, as shown in supplementary material and Fig. 6. In this work, gKS-PBE and $G_{0}W_{0}@\text{PBE}$ were found to have the same orbital order.

IV.1 $G_{0}W_{0}$ Starting Points

Figs. 5 and 6 show PES and $G_{0}W_{0}@\text{PBE}\alpha(0.00\leq\alpha\leq 1.00)$ QP spectra of ScO*-, TiO-, CuO-, and ZnO-*. In PES spectra, vertial dashed and solid lines represent experimental $sp$- and $d$-electron binding energies, respectively. In $GW$ spectra, oblique dashed and solid lines track calculated $sp$- and $d$-electron binding energies, respectively. In Figs. 5 and 6, we find a few general trends common in all TMO anions considered in this work: (i) no $G_{0}W_{0}@\text{PBE}\alpha(0.00\leq\alpha\leq 1.00)$ results are in perfect agreement with experiment, (ii) $G_{0}W_{0}@\text{PBE}$ underestimates the IE of TMO anions by $\sim$1 eV, which is larger than the typical underestimation for $sp$ molecules ($\sim$0.5 eV), van Setten et al. (2015); Maggio et al. (2017); Govoni and Galli (2018) and (iii) $G_{0}W_{0}@\text{PBE}\alpha(0.25\leq\alpha\leq 0.50)$ reduces it to $\sim$0.1 eV. In the following, we analyze each TMO anion individually.

IV.1.1 ScO*-*

Scandium is the first transition metal and has only one 3$d$ electron. DFT and CCSD(T) calculations in Refs. Gonzales, King, and Schaefer, 2000; Gutsev, Rao, and Jena, 2000 confirmed the ground state of ScO*-* as ${}^{1}\Sigma^{+}$ ($8\sigma^{2}3\pi^{4}9\sigma^{2}$), correcting the wrongly assumed state ${}^{3}\Delta^{-}$ ($8\sigma^{2}3\pi^{4}9\sigma^{1}1\delta$) in Ref. Wu and Wang, 1998. There is no 3$d$ peak or band in the PES spectrum of ScO*-*, and the top three valence molecular orbitals have weak Sc 3$d$ character (6%, 20%, and 18%, respectively), as shown in Table 2.

In the left panel of Fig. 5, we see that for HOMO-1 and HOMO-2 of ScO*-, $G_{0}W_{0}@\text{PBE}$ binding energies slightly overestimate PES ones (by 0.20 and 0.02 eV for HOMO-1 and HOMO-2, respectively), whereas $G_{0}W_{0}@\text{PBE}\alpha(0.25\leq\alpha\leq 1.00)$ binding energies significantly overestimate PES ones by $\sim$2 eV [e.g. for HOMO-1, $G_{0}W_{0}@\text{PBE}\alpha$($\alpha$=0.25) and $G_{0}W_{0}@\text{PBE}\alpha$($\alpha$=1.00) binding energies are greater than PES ones by 1.67 and 2.72 eV, respectively]. This seems to suggest that for HOMO-1 and HOMO-2 of ScO-, $G_{0}W_{0}@\text{PBE}$ performs better than $G_{0}W_{0}@\text{PBE}\alpha(0.25\leq\alpha\leq 1.00)$, but this is not the case due to the nature of the corresponding peaks in the PES experiment. Ref. Gonzales, King, and Schaefer, 2000 suggests that the second and third peaks in the PES epectrum of ScO-* are likely due to two-electron transitions from $8\sigma^{2}3\pi^{4}9\sigma^{2}$ (${}^{1}\Sigma^{+}$ ScO*-) to $8\sigma^{2}3\pi^{4}10\sigma$ ($B^{2}\Sigma^{+}$ ScO) and to $8\sigma^{2}3\pi^{4}1\delta$ ($A^{\prime 2}\Delta$ ScO) states, respectively, which $GW$ calculations for quasiparticle excitations cannot account for. In other words, the seemingly excellent agreement between $G_{0}W_{0}@\text{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 $\alpha$ increases, $G_{0}W_{0}@\text{PBE}\alpha$ IE always increases, but this happens at different rates: As $\alpha$ increases from 0.00 to 0.25, it increases rapidly, whereas as $\alpha$ increases from 0.25 to 1.00, it increases slowly. The weak sensitivity of $G_{0}W_{0}@\text{PBE}\alpha(0.25\leq\alpha\leq 1.00)$ IE to a change in $\alpha$ gives a large margin for an optimal amount of EXX: 25%–100%.

IV.1.2 TiO*-*

Titanium is the second transition metal and has two 3$d$ electrons. Several theoretical studies in Refs. Walsh, King, and Schaefer, 1999; Gutsev, Rao, and Jena, 2000; Gutsev, Andrews, and Bauschlicher, 2003; Dai et al., 2003; Miliordos and Mavridis, 2010 confirmed $9\sigma^{2}\delta^{1}$ (${}^{2}\Delta$) as the ground-state electron configuration of TiO*-, correcting the wrongly assigned configuration $9\sigma^{1}\delta^{2}$ (${}^{1}\Sigma^{+}$) in Ref. Wu and Wang, 1997. Unlike ScO-, which has an empty $\delta$ shell, TiO-* has one 3$d$ electron in the $\delta$ shell. The transition of the 3$d$ electron from $9\sigma^{2}\delta^{1}$ (${}^{2}\Delta$ TiO*-) to $9\sigma^{2}$ (${}^{1}\Sigma^{+}$ TiO) states produces the third peak in the PES spectrum of TiO-* at 2.0 eV. In the $G_{0}W_{0}@\text{PBE}\alpha$ QP spectrum of TiO*-*, $\uparrow$-HOMO is of entirely Ti 3$d$ character, as shown in Table 2.

The right panel of Fig. 5 clearly shows that $G_{0}W_{0}@\text{PBE}\alpha(0.00\leq\alpha\leq 1.00)$ binding energy for $\uparrow$-HOMO of TiO*-* is much more sensitive to a change in $\alpha$ than those for other occupied molecular orbitals with mainly $sp$ character, as shown in Table 2. The orbital-character-dependent sensitivity of $G_{0}W_{0}@\text{PBE}\alpha$ binding energy to a change in $\alpha$ causes a couple of problems. First, $G_{0}W_{0}@\text{PBE}$ underestimates the IE and the 3$d$-electron binding energy of TiO*-* non-uniformly (by 0.99 and 1.74 eV, respectively), leading to the wrong orbital order. In other words, $G_{0}W_{0}$ does not correct the wrong orbital order produced by PBE. Second, the $G_{0}W_{0}$ starting-point approach does not give accurate results for both the IE and the 3$d$-electron binding energy of TiO*-* at the same time. For example, $G_{0}W_{0}@\text{PBE}\alpha$($\alpha$=0.50) gives a better result for the IE of TiO*-* by 0.29 eV, but a worse result for the 3$d$-electron binding energy of TiO*-* by 1.46 eV, than $G_{0}W_{0}@\text{PBE}\alpha$($\alpha$=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. Shishkin and Kresse (2007); Fuchs et al. (2007); Klimeš, Kaltak, and Kresse (2014)

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

IV.1.3 CuO*-*

Copper is the 11th transition metal and has ten 3$d$ electrons. DFT calculations in Ref. Gutsev, Rao, and Jena, 2000 confirmed the ground state of CuO*-* as ${}^{1}\Sigma^{+}$ with the electron configuration of $3d^{10}2p\sigma^{2}2p\pi^{4}$. There are three bands (named as X, Y, and Z in Ref. Wu, Desai, and Wang, 1997) in the PES spectrum of CuO*-, as shown in the top of the left panel of Fig. 6. Ref. Wu, Desai, and Wang, 1997 suggested that the photodetachment transition of 3$d$ electrons ($3d\delta^{4}3d\pi^{4}3d\sigma^{2}$) from ${}^{1}\Sigma^{+}$ CuO-* $3d^{10}2p\sigma^{2}2p\pi^{4}$ to Z CuO $3d^{9}2p\sigma^{2}2p\pi^{4}$ states produces the broad Z band in the PES spectrum of CuO*-* at $\sim$4.5 eV (which we selected from the position of the highest peak in the Z band) and assumed that the Z band is unusually broad likely due to a large geometry change from the anion to the neutral. In the $G_{0}W_{0}@\text{PBE}\alpha$ QP spectrum of CuO*-*, HOMO-2 is of entirely Cu 3$d$ character, as shown in Table 2.

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

IV.1.4 ZnO*-*

Zinc is the 12th transition metal and has ten 3$d$ electrons. Zinc is rather distinct from other first row transition metals due to its closed-shell electron configuration. In other words, zinc is more similar to alkaline earth metals than other transition metals because Zn 3$d$ electrons generally do not participate in bonding. Fancher et al. (1998) DFT calculations in Ref. Gutsev, Rao, and Jena, 2000 confirmed the ground-state electron configuration of ZnO*-* as ${}^{2}\Sigma^{+}$ $10\sigma^{1}9\sigma^{2}4\pi^{4}\delta^{4}$. There are four bands in the PES spectrum of ZnO*-, as shown in the top of the right panel of Fig. 6. The photodetachment of 3$d$ electrons is not measured in the PES experiment due to insufficient photon energy of 4.66 eV. Moravec et al. (2001) Unlike CuO-, all valence molecular orbitals in the $G_{0}W_{0}@\text{PBE}\alpha$ QP spectrum of ZnO-* have weak Zn 3$d$ character, as shown in Table 2.

The right panel of Fig. 6 shows that unlike CuO*-, $G_{0}W_{0}@\text{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 $\uparrow$-HOMO, $\uparrow$-HOMO-1, $\downarrow$-HOMO, and $\downarrow$-HOMO-1, respectively) possibly because all valence molecular orbitals have similar Zn 3$d$ character and thus their $G_{0}W_{0}@\text{PBE}\alpha$ binding energies have similar sensitivity to a change in $\alpha$. $G_{0}W_{0}@\text{PBE}\alpha$($\alpha$=0.50) gives good results for the IE and the orbital order of ZnO*-* at the same time. Unlike CuO*-*, $G_{0}W_{0}@\text{PBE}\alpha(0.50\leq\alpha\leq 1.00)$ binding energies for all valence molecular orbitals, except for $\uparrow$-HOMO and $\downarrow$-HOMO-1, are strongly sensitive to a change in $\alpha$. This trend can be explained in terms of the importance of spin-splitting in open-shell molecules (as has been discussed in Ref. Shi et al., 2018 for CuO${}_{2}^{-}$ molecule).

IV.2 ev$GW$ Self-Consistency Levels

Fig. 7 shows PES and ev$GW$ QP spectra of ScO*-, TiO-, CuO-, and ZnO-*. In the following, we analyze $G_{n}W_{0}@\text{PBE}$ and $G_{n}W_{n}@\text{PBE}$ results individually.

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

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

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

IV.3 Comparison of $G_{0}W_{0}$ starting-point and ev$GW$ self-consistency approaches

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

First, $G_{n}W_{0}@\text{PBE}$ does not contain system-dependent adjustable parameters. Unlike extended systems, there is no unique amount of EXX for the $G_{0}W_{0}$ starting point, which works well for all finite systems. For example, we showed in Section IV.1 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. Blase et al. (2016) Second, $G_{n}W_{0}@\text{PBE}$ is transferable between finite and extended systems. $G_{n}W_{0}@\text{PBE}$ works well for both molecules and solids (e.g. ZnO-* anion and bulk ZnO, respectively). Shishkin and Kresse (2007) This greatly extends the range of applicability of the $GW$ method. For example, $G_{n}W_{0}@\text{PBE}$ may be applicable to solid-molecule hybrid systems such as molecular junctions and molecules adsorbed on solid surfaces. Strange et al. (2011) Also, $G_{n}W_{0}@\text{PBE}$ may be used for the study of quantum size effects in clusters because it is independent of the cluster size. Third, $G_{n}W_{0}@\text{PBE}$ is easy to use and reliable. Unlike PBE$\alpha(0.00<\alpha\leq 1.00)$, PBE is safe from the SCF convergence issue, and unlike $G_{0}W_{0}$, ev$GW$ with $Z=1$ is immune to the $GW$ multi-solution issue. Therefore, $G_{n}W_{0}@\text{PBE}$ does not need manual, time-consuming, and error-prone tests to address the two issues, which are explained in detail in Section III.

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

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

IV.4 Comparison with results in the literature

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

We begin with the $G_{0}W_{0}$ starting-point approach. Table 4 summarizes a few selected results for the optimal amount of EXX in the $G_{0}W_{0}$ starting point out of numerous results, such as Refs. Marom et al., 2012; Caruso et al., 2016, in the literature. Interestingly, we see that there is a wide range of EXX amounts from 25% to 100%, and Refs. Bruneval and Marques, 2013; Rostgaard, Jacobsen, and Thygesen, 2010 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. Maggio et al., 2017; Govoni and Galli, 2018 showed that such implementation differences have little effect on $G_{0}W_{0}$ IE ($\sim$0.06 eV). There are a couple of other factors that have a stronger effect on $G_{0}W_{0}$ results than implementation differences. One factor is the choice of system and property. As shown in Section IV.1, $G_{0}W_{0}@\text{PBE}\alpha(0.25\leq\alpha\leq 1.00)$ IEs of $sp$ systems are slightly different (by $\sim$0.1 eV). Most existing $G_{0}W_{0}$ studies used the IE of $sp$-bonded systems to determine the optimal amount of EXX in the $G_{0}W_{0}$ starting point. The other factor is that the choice of QP equation solver and CBS extrapolation method. As shown in Section III.3.1 and Section III.3.3, the linearization method and the CBS extrapolation method (e.g. whether to extrapolate or not and which fitting function and basis set to use for extrapolation) can cause a difference in $G_{0}W_{0}$ IE on the order of $\sim$0.1 eV. Overall, the combination of the two factors gives a large margin for the optimal amount of EXX in the $G_{0}W_{0}$ starting point, and thus is likely to produce the wide range of amounts that exist in the literature.

Next, we move on to the ev$GW$ self-consistency approach, and discuss the origin of apparently conflicting ev$GW$ results for IE and starting-point dependency. First, Ref. Blase, Attaccalite, and Olevano, 2011 reported that the $G_{n}W_{n}$ approach with a local-density approximation starting point gives good results for the IE of large $sp$ molecules, whereas we found in Section IV.2 that $G_{n}W_{0}@\text{PBE}$ gives satisfactory results for the electronic structure (including the IE) of small 3$d$ molecules. A comparison of ev$GW$ implementations in Ref. Blase, Attaccalite, and Olevano, 2011 and this work is provided in supplementary material. We believe that the main origin of the different results is the orbital-character-dependent sensitivity of ev$GW$ binding energy to a change in ev$GW$ self-consistency level. As shown in Section IV.2, $G_{n}W_{0}@\text{PBE}$ and $G_{n}W_{n}@\text{PBE}$ binding energies are slightly different for delocalized HOMO with weak 3$d$ character by $\sim$0.1 eV, but significantly different for localized HOMO with strong 3$d$ character by $\sim$1 eV. Unlike this work, Ref. Blase, Attaccalite, and Olevano, 2011 used the linearization method, employed pseudopotentials and RI, and did not use the CBS limit, but these cause small ($\sim$0.1 eV) differences in ev$GW$ IE, as shown in Section III.3. Accordingly, they are most likely not the reason for the large (0.98 eV) difference between $G_{n}W_{0}@\text{PBE}$ and $G_{n}W_{n}@\text{PBE}$ IEs of CuO*-. Second, Ref. Blase et al., 2016 reported that in small water clusters, as the ev$GW$ self-consistency level increases, the ev$GW$ starting-point dependency decreases, whereas we found in Section III.3.4 that in TMO anions, $G_{n}W_{n}$ sometimes depends more strongly on the starting point than $G_{n}W_{0}$. As mentioned in Section III.3.4, we believe that the orbital character influences the ev$GW$ starting-point dependency: for molecular orbitals with strong (weak) 3$d$ character, $G_{n}W_{n}$ depends more strongly (weakly) on the starting point than $G_{n}W_{0}$. Overall, without molecular orbitals with strong 3$d$ character (e.g. HOMO of CuO-*), our ev$GW$ results for IE and starting-point dependency in this work are consistent with those in Refs. Blase, Attaccalite, and Olevano, 2011; Blase et al., 2016.

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 $\uparrow$-HOMO of TiO*-* with entirely Ti 3$d$ character, (ii) we applied 15 different starting-point–self-consistency hybrid $GW$ schemes ($G_{0}W_{0}$, $G_{n}W_{0}$, and $G_{n}W_{n}$; 0%, 25%, 50%, 75%, and 100% EXX) to them, and (iii) we searched for $GW$ schemes that give a reasonably small error of less than 0.5 eV in the IE and the 3$d$-electron binding energy with respect to experiment. Fig. 8 shows the results of the test. We see that $GW$ IEs of ScO*-* and TiO*-* (red dashed lines) depend weakly on the starting point and the self-consistency level, giving a large margin for the choice of $GW$ schemes. 14 $GW$ schemes out of 15 ($G_{0}W_{0}@\text{PBE}$ is an exception as expected) give a small error (less than 0.5 eV), which explains why there are a large number of different good $GW$ schemes for the IE of $sp$ molecules in the literature. We also see that the $GW$ 3$d$-electron binding energy of TiO*-* (green solid lines) depends strongly on the starting point and the self-consistency level, yielding a small margin for the choice of $GW$ schemes. Only two $GW$ schemes ($G_{0}W_{0}@\text{PBE}$0 and $G_{n}W_{0}@\text{PBE}$) out of 15 give a small error (less than 0.5 eV), which is why we obtained a small number of good $GW$ schemes for the electronic structure of $d$ molecules in this work. Overall, we confirm that evaluation results for the performance of $GW$ schemes depend strongly on the choice of system and property (e.g. the IE with mainly $sp$ character vs the electronic structure containing $d$ states).

V Summary and Conclusions

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

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

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

We found that the eigenvalue self-consistency approaches, $G_{n}W_{0}@\text{PBE}$ and $G_{n}W_{n}@\text{PBE}$, can improve $G_{0}W_{0}@\text{PBE}$, too. Especially, $G_{n}W_{0}@\text{PBE}$ gives as good results for the IE and the $3d$-electron binding energy as $G_{0}W_{0}@\text{PBE}\alpha$ without suffering from $GW$ multi-solution and SCF convergence issues.

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

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

Our results in this work – (i) $G_{0}W_{0}@\text{PBE}\alpha$($0.25\leq\alpha\leq 0.50$) and $G_{n}W_{0}@\text{PBE}$ give good QP energies for molecular orbitals with both weak and strong 3$d$ character, and (ii) the ev$GW$ starting-point dependency is more related to the orbital character than the self-consistency level – may seem to disagree with some results in the literature, but this is not the case. The origin of the seeming disagreement is that except for $G_{0}W_{0}@\text{PBE}\alpha$($0.00\leq\alpha\leq 0.25$), varying the self-consistency level and the starting point generally makes a small ($\sim$0.1 eV) change in QP energy for HOMO with mainly $sp$ character, which is accidentally comparable to individual or combined errors from multiple sources, such as the incomplete basis set, the linearization method in $G_{0}W_{0}$, and the insufficient number of eigenvalues to update in ev$GW$.

$G_{n}W_{0}@\text{PBE}$ is not a conserving and starting-point-independent $GW$ scheme. It is not the most accurate or efficient $GW$ scheme, either. However, $G_{n}W_{0}@\text{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.

VI Supplementary Material

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

Acknowledgements.

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 also would like to thank Bin Shi and Meisam Rezaei for useful discussions in the earlier stage of this work.