Photoelectron spectra of early 3d-transition metal dioxide anions, ScO2, TiO2, VO2, CrO2, and MnO2, are calculated using semilocal and hybrid density functional theory (DFT) and many-body perturbation theory within the GW approximation using one-shot perturbative and eigenvalue self-consistent formalisms. Different levels of theory are compared with each other and with available photoelectron spectra. We show that one-shot GW with a PBE0 starting point (G0W0@PBE0) consistently provides very good agreement for all experimentally measured binding energies (within 0.1 eV–0.2 eV or less). We attribute this to the success of PBE0 in mitigating self-interaction error and providing good quasiparticle wave functions, which renders a first-order perturbative GW correction effective. One-shot GW calculations with a Perdew–Burke–Ernzerhof (PBE) starting point do poorly in predicting electron removal energies by underbinding orbitals with typical errors near 1.5 eV. A higher exact exchange amount of 50% in the DFT starting point of one-shot GW does not provide very good agreement with experiment by overbinding orbitals with typical errors near 0.5 eV. While not as accurate as G0W0@PBE0, the G-only eigenvalue self-consistent GW scheme with W fixed to the PBE level provides a reasonably predictive level of theory (typical errors near 0.3 eV) to describe photoelectron spectra of these 3d-transition metal dioxide anions. Adding eigenvalue self-consistency also in W, on the other hand, worsens the agreement with experiment overall. Our findings on the performance of various GW methods are discussed in the context of our previous studies on other transition metal oxide molecular systems.

The last decade has witnessed a growing number of computational studies that have benchmarked Green’s function methods, such as the GW approximation1 and the Bethe–Salpeter equation, for excited state properties of bulk and molecular systems.2–44 A large majority of these studies have focused on the performance of various flavors of the GW approximation in predicting the electron removal energies in sp-bonded molecules and clusters by comparing their predictions with accurate quantum chemistry calculations and experimental photoelectron spectroscopy data. Modeling excited states of systems containing transition metal elements in general, and transition metal oxides in particular, within the GW theory has faced additional theoretical and computational challenges.45 Enhanced electron correlations inherent in these systems and their propensity to having open-shell electronic configurations necessitate the use of more sophisticated approaches beyond simple perturbative implementations on top of density functional theory (DFT) with semi-local exchange-correlation functionals. Furthermore, convergence issues with respect to basis set size and other implementation parameters present significant computational bottlenecks for maintaining a desirable level of accuracy comparable to what can be achieved for sp-bonded systems that is typically ∼0.1 eV. Accordingly, there have been much fewer studies on the performance of the GW approximation for systems that contain transition metal elements. Motivated by this observation, here we continue with our recent benchmark studies37,38,42,43 by focusing on the electronic structure of negatively charged 3d-transition metal dioxide molecules TMO2 for TM = Sc, Ti, V, Cr, and Mn.

Photoelectron spectra (PES) of early 3dTMO2 molecules have been available up to photon energies of 5 eV–6.5 eV since the pioneering studies of Wang and collaborators from more than two decades ago.46–52 Their structural and electronic properties have been investigated53–72 in various computational studies using methods based on DFT and quantum chemistry. While these TMO2 molecules are isostructural with few changes in the bond angles and lengths upon changing the TM element, their frontier molecular orbitals display a wide range of spatial localization properties and have significantly varying amounts of TM 3d and O 2p contents. Therefore, their electronic structures provide a diverse set of challenges to DFT within the generalized Kohn–Sham scheme and to many-body perturbation theory within the GW approximation due to the need to mitigate self-interaction error (SIE) and the importance of a suitable DFT starting point or a form of self-consistency in their GW treatment. In this work, we model the photoelectron spectra (PES) of early TMO2 molecules using DFT with semi-local and hybrid exchange–correlation functionals and within the GW approximation using one-shot perturbative schemes with different DFT starting points as well as eigenvalue self-consistent formalisms.

This article is organized as follows: We begin with an overview of the computational methods and parameters used in our calculations in Sec. II. This is followed in Sec. III with a brief discussion of the structural and vibrational properties of the five TMO2 molecules considered in this study and detailed analyses of their PES computed with DFT and GW methods in comparison to each other and to experimental data. The trends in the performance of DFT and variants of the GW approximation are discussed in Sec. IV. Finally, we summarize our findings and analyses in Sec. V.

DFT computations to find the optimized structures were carried out with the NWCHEM code,73 version 6.6, using the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional74 and aug-cc-pVTZ basis sets. The photoelectron spectra were simulated by broadening the eigenvalue spectra with a Gaussian distribution function of 0.1 eV smearing width, without taking into account the photoionization cross sections. For DFT spectra obtained with PBE and PBE075 functionals, the Kohn–Sham eigenvalues were shifted to align the first peak with the vertical ionization potential (IP) of the anion, which is computed as the total energy difference between the anionic and neutral molecules at the fixed anionic geometry.

GW calculations were carried out using the MOLGW76 software within both the one-shot (perturbative) G0W0 scheme and the eigenvalue self-consistent GW (evGW) scheme with two types of self-consistency, GnW0 and GnWn, which update the eigenvalues only in G and in both G and W, respectively. In order to investigate the starting point dependency of the GW approximation, we used global hybrid functionals

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

where ExHF, ExPBE, EcPBE, and α are the Fock exact exchange, PBE exchange, PBE correlation energies, and amount of exact exchange, respectively. For the starting point dependency of the G0W0 scheme, we tested three different values for α: 0 (PBE starting point), 0.25 (PBE0 starting point), and 0.50. In these full-frequency G0W0 calculations, we solved the nonlinear quasiparticle equation for each state graphically using the secant (quasi-Newton) method. The numerical parameter η used to broaden the self-energy poles was chosen as η = 0.001 hartree, and the self-energy was evaluated on a frequency grid of spacing Δω = 0.001 hartree. As discussed in detail in Ref. 43, the full-frequency G0W0 method leads to complicated self-energy pole structures, typically at nonfrontier orbitals, which results in multiple solutions of the quasiparticle equation and makes determining the correct and accurate quasiparticle energies difficult. This multisolution issue of G0W0 is particularly prevalent in transition metal oxide systems with a PBE starting point. As such, we checked the accuracy of our predictions for the quasiparticle energies by (i) also using the spectral-function method, which locates the quasiparticle peak with the highest weight in the spectral function, (ii) by looking at the trends as a function of the basis set size and the amount of exact exchange, and (iii) by varying the parameters η and Δω, as recommended in Ref. 43.

For calculations performed with the GnW0 and GnWn self-consistency methods, we used the iterative scheme employed in MOLGW, described in Ref. 43, where the quasiparticle renormalization factor Z is set to Z = 1. In these calculations, we used a larger broadening parameter η = 0.01 hartree for the self-energy poles (and accordingly, a larger frequency grid spacing of Δω = 0.01 hartree) in order to avoid the oscillations typically observed in the quasiparticle energies of nonfrontier orbitals as a function of the iteration index. With these parameters, self-consistency to 0.01 eV was achieved within 3–9 iterations depending on the starting point, basis set size, and particular orbital considered. In this paper, we only present results from GnW0 and GnWn calculations with a PBE starting point, as we observed that GnW0 and GnWn schemes with PBE0 (or larger values of exact exchange) lead to worse agreement with experiment. The results for GnW0 with PBE0 and hybrid functional with α = 0.5 starting points are presented for completeness in the supplementary material.

Our GW calculations were performed using aug-cc-pVTZ and aug-ccpVQZ basis sets and extrapolated to the complete basis set (CBS) limit with the following function:77 

E=a+bNBF,
(2)

where E is the quasiparticle energy, a and b are the fitting parameters, and NBF is the number of basis functions. We tested the accuracy of these fits by performing the G0W0 calculations with various starting points for a closed-shell (ScO2) and an open-shell (TiO2) molecule with the aug-cc-pV5Z basis sets. We observed that including aug-cc-pV5Z basis sets only had a negligible effect (few tens of meVs) on the CBS limit. We note, however, that extrapolation to the CBS limit is a must for accurate predictions of quasiparticle energies in transition metal oxide molecular systems, as the CBS limit is typically lower than the value obtained with aug-cc-pVQZ basis sets by appreciable amounts ranging from ∼0.2 eV to 0.6 eV depending on the molecule and the particular orbital considered. These findings are consistent with the trends observed in Ref. 43. Finally, the Coulomb interaction terms were evaluated using the resolution-of-identity (RI) approximation.14,78 We checked the accuracy of the RI approximation at both the DFT and G0W0 levels, and the differences were found to be negligible, typically in the 1 meV–10 meV range.

The optimized structural parameters of the molecular anions considered in this study are shown in Table I. All molecules have the C2v symmetry. Our results for the TM–O bond lengths and the O–TM–O bond angles (obtained with the PBE functional) are in excellent agreement with the results of Gutsev et al.59 Also shown in Table I are the harmonic vibrational frequencies, which again agree very well with the results of Gutsev et al., and the ground state symmetries and spin multiplicities of the molecules. ScO2 is the only molecule with a closed-shell ground state; all others are open-shell ranging from the doublet for TiO2 to quintet ground state for MnO2. Next, we discuss the electronic structures and the photoelectron spectra obtained at various levels of DFT and GW theory for each of the five transition metal dioxide anions. For this discussion, we place the molecules in the yz plane and choose z as the C2 symmetry axis. We use the convention in which the b1 orbital is antisymmetric with respect to reflection in the yz plane and b2 is symmetric. We employ Mulliken population analysis to quantify the 3d characters of the molecular orbitals.

TABLE I.

Computed TM–O bond lengths, O–TM–O bond angles, vibrational frequencies, and ground state symmetries for the TMO2 molecules considered in this study.

BondBondGround
length (Å)angle (°)ω1 (cm−1)ω2 (cm−1)ω3 (cm−1)state
ScO2 1.81 124.1 184 769 671 1A1 
TiO2 1.682 111.4 321 924 888 2A1 
VO2 1.653 117.9 293 908 907 3B1 
CrO2 1.649 134.0 242 867 921 4B1 
MnO2 1.662 126.0 250 879 835 5B2 
BondBondGround
length (Å)angle (°)ω1 (cm−1)ω2 (cm−1)ω3 (cm−1)state
ScO2 1.81 124.1 184 769 671 1A1 
TiO2 1.682 111.4 321 924 888 2A1 
VO2 1.653 117.9 293 908 907 3B1 
CrO2 1.649 134.0 242 867 921 4B1 
MnO2 1.662 126.0 250 879 835 5B2 

Figure 1 shows the experimental photoelectron spectrum47 of ScO2 along with the spectra calculated at various levels of theory. The experimental spectrum obtained with a 266 nm (4.66 eV) laser consists of three main features, starting with a peak centered at 2.32 eV, followed by a more intense and broad peak near 2.9 eV (attributed in Ref. 47 to two overlapping features at 2.89 eV and 2.95 eV) and a higher energy broad peak centered at 3.68 eV. With the PBE (as well as PBE0) functional, we find the ground state of ScO2 to be a singlet with a valence configuration (8a1)2(5b2)2(1a2)2(3b1)2(9a1)2(6b2)2. We find the ground state of the neutral molecule to be a doublet with the electron removed from the 6b2 orbital. The vertical IPs are calculated as 2.10 eV and 2.17 eV with PBE and PBE0 functionals, respectively, which compare reasonably well with the experimental value of 2.32 eV. Both the shifted PBE and PBE0 spectra agree quite well with experimental data for the first three peaks, while they underestimate the fourth peak by 0.44 eV, which corresponds to the removal of 1a2 (or 5b2) electron with both functionals. In transition metal oxide molecules, shifted PBE spectra do not typically lead to good agreement with experimental data for localized orbitals or orbitals with considerable 3d character, as shown for the case of copper oxide molecular anions in Ref. 42. We attribute the apparent success of the shifted spectra of ScO2 computed with PBE for the first three peaks to the observation that the three highest occupied orbitals (3b1, 9a1, 6b2) are delocalized and primarily due to O 2px, 2py, and 2pz states, respectively, with negligible (<6%) Sc 3d character. The orbitals with the largest (∼20%) 3d content for this molecule are the 1a2 (with dxy character) and 5b2 (with dxz character) orbitals, for which the agreement with experimental data is not as satisfactory.47 Both the shifted PBE and shifted PBE0 spectra have the same mean absolute error (MAE) of 0.22 eV.

FIG. 1.

Experimental photoelectron spectrum of ScO2 (Ref. 47) along with spectra computed with shifted PBE, shifted PBE0, G0W0@PBE, GnW0@PBE, GnWn@PBE, G0W0@PBE0, and G0W0@PBEα = 0.5. Contour plots for some of the occupied molecular orbitals are shown on the right, with matching color codes displayed in the spectra.

FIG. 1.

Experimental photoelectron spectrum of ScO2 (Ref. 47) along with spectra computed with shifted PBE, shifted PBE0, G0W0@PBE, GnW0@PBE, GnWn@PBE, G0W0@PBE0, and G0W0@PBEα = 0.5. Contour plots for some of the occupied molecular orbitals are shown on the right, with matching color codes displayed in the spectra.

Close modal

Among the variants of the GW approximation, the predictions from the one-shot GW with a PBE starting point (G0W0@PBE) are particularly poor. All quasiparticle energies are underestimated significantly (by more than ∼1.1 eV) compared to experiment, leading to a large MAE of 1.24 eV at this level of theory. Predictions from G0W0 calculations with hybrid functional starting points are much better, as shown in Fig. 1. In fact, both G0W0@PBE0 and G0W0@PBEα = 0.5 levels of theory have the lowest MAE of 0.21 eV averaged over the four experimentally measured peaks. It is important, however, to note that not all peaks are predicted equally well at these two levels of theory: G0W0@PBE0 predictions for the first three peaks are quite accurate (within ∼0.15 eV), while the fourth peak is underestimated by 0.45 eV, very similar to the case of the shifted PBE0. G0W0@PBEα = 0.5 predictions, on the other hand, are not as good for the first three peaks (off by ∼0.3 eV), while the fourth peak is predicted within 0.01 eV. These findings correlate with the amount of 3d character of the relevant orbitals: Typically, orbitals with larger 3d content require a larger amount of exact exchange in the G0W0 starting point, as also observed for the case of copper oxide molecular anions.42 Another level of theory that leads to good agreement with experimental data is GnW0@PBE, which has been argued to be a practical GW scheme for the electronic structure of transition metal oxide molecular systems as a good compromise between computational efficiency and accuracy.43 The MAE for GnW0@PBE is 0.24 eV. Applying eigenvalue self-consistency also in W (GnWn), on the other hand, deteriorates the agreement with experimental data, and the MAE for this level of theory is 0.76 eV.

Figure 2 shows the experimental photoelectron spectrum of TiO2 and the spectra calculated with shifted DFT and various GW flavors. The most recent experimental spectrum obtained with a 193 nm (6.42 eV) laser consists of four main peaks at 1.60 eV, 3.90 eV, 4.72 eV, and 5.38 eV.51 With both PBE and PBE0 functionals, we find the ground state of TiO2 to be a doublet (2A1) with a valence configuration (1a2)2(3b1)2(9a1)2(6b2)2(10a1)1, and the ground state of the neutral molecule is a singlet with the electron removed from the very delocalized 10a1 orbital with large Ti 4s content and small Ti pz and 3dx2z2 admixtures. As in ScO2, the 6b2 state is entirely due to O 2pz states, but it is considerably more localized in TiO2. The higher binding energy (BE) orbitals 9a1 and 3b1 are also more localized in TiO2, and they have larger 3d content (∼15%) compared to ScO2.

FIG. 2.

Experimental photoelectron spectrum of TiO2 (Ref. 51) along with spectra computed with shifted PBE, shifted PBE0, G0W0@PBE, GnW0@PBE, GnWn@PBE, G0W0@PBE0, and G0W0@PBEα = 0.5. Contour plots for some of the occupied molecular orbitals, in both majority (↑) and minority (↓) spin channels, are shown on the right, with matching color codes displayed in the spectra. For doubly occupied orbitals, the calculated energy displayed with the higher (lower) height refers to the orbital in the majority (minority) spin channel.

FIG. 2.

Experimental photoelectron spectrum of TiO2 (Ref. 51) along with spectra computed with shifted PBE, shifted PBE0, G0W0@PBE, GnW0@PBE, GnWn@PBE, G0W0@PBE0, and G0W0@PBEα = 0.5. Contour plots for some of the occupied molecular orbitals, in both majority (↑) and minority (↓) spin channels, are shown on the right, with matching color codes displayed in the spectra. For doubly occupied orbitals, the calculated energy displayed with the higher (lower) height refers to the orbital in the majority (minority) spin channel.

Close modal

The vertical IPs calculated with PBE and PBE0 functionals are 1.54 eV and 1.66 eV, bracketing the experimental value of 1.60 eV. In the following comparisons, we compare the experimental values for the higher BE peaks with the average of the calculated spin-up and spin-down eigenvalues of a given doubly occupied orbital. While the vertical IPs calculated with both PBE and PBE0 are in very good agreement with experiment, shifted PBE clearly fails to provide satisfactory predictions when higher BE peaks are taken into account, as shown in Fig. 2. Shifted PBE0, on the other hand, performs quite well for the first three peaks and underestimates the energy of the fourth. This observation is another example of the well-known failure of semi-local functionals, such as PBE, as they underbind localized orbitals due to SIE, while the addition of a fraction of exact exchange globally, as in PBE0 or in range-separated form, helps mitigate this error.42 The overall MAE of shifted PBE and PBE0 is found to be 0.77 eV and 0.21 eV, respectively. These observations are in good agreement with the results of Marom et al.,68 who reported similar findings.

Among the variants of the GW approximation, G0W0@PBE has the worst performance with a MAE of 1.68 eV. Similar to ScO2, the best overall agreement with experiment is achieved with the G0W0@PBE0 level of theory (with a MAE of 0.16 eV), while adding more exact exchange worsens the agreement with experiment as the MAE of G0W0@PBEα = 0.5 is 0.51 eV, significantly larger than that of ScO2. The effect of eigenvalue self-consistency in predicting quasiparticle energies of TiO2 is similar to that of ScO2. GnW0 provides fairly good agreement with experiment (MAE of 0.25 eV), while iterating eigenvalues in W worsens the agreement by overbinding all orbitals other than the highest occupied molecular orbital (HOMO), leading to a MAE of 0.77 eV.

Previous studies on the electronic structure of VO2 have revealed many possible low energy configurations as candidates for its ground state.52,54,58,59,63,71 Based on an earlier study on the electronic structure of the neutral molecule,53 Wu and Wang interpreted their experimental photoelectron spectrum of VO2 by assuming that the ground state of the anion is a singlet of 1A1 symmetry.48 Later computational studies showed that the 1A1 state is higher in energy than the triplet states 3B1 and 3A1 by 0.3 eV–0.5 eV. However, depending on the electronic structure method and computational details, such as the basis set size and the choice of the exchange-correlation functional in DFT treatments, different studies have found the ground state to be either 3B1 or 3A1, typically by a very small energy difference. For example, the most recent DFT calculations of Kim et al.52 and quantum chemistry calculations of Hendrickx and Tran71 have found 3B1 as the ground state by small differences of 0.02 eV–0.07 eV relative to 3A1.

In our DFT calculations with the PBE functional, we find the ground state of VO2 to be the 3B1 triplet with the valence configuration (3b1)2(9a1)2(6b2)2(10a1)1(4b1)1. The ground state of the neutral molecule is found to be a doublet with the electron removed from the 4b1 orbital, which has a large (∼70%) 3dxz character and is significantly more localized than the 10a1 orbital. As a result, with the addition of a fraction of exact exchange, which reduces SIE, the 4b1 orbital moves down in energy much more than the 10a1 orbital. Accordingly, in our calculations with the PBE0 functional, while the ground state of VO2 is still found to be the 3B1 triplet, the HOMO is the 10a1 orbital. The downward shifts in the Kohn–Sham eigenvalues of the 4b1 and 10a1 orbitals in going from PBE to PBE0 are 1.83 eV and 0.90 eV, respectively, which leads to the observed reordering of the orbitals. We note that our PBE calculations find the 1A1 singlet state to be 0.4 eV higher in energy, in agreement with previous DFT and quantum chemistry calculations. Below, we only discuss the spectra obtained for the triplet configuration. As we show in the supplementary material, none of the theoretical spectra obtained for the singlet configuration resembles the experimental spectrum, which provides more evidence that the molecule sampled in the experiments is a triplet.

Figure 3 shows the experimental photoelectron spectrum of VO2 obtained with a 193 nm laser48 and our spectra calculated at various levels of theory. The experimental spectrum consists of four main peaks at 2.03 eV, 2.75 eV, 4.10 eV, and 4.85 eV. In our comparisons with the third and fourth experimental peaks, we take the average of our calculated spin-up and spin-down eigenvalues of doubly occupied orbitals. The vertical IPs calculated with PBE and PBE0 are 1.84 eV and 2.04 eV, the latter being in nearly perfect agreement with experiment. The other most obvious difference between shifted PBE and PBE0 spectra is the splitting of the first two peaks, which are calculated as 0.26 eV and 0.67 eV, respectively, compared to the experimental value of 0.72 eV. In fact, the whole shifted PBE0 spectrum has excellent agreement with experimental data with a MAE of 0.06 eV. The agreement with shifted PBE, on the other hand, is not as good with a MAE of 0.34 eV.

FIG. 3.

Experimental photoelectron spectrum of VO2 (Ref. 48) along with spectra computed with shifted PBE, shifted PBE0, G0W0@PBE, GnW0@PBE, GnWn@PBE, G0W0@PBE0, and G0W0@PBEα = 0.5. Contour plots for some of the occupied molecular orbitals, in both majority (↑) and minority (↓) spin channels, are shown on the right, with matching color codes displayed in the spectra. For doubly occupied orbitals, the calculated energy displayed with the higher (lower) height refers to the orbital in the majority (minority) spin channel.

FIG. 3.

Experimental photoelectron spectrum of VO2 (Ref. 48) along with spectra computed with shifted PBE, shifted PBE0, G0W0@PBE, GnW0@PBE, GnWn@PBE, G0W0@PBE0, and G0W0@PBEα = 0.5. Contour plots for some of the occupied molecular orbitals, in both majority (↑) and minority (↓) spin channels, are shown on the right, with matching color codes displayed in the spectra. For doubly occupied orbitals, the calculated energy displayed with the higher (lower) height refers to the orbital in the majority (minority) spin channel.

Close modal

The predictions from the one-shot GW approximation with PBE, PBE0, and PBEα = 0.5 starting points follow similar trends discussed earlier. G0W0@PBE predictions are very poor, while G0W0@PBE0 provides excellent agreement with experiment with a MAE of 0.09 eV, very similar to shifted PBE0. More exact exchange in the starting point worsens the agreement with experiment by overbinding all orbitals significantly, leading to a MAE of 0.58 eV at the G0W0@PBEα = 0.5 level of theory. Eigenvalue self-consistency in G with a PBE starting point interestingly does not provide as good agreement with experiment compared to other molecules. While the GnW0@PBE prediction for the IP is very good (1.98 eV), the second peak predicted at 2.00 eV is underestimated by 0.75 eV. Here, we should also note that the HOMO eigenvalue calculated with GnW0@PBE does not correspond to the 4b1 orbital (which is the HOMO at PBE and G0W0@PBE levels) but to the 10a1 orbital. In particular, even though GnW0 predicts the quasiparticle energies of the 10a1 orbital to be lower than those of the 4b1 orbital calculated with aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets, the eigenvalue difference gets smaller as the basis set size increases such that when the eigenvalues are extrapolated, the 10a1 orbital becomes the HOMO. This difference in the convergence behavior of 10a1 and 4b1 orbitals is related to the fact that the 4b1 orbital has a much higher 3d content, and increasing the basis set size for such localized orbitals typically leads to a much larger lowering of the corresponding quasiparticle energies.37 At the GnWn@PBE level, the predicted splitting between the 10a1 (now clearly the HOMO) and 4b1 (HOMO-1) orbitals significantly improves; however, the peaks with higher BE are overbound at this level of theory, leading to a MAE of 0.64 eV.

Figure 4 shows the experimental photoelectron spectrum of CrO2 obtained with a 193 nm laser50 along with the computed spectra. There are three main peaks in the experimental spectrum at energies 2.43 eV, 3.41 eV, and 4.25 eV. With both PBE and PBE0 functionals, we find the ground state of CrO2 to be the 4B1 quartet with the valence configuration (6b2)2(10a1)1(4b1)1(11a1)1. The ground state of the neutral molecule is the 3B1 triplet with the electron removed from the 11a1 orbital. Unlike the HOMO of VO2, which is a localized orbital with a large 3dxz character, the HOMO (11a1) of CrO2 is a delocalized orbital of primarily Ti 4s character with some Ti pz and 3dx2z2 admixtures, similar to the 10a1 orbitals in TiO2 and VO2. The frontier orbitals of large (∼70%) 3d character in CrO2 are HOMO-1 and HOMO-2, whose DFT eigenvalues undergo a large downward shift (∼1.8 eV) in going from PBE to PBE0 compared to the HOMO eigenvalue, which undergoes a shift of ∼1.2 eV. As a result, the 11a1 orbital remains as the HOMO with PBE0 also, and there is no orbital reordering as observed in VO2. We also note that in the computed photoelectron spectra, since the Kohn–Sham eigenvalues or the GW quasiparticle energies of the 10a1 and 4b1 states are quite close to each other, we assign their average to the measured peak at 3.41 eV (a similar assignment was done in Ref. 50), and the measured peak at 4.25 is compared with the average of the spin-up and spin-down 6b2 orbital energies.

FIG. 4.

Experimental photoelectron spectrum of CrO2 (Ref. 50) along with spectra computed with PBE, PBE0, G0W0@PBE, GnW0@PBE, GnWn@PBE, G0W0@PBE0, and G0W0@PBEα = 0.5. Contour plots for some of the occupied molecular orbitals, in both majority (↑) and minority (↓) spin channels, are shown on the right, with matching color codes displayed in the spectra. For a doubly occupied 6b2 orbital, the calculated energy displayed with the higher (lower) height refers to the orbital in the majority (minority) spin channel.

FIG. 4.

Experimental photoelectron spectrum of CrO2 (Ref. 50) along with spectra computed with PBE, PBE0, G0W0@PBE, GnW0@PBE, GnWn@PBE, G0W0@PBE0, and G0W0@PBEα = 0.5. Contour plots for some of the occupied molecular orbitals, in both majority (↑) and minority (↓) spin channels, are shown on the right, with matching color codes displayed in the spectra. For a doubly occupied 6b2 orbital, the calculated energy displayed with the higher (lower) height refers to the orbital in the majority (minority) spin channel.

Close modal

The vertical IPs calculated with PBE and PBE0 functionals are 2.30 eV and 2.48 eV, respectively, which are in reasonably good agreement with experiment. However, as shown in Fig. 4, the shifted PBE does not have a very good agreement with experiment overall (MAE of 0.36 eV) since the second peak is underestimated significantly. This is not surprising as this peak is due to 10a1 and 4b1 orbitals with large 3d character, which both suffer from SIE. As expected, the addition of a fraction of exact exchange mitigates this problem, and the shifted PBE0 spectrum has a much better agreement with experiment with a very low MAE of 0.06 eV. The GW predictions are overall similar to the general trends discussed earlier. One-shot GW with the PBE0 starting point has the best overall agreement with experiment with a MAE of 0.09 eV, and each predicted peak matches almost perfectly with the shifted PBE eigenvalues. Different from VO2 (but in line with other molecules), the GnW0@PBE predictions are quite good with a MAE of 0.18 eV. As before, adding self-consistency in W eigenvalues or one-shot GW with a larger exact exchange fraction worsens the agreement with experiment by overbinding the orbitals, and G0W0@PBE predicts orbitals that are too underbound, resulting in a very poor MAE of 1.66 eV.

In their combined experimental and theoretical study of the photoelectron spectra of MnOx molecules,49 Gutsev et al. found a strong temperature dependence of the spectrum for MnO2. The experimental data displayed in Fig. 5 show the spectrum obtained with a 266 nm laser under cold conditions. The small peak near 1.89 eV was attributed to a higher energy isomer and will not be considered in the following discussion. The vertical IP of MnO2 was found to be 2.26 eV based on the analysis of the vibrational progression between 2.0 eV and 2.7 eV, followed by a broad peak at 3.09 eV. This is followed by a congested spectrum with many overlapping features between 3.6 eV and 4.3 eV, from which we have identified peaks at 3.67 eV, 3.77 eV, 3.85 eV, 3.96 eV, and 4.11 eV, giving a total of seven peaks for comparison with our calculated spectra.

FIG. 5.

Experimental photoelectron spectrum of MnO2 (Ref. 49) along with spectra computed with PBE, PBE0, G0W0@PBE, GnW0@PBE, GnWn@PBE, G0W0@PBE0, and G0W0@PBEα = 0.5. Contour plots for some of the occupied molecular orbitals, in both majority (↑) and minority (↓) spin channels, are shown on the right, with matching color codes displayed in the spectra. For a doubly occupied 6b2 orbital, the calculated energy displayed with the higher (lower) height refers to the orbital in the majority (minority) spin channel.

FIG. 5.

Experimental photoelectron spectrum of MnO2 (Ref. 49) along with spectra computed with PBE, PBE0, G0W0@PBE, GnW0@PBE, GnWn@PBE, G0W0@PBE0, and G0W0@PBEα = 0.5. Contour plots for some of the occupied molecular orbitals, in both majority (↑) and minority (↓) spin channels, are shown on the right, with matching color codes displayed in the spectra. For a doubly occupied 6b2 orbital, the calculated energy displayed with the higher (lower) height refers to the orbital in the majority (minority) spin channel.

Close modal

With both PBE and PBE0 functionals, we find the ground state of MnO2 to be the 5B2 quintet with a valence configuration (3b1)2(6b2)2(10a1)1(4b1)1(11a1)1(2a2)1. We find the ground state of the neutral molecule to be the 4B1 quartet with the electron removed from the 2a2 orbital, which is a localized molecular orbital derived primarily from the Mn 3dxy atomic orbital hybridized with O 2px orbitals. Similar to the case of CrO2, 11a1 is rather delocalized with significant 4s character, while 4b1 and 10a1 are both localized and derived primarily from Mn 3dxz and 3dx2z2, respectively. In the following discussion, we compare the energies of five majority-spin orbital energies and two minority-spin orbital (6b2↓ and 3b1↓) energies with the seven peaks from the experimental spectrum.

The vertical IPs calculated with PBE and PBE0 functionals are 2.13 eV and 2.54 eV, respectively, compared to the experimental value of 2.26 eV. The deviation of the IP calculated with PBE0 from the experimental value is relatively large in comparison to other molecules considered here; however, the rest of the shifted PBE0 spectrum is in very good agreement with experiment, leading to a MAE of 0.19 eV. While the IP calculated with PBE is very close to the experimental value, overall the shifted PBE spectrum with a MAE of 0.30 eV does not have as good agreement with experiment as PBE0. The most obvious difference between the two spectra is the splitting between the first two peaks, which are found to be 1.19 eV and 0.61 eV with PBE and PBE0, respectively, compared to the experimental value of 0.83 eV. This can be understood in terms of the difference in the spatial extents of the relevant orbitals, 2a2 and 11a1, and the downward shift in energy they experience with the addition of a fraction of exact exchange, as discussed for several cases earlier. In the case of MnO2, the DFT eigenvalue of the more localized 2a2 orbital decreases by 1.87 eV, while that of the less localized 11a1 orbital decreases by 1.29 eV in going from PBE to PBE0 (keep in mind that the PBE and PBE0 spectra plotted in Fig. 5 are shifted spectra).

The predictions from the variants of the GW approximation are consistent with the overall trends discussed so far. Best agreement is achieved with G0W0@PBE0 (MAE of 0.12 eV), followed by GnW0@PBE with a MAE of 0.37 eV. While the first two peaks are predicted well with GnW0, states with higher BEs are overbound by 0.4 eV–0.5 eV. Adding eigenvalue self-consistency in W overbinds these orbitals more, and MAE increases to 0.97 eV. The PBE starting point for one-shot GW again leads to the worst agreement with experiment (MAE of 1.23 eV), and G0W0@PBEα = 0.5 overbinds all states by 0.3 eV–0.6 eV with a MAE of 0.43 eV.

The overall trends for how the predictions from various levels of DFT and GW theories compare with experimental data for all TMO2 anions are summarized in Fig. 6, which shows the differences between the computed and experimental vertical IPs (left panel) as well as the MAEs (right panel). We notice that the vertical IPs computed with the PBE functional are quite good, underestimated slightly (0.1 eV–0.2 eV) with respect to experiment. However, when higher BE peaks are taken into account, the agreement with experiment worsens, particularly for TiO2 with a MAE of 0.77 eV, which highlights the importance of comparing not just the vertical IP but also all the higher BE peaks in available PES in assessing the performance of a particular level of theory. For shifted PBE0, on the other hand, not only are the predictions for vertical IPs in reasonably good agreement with experiment (typically less than 0.15 eV, but with an outlier for the case of MnO2 where the deviation from experiment is 0.28 eV), but also the overall MAEs (less than ∼0.2 eV) are quite satisfactory.

FIG. 6.

Difference between experimental and our calculated vertical IPs (EbindexpEbindcal) (left) and mean absolute error (MAE) of calculated electron binding energies (right) of TMO2 molecules at different levels of theory.

FIG. 6.

Difference between experimental and our calculated vertical IPs (EbindexpEbindcal) (left) and mean absolute error (MAE) of calculated electron binding energies (right) of TMO2 molecules at different levels of theory.

Close modal

For the case of the one-shot G0W0@PBEα level of theory, the starting point makes a striking difference. As highlighted for all TMO2 anions in Sec. III, a PBE starting point leads to very poor agreement with experiment for both the vertical IPs (average underestimate compared to experiment for all cases considered being 1.36 eV) and the MAEs, which are slightly higher. The PBE0 starting point, on the other hand, leads to the best agreement overall with experimental data among all levels of theory considered in this study, not only for vertical IPs (within ∼0.1 eV) but also for higher BE peaks (with a similar MAE of 0.1 eV–0.2 eV). With an increasing amount of exact exchange in the DFT starting point, however, the agreement with experiment deteriorates: With G0W0@PBEα = 0.5, all orbitals are overbound by ∼0.5 eV, with the exception of ScO2 where MAE remains near 0.2 eV, close to the value obtained with α = 0.25. This is consistent with the results obtained for ScO in Ref. 43, where any value of α in the DFT starting point in the range 0.25 ≤ α ≤ 1 leads to good agreement with experimental data and does not shift the predicted quasiparticle energies significantly for this molecule, which was attributed to the lack of significant 3d character in the relevant orbitals. We also note that our finding of nearly excellent agreement of G0W0@PBEα at the “sweet spot” value of α = 0.25 should not be taken as a general trend for the G0W0@PBE0 level of theory for TMO molecular systems, as this work has focused particularly on early TM dioxide systems. Typically, the amount of α needed for good agreement with experimental data in the DFT starting point of one-shot GW calculations for TMO molecules is dependent on the amount of the 3d character of the orbitals43 as well as the content of TM in the molecule. For example, Shi et al.42 showed in their study of small copper oxide molecular anions that while α = 0.5 worked well for molecules such as Cu2O and CuO, values near α = 0.25 were optimal for molecules with less Cu content, such as CuO2 and CuO3.

For the case of evGW, updating eigenvalues only in G (GnW0) with a PBE starting point leads to very good agreement with experiment for the vertical IPs, typically within 0.05 eV, except for the two slight outliers of ScO2 and MnO2 where the deviation is ∼0.25 eV. When higher BE peaks are considered, the agreement with experiment slightly worsens, but the MAE value of 0.3 eV averaged over the five TMO2 anions is still reasonable. Adding self-consistency in W, on the other hand, overbinds all orbitals as the poles of the self-energy become more negative (for occupied orbitals) with the increase in the neutral excitation energies. For vertical IPs of TiO2, VO2, and CrO2, GnWn@PBE predictions are in good agreement with experiment (overestimated by 0.15 eV, 0.23 eV, and 0.40 eV, respectively), but the predicted IPs for ScO2 and MnO2 are not as they are overestimated by ∼0.9 eV. When all peaks are considered, GnWn@PBE is clearly not a predictive level of theory as the MAE averaged over the five TMO2 anions is ∼0.75 eV.

In order to gain more insight into the successes and failures of various levels of GW theories for the TMO2 anions considered, we now focus on the (averaged) eigenvalue shifts for each molecule, as shown in Fig. 7. For the GW levels of theory, the bars show for each TMO2 anion the magnitudes of the shifts of the GW eigenvalues from the corresponding DFT starting point eigenvalues averaged over the measured peaks. The standard deviations (SDs) of the shifts are shown with error bars for each level of GW theory. For PBE and PBE0, the bars show the magnitudes of the shifts needed to align the HOMO eigenvalue with negative of the computed vertical IP of the anion discussed earlier. One immediate observation from Fig. 7 is that the eigenvalue shifts for the G0W0@PBE0 (blue bars) are very close to the PBE0 shifts (cyan bars) for all molecules considered. This, combined with the observation that the SDs of the shifts for G0W0@PBE0 are quite small (less than 0.1 eV), shows that the G0W0@PBE0 spectra can be obtained by almost a uniform shift of the PBE0 spectra by the amount needed to align HOMO with negative of the vertical IP of the anion, irrespective of particular nature (localized, delocalized, large, or small 3d character) of the orbitals, resulting also in similar MAEs with respect to experimental data. The very good agreement of the G0W0@PBE0 level of theory with experimental PES can then be attributed to the correct energetic positioning of the orbitals obtained with the PBE0 and PBE0 wave functions being good approximations to the true quasiparticle wave functions so that a simple first-order G0W0 perturbative correction is enough to lead to accurate quasiparticle energies.

FIG. 7.

Calculated eigenvalue shifts for TMO2 molecules at each level of theory. For PBE and PBE0, the bars show the magnitudes of shifts needed to align the HOMO eigenvalues with negative of the computed vertical IPs. For GW levels of theory, the bars show the magnitudes of the shifts of the GW eigenvalues from the eigenvalues of the corresponding DFT starting points averaged over the measured peaks. The error bars are the standard deviations of the shifts for each level of GW theory. The black dots show the shifts for the delocalized 10a1 and 11a1 orbitals discussed in the text. This orbital is unoccupied in ScO2.

FIG. 7.

Calculated eigenvalue shifts for TMO2 molecules at each level of theory. For PBE and PBE0, the bars show the magnitudes of shifts needed to align the HOMO eigenvalues with negative of the computed vertical IPs. For GW levels of theory, the bars show the magnitudes of the shifts of the GW eigenvalues from the eigenvalues of the corresponding DFT starting points averaged over the measured peaks. The error bars are the standard deviations of the shifts for each level of GW theory. The black dots show the shifts for the delocalized 10a1 and 11a1 orbitals discussed in the text. This orbital is unoccupied in ScO2.

Close modal

With a PBE starting point, on the other hand, the average shifts for G0W0@PBE (yellow bars in Fig. 7) are much smaller than the PBE eigenvalue shifts required to align HOMO with negative of the vertical IP (red bars in Fig. 7). That is, unlike the very good agreement of shifted PBE0 and G0W0@PBE0 with each other (as well as with experiment), application of first-order G0W0 perturbative correction to the PBE eigenvalues using PBE wave functions fails to provide the required amount of downward shift to achieve a reasonably good level of agreement with experimental data. To quantify this further, we compared the PBE and PBE0 eigenvalues for all TMO2 anions and considered the amounts of shifts that would be needed to achieve perfect agreement with experimental data. Averaged over all experimental peaks and all TMO2 anions, the required shift is 3.09 ± 0.27 eV and 1.79 ± 0.13 eV for PBE and PBE0, respectively, showing (i) how underbound PBE orbitals are for perturbative corrections to be effective and (ii) that the required shifts for PBE have a much larger variation than those for PBE0 due to the significant difference in the SDs.

These last two observations are also consistent with the trends of eigenvalue shifts for GnW0@PBE shown with green bars in Fig. 7. In this case, the shifts are significantly larger, correcting the eigenvalues of the underbound PBE orbitals, and the amount of shift depends on the nature of the particular orbital due to the large SD observed for most TMO2 molecular anions, with the exception of ScO2, where the occupied frontier orbitals are composed of primarily O 2p states and have similar spatial extents. The large SD observed in GnW0@PBE eigenvalue shifts is primarily due to a particular orbital with a1 character (10a1 in TiO2 and VO2, and 11a1 in CrO2 and MnO2), which is either HOMO or HOMO-1 in the majority spin channel depending on the level of theory (see Fig. 3 for VO2) or the molecule (in MnO2, HOMO is the 2a2 orbital). The shifts for this particular orbital also displayed in Fig. 7 are significantly less than the averaged shift, e.g., the GnW0@PBE eigenvalue for the 10a1 state of TiO2 is 2.2 eV lower than the PBE eigenvalue, which is much smaller than the average downward shift of 3.27 eV. As a result, the SD is particularly high at the GnW0@PBE level of theory as different orbitals are shifted by varying amounts to achieve a reasonably good level of agreement with experiment. As mentioned earlier, this particular a1 orbital has a large TM 4s character with some TM pz and 3dx2z2 admixtures. It is important to note that while the 3d character is not very large, it is still appreciable, ranging from 15% to 30% in going from TiO2 to MnO2. The primary reason why this orbital does not undergo a large downward shift is its very delocalized nature despite having an appreciable 3d content, consistent with the observation that GnW0@PBE does a particularly good job for the vertical IPs, and the agreement deteriorates slightly for the more localized higher BE peaks.

In summary, we provided a detailed comparison of the DFT- and GW-computed eigenvalue spectra to experimental photoelectron spectra of five 3d-transition metal dioxide molecular anions, ScO2, TiO2, VO2, CrO2, and MnO2. We focused our study on the comparison between semi-local and hybrid functional predictions (within DFT and as a starting point for G0W0 calculations) as well as the effects of eigenvalue self-consistency with a semi-local DFT starting point. Overall, both the shifted PBE0 and G0W0@PBE0 appear to stand out as the best levels of theory as they are able to reproduce the experimental BEs with 0.1 eV–0.2 eV accuracy (for some peaks, especially the vertical IPs, with much better accuracy). With shifted PBE, on the other hand, while the total energy differences between the anion and the neutral molecule at the geometry of the anion lead to similarly accurate vertical IPs, the agreement with experiment quickly deteriorates for higher BE peaks. A one-shot GW calculation with a semi-local starting point such as PBE is clearly a poor choice leading to the worst agreement with experiment (average error ∼1.5 eV or more) as the PBE orbitals are too underbound and the PBE wave functions are too delocalized for a perturbative self-energy correction to be effective. These observations highlight the importance of adding a fraction of Fock exchange to mitigate self-interaction errors in predicting the quasiparticle energies of transition metal oxide molecular systems within DFT or GW levels of theory. However, too much Fock exchange (e.g., α = 0.5) also worsens the agreement with experiment by overbinding orbitals. We found that G-only eigenvalue self-consistent GW with W computed from PBE (GnW0@PBE), while not as accurate as shifted PBE0 or G0W0@PBE0, still provides a reasonably good description of the quasiparticle energies for occupied frontier orbitals in these transition metal oxide molecular systems. Updating the eigenvalues in W (GnWn@PBE), on the other hand, leads to poor agreement with experiment overall, even though it may be reasonably accurate for vertical IPs of some of the molecular anions considered. We caution the reader that our finding of very good agreement of G0W0@PBE0 predictions with experimental data for the particular case of these early 3d-transition metal dioxide molecular anions should be viewed in the context of our previous studies on other transition metal oxide molecular systems,42,43 which showed that higher values of Fock exchange might be needed in the DFT starting point of GW calculations to achieve good agreement with experiment in systems with more localized orbitals and higher transition metal content. When viewed from this perspective, the results presented in this paper still lead us to advocate using GnW0@PBE as a practical GW scheme that presents a good balance between accuracy and efficiency in predicting quasiparticle energies of transition metal oxide molecular systems without the need for a system-dependent parameter in the starting DFT description. Future studies on the more challenging late 3d-transition metal dioxide molecular anions, such as FeO2, NiO2, and CuO2, and other transition metal oxide molecules with varying amounts of transition metal content are needed to further assess and understand the predictive capabilities of GW methods for electronic excitations in molecular systems with moderate electron correlation.

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

This work was supported by the U.S. Department of Energy under Grant No. DE-SC0017824. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We would also like to thank Young-Moo Byun for useful discussions in the earlier stages of this work.

The data that support the findings of this study are available within this article and its supplementary material.

1.
L.
Reining
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1344
(
2018
).
2.
P.
Rinke
,
A.
Qteish
,
J.
Neugebauer
,
C.
Freysoldt
, and
M.
Scheffler
,
New J. Phys.
7
,
126
(
2005
).
3.
F.
Bruneval
,
N.
Vast
, and
L.
Reining
,
Phys. Rev. B
74
,
045102
(
2006
).
4.
M. L.
Tiago
and
J. R.
Chelikowsky
,
Phys. Rev. B
73
,
205334
(
2006
).
5.
M.
Shishkin
and
G.
Kresse
,
Phys. Rev. B
75
,
235102
(
2007
).
6.
F.
Fuchs
,
J.
Furthmüller
,
F.
Bechstedt
,
M.
Shishkin
, and
G.
Kresse
,
Phys. Rev. B
76
,
115109
(
2007
).
7.
C.
Rostgaard
,
K. W.
Jacobsen
, and
K. S.
Thygesen
,
Phys. Rev. B
81
,
085103
(
2010
).
8.
Y.
Ma
,
M.
Rohlfing
, and
C.
Molteni
,
J. Chem. Theory Comput.
6
,
257
265
(
2010
).
9.
X.
Blase
and
C.
Attaccalite
,
Appl. Phys. Lett.
99
,
171909
(
2011
).
10.
X.
Blase
,
C.
Attaccalite
, and
V.
Olevano
,
Phys. Rev. B
83
,
115103
(
2011
).
11.
C.
Faber
,
C.
Attaccalite
,
V.
Olevano
,
E.
Runge
, and
X.
Blase
,
Phys. Rev. B
83
,
115123
(
2011
).
12.
X.
Qian
,
P.
Umari
, and
N.
Marzari
,
Phys. Rev. B
84
,
075103
(
2011
).
13.
14.
X.
Ren
,
P.
Rinke
,
V.
Blum
,
J.
Wieferink
,
A.
Tkatchenko
,
A.
Sanfilippo
,
K.
Reuter
, and
M.
Scheffler
,
New J. Phys.
14
,
053020
(
2012
).
15.
N.
Marom
,
F.
Caruso
,
X.
Ren
,
O. T.
Hofmann
,
T.
Körzdörfer
,
J. R.
Chelikowsky
,
A.
Rubio
,
M.
Scheffler
, and
P.
Rinke
,
Phys. Rev. B
86
,
245127
(
2012
).
16.
B.
Baumeier
,
D.
Andrienko
,
Y.
Ma
, and
M.
Rohlfing
,
J. Chem. Theory Comput.
8
,
997
1002
(
2012
).
17.
M. J.
van Setten
,
F.
Weigend
, and
F.
Evers
,
J. Chem. Theory Comput.
9
,
232
246
(
2013
).
18.
F.
Bruneval
and
M. A. L.
Marques
,
J. Chem. Theory Comput.
9
,
324
329
(
2013
).
19.
T. A.
Pham
,
H.-V.
Nguyen
,
D.
Rocca
, and
G.
Galli
,
Phys. Rev. B
87
,
155148
(
2013
).
20.
F.
Caruso
,
P.
Rinke
,
X.
Ren
,
A.
Rubio
, and
M.
Scheffler
,
Phys. Rev. B
88
,
075105
(
2013
).
21.
V.
Atalla
,
M.
Yoon
,
F.
Caruso
,
P.
Rinke
, and
M.
Scheffler
,
Phys. Rev. B
88
,
165122
(
2013
).
22.
J.
Klimeš
,
M.
Kaltak
, and
G.
Kresse
,
Phys. Rev. B
90
,
075125
(
2014
).
23.
C.
Faber
,
P.
Boulanger
,
C.
Attaccalite
,
I.
Duchemin
, and
X.
Blase
,
Philos. Trans. R. Soc., A
372
,
20130271
(
2014
).
24.
P.
Koval
,
D.
Foerster
, and
D.
Sánchez-Portal
,
Phys. Rev. B
89
,
155417
(
2014
).
25.
P.
Boulanger
,
D.
Jacquemin
,
I.
Duchemin
, and
X.
Blase
,
J. Chem. Theory Comput.
10
,
1212
1218
(
2014
).
26.
S.
Körbel
,
P.
Boulanger
,
I.
Duchemin
,
X.
Blase
,
M. A. L.
Marques
, and
S.
Botti
,
J. Chem. Theory Comput.
10
,
3934
3943
(
2014
).
27.
L.-W.
Wang
,
Phys. Rev. B
91
,
125135
(
2015
).
28.
D.
Hirose
,
Y.
Noguchi
, and
O.
Sugino
,
Phys. Rev. B
91
,
205111
(
2015
).
29.
F.
Bruneval
,
S. M.
Hamed
, and
J. B.
Neaton
,
J. Chem. Phys.
142
,
244101
(
2015
).
30.
D.
Jacquemin
,
I.
Duchemin
, and
X.
Blase
,
J. Chem. Theory Comput.
11
,
3290
3304
(
2015
).
31.
M. J.
van Setten
,
F.
Caruso
,
S.
Sharifzadeh
,
X.
Ren
,
M.
Scheffler
,
F.
Liu
,
J.
Lischner
,
L.
Lin
,
J. R.
Deslippe
,
S. G.
Louie
,
C.
Yang
,
F.
Weigend
,
J. B.
Neaton
,
F.
Evers
, and
P.
Rinke
,
J. Chem. Theory Comput.
11
,
5665
5687
(
2015
).
32.
J.
Wilhelm
,
M.
Del Ben
, and
J.
Hutter
,
J. Chem. Theory Comput.
12
,
3623
3635
(
2016
).
33.
F.
Kaplan
,
M. E.
Harding
,
C.
Seiler
,
F.
Weigend
,
F.
Evers
, and
M. J.
van Setten
,
J. Chem. Theory Comput.
12
,
2528
2541
(
2016
).
34.
F.
Caruso
,
M.
Dauth
,
M. J.
van Setten
, and
P.
Rinke
,
J. Chem. Theory Comput.
12
,
5076
5087
(
2016
).
35.
J. W.
Knight
,
X.
Wang
,
L.
Gallandi
,
O.
Dolgounitcheva
,
X.
Ren
,
J. V.
Ortiz
,
P.
Rinke
,
T.
Körzdörfer
, and
N.
Marom
,
J. Chem. Theory Comput.
12
,
615
626
(
2016
).
36.
L.
Hung
,
F. H.
da Jornada
,
J.
Souto-Casares
,
J. R.
Chelikowsky
,
S. G.
Louie
, and
S.
Öğüt
,
Phys. Rev. B
94
,
085125
(
2016
).
37.
L.
Hung
,
F.
Bruneval
,
K.
Baishya
, and
S.
Öğüt
,
J. Chem. Theory Comput.
13
,
2135
2146
(
2017
).
38.
L.
Hung
,
F.
Bruneval
,
K.
Baishya
, and
S.
Öğüt
,
J. Chem. Theory Comput.
13
,
5820
5821
(
2017
).
39.
E.
Maggio
,
P.
Liu
,
M. J.
van Setten
, and
G.
Kresse
,
J. Chem. Theory Comput.
13
,
635
648
(
2017
).
40.
M.
Govoni
and
G.
Galli
,
J. Chem. Theory Comput.
14
,
1895
1909
(
2018
).
41.
M.
Grumet
,
P.
Liu
,
M.
Kaltak
,
J.
Klimeš
,
G.
Kresse
,
M.
Kaltak
, and
G.
Kresse
,
Phys. Rev. B
98
,
155143
(
2018
).
42.
B.
Shi
,
S.
Weissman
,
F.
Bruneval
,
L.
Kronik
, and
S.
Öğüt
,
J. Chem. Phys.
149
,
064306
(
2018
).
43.
Y.-M.
Byun
and
S.
Öğüt
,
J. Chem. Phys.
151
,
134305
(
2019
).
44.
W.
Gao
and
J. R.
Chelikowsky
,
J. Chem. Theory Comput.
15
,
5299
5307
(
2019
).
45.
D.
Golze
,
M.
Dvorak
, and
P.
Rinke
,
Front. Chem.
7
,
377
(
2019
).
46.
H.
Wu
and
L. S.
Wang
,
J. Chem. Phys.
107
,
8221
(
1997
).
47.
H.
Wu
and
L.-S.
Wang
,
J. Phys. Chem. A
102
,
9129
9135
(
1998
).
48.
H.
Wu
and
L.-S.
Wang
,
J. Chem. Phys.
108
,
5310
(
1998
).
49.
G. L.
Gutsev
,
B. K.
Rao
,
P.
Jena
,
X.
Li
, and
L.-S.
Wang
,
J. Chem. Phys.
113
,
1473
(
2000
).
50.
G. L.
Gutsev
,
P.
Jena
,
H.-J.
Zhai
, and
L.-S.
Wang
,
J. Chem. Phys.
115
,
7935
(
2001
).
51.
H.-J.
Zhai
and
L.-S.
Wang
,
J. Am. Chem. Soc.
129
,
3022
3026
(
2007
).
52.
J. B.
Kim
,
M. L.
Weichman
, and
D. M.
Neumark
,
J. Chem. Phys.
140
,
034307
(
2014
).
53.
L. B.
Knight
,
R.
Babb
,
M.
Ray
,
T. J.
Banisaukas
,
L.
Russon
,
R. S.
Dailey
, and
E. R.
Davidson
,
J. Chem. Phys.
105
,
10237
(
1996
).
54.
G. V.
Chertihin
,
W. D.
Bare
, and
L.
Andrews
,
J. Phys. Chem. A
101
,
5090
5096
(
1997
).
55.
M. B.
Walsh
,
R. A.
King
, and
H. F.
Schaefer
,
J. Chem. Phys.
110
,
5224
(
1999
).
56.
M.
Zhou
and
L.
Andrews
,
J. Chem. Phys.
111
,
4230
(
1999
).
57.
J. M.
Gonzales
,
R. A.
King
, and
H. F.
Schaefer
,
J. Chem. Phys.
113
,
567
(
2000
).
58.
S. F.
Vyboishchikov
and
J.
Sauer
,
J. Phys. Chem. A
104
,
10913
10922
(
2000
).
59.
G. L.
Gutsev
,
B. K.
Rao
, and
P.
Jena
,
J. Phys. Chem. A
104
,
11961
11971
(
2000
).
60.
J.
Dong
,
Y.
Wang
, and
M.
Zhou
,
Chem. Phys. Lett.
364
,
511
516
(
2002
).
61.
F.
Grein
,
J. Chem. Phys.
126
,
034313
(
2007
).
62.
S.
Li
and
D. A.
Dixon
,
J. Phys. Chem. A
112
,
6646
6666
(
2008
).
63.
E. L.
Uzunova
,
H.
Mikosch
, and
G. S.
Nikolov
,
J. Chem. Phys.
128
,
094307
(
2008
).
64.
Y.
Gong
,
M.
Zhou
, and
L.
Andrews
,
Chem. Rev.
109
,
6765
6808
(
2009
).
65.
Y.
Liu
,
Y.
Yuan
,
Z.
Wang
,
K.
Deng
,
C.
Xiao
, and
Q.
Li
,
J. Chem. Phys.
130
,
174308
(
2009
).
66.
E. P. F.
Lee
,
D. K. W.
Mok
,
F.-T.
Chau
, and
J. M.
Dyke
,
J. Comput. Chem.
30
,
337
345
(
2009
).
67.
Z.-W.
Qu
and
H.
Zhu
,
J. Comput. Chem.
31
,
2038
2045
(
2010
).
68.
N.
Marom
,
J. E.
Moussa
,
X.
Ren
,
A.
Tkatchenko
, and
J. R.
Chelikowsky
,
Phys. Rev. B
84
,
245115
(
2011
).
69.
D. J.
Taylor
and
M. J.
Paterson
,
J. Chem. Phys.
408
,
1
10
(
2012
).
70.
J. B.
Kim
,
M. L.
Weichman
, and
D. M.
Neumark
,
Phys. Chem. Chem. Phys.
15
,
20973
(
2013
).
71.
M. F. A.
Hendrickx
and
V. T.
Tran
,
J. Chem. Theory Comput.
10
,
4037
4044
(
2014
).
72.
G.
Xu
and
F.
Yu
,
J. Electron Spectrosc. Relat. Phenom.
205
,
10
22
(
2015
).
73.
M.
Valiev
,
E. J.
Bylaska
,
N.
Govind
,
K.
Kowalski
,
T. P.
Straatsma
,
H. J. J.
Van Dam
,
D.
Wang
,
J.
Nieplocha
,
E.
Apra
,
T. L.
Windus
, and
W. A.
de Jong
,
Comput. Phys. Commun.
181
,
1477
(
2010
).
74.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
75.
J. P.
Perdew
,
M.
Ernzerhof
, and
K.
Burke
,
J. Chem. Phys.
105
,
9982
(
1996
).
76.
F.
Bruneval
,
T.
Rangel
,
S. M.
Hamed
,
M.
Shao
,
C.
Yang
, and
J. B.
Neaton
,
Comput. Phys. Commun.
208
,
149
(
2016
).
77.
C.
Hättig
,
W.
Klopper
,
A.
Köhn
, and
D. P.
Tew
,
Chem. Rev.
112
,
4
74
(
2012
).
78.
J. G.
Hill
and
J. A.
Platts
,
J. Chem. Phys.
128
,
044104
(
2008
).

Supplementary Material