Photoelectron spectra of early 3d-transition metal dioxide cluster anions from GW calculations
PPhotoelectron spectra of early d -transition metal dioxide cluster anionsfrom GW calculations Meisam Rezaei and Serdar Öğüt Department of Physics, University of Illinois at Chicago, Chicago, IL 60607, USA a) (Dated: 6 January 2021) Photoelectron spectra of early 3 d − transition metal dioxide anions, ScO − , TiO − , VO − , CrO − , MnO − , are calculatedusing semilocal and hybrid density functional theory (DFT) and many-body perturbation theory within the GW approx-imation using one-shot perturbative and eigenvalue self-consistent formalisms. Different levels of theory are comparedwith each other and with available photoelectron spectra. We show that one-shot GW with a PBE0 starting point( G W @PBE0) consistently provides very good agreement for all experimentally measured binding energies (within0.1-0.2 eV or less), which we attribute to the success of PBE0 in mitigating self-interaction error and providing goodquasiparticle wave functions, which renders a first-order perturbative GW correction effective. One-shot GW calcula-tions with semilocal exchange in the DFT starting point ( e.g. G W @PBE) do poorly in predicting electron removalenergies by underbinding orbitals with typical errors near 1.5 eV. Higher amounts of exact exchange ( e.g. GW do not provide very good agreement with experiment by overbinding orbitals withtypical errors near 0.5 eV. While not as accurate as G W @PBE0, the G -only eigenvalue self-consistent GW schemewith W fixed to the PBE level ( G n W @PBE) provides a reasonably predictive level of theory (typical errors near 0.3eV) to describe photoelectron spectra of these 3 d − transition metal dioxide anions. Adding eigenvalue self-consistencyalso in W ( G n W n @PBE), on the other hand, worsens the agreement with experiment overall. Our findings on the per-formance of various GW methods are discussed in the context of our previous studies on other transition metal oxidemolecular systems. I. INTRODUCTION
The last decade has witnessed a growing number of com-putational studies that have benchmarked Green’s functionmethods, such as the GW approximation and the Bethe-Salpeter equation, for excited state properties of bulk andmolecular systems. A large majority of these studies havefocused on the performance of various flavors of the GW ap-proximation in predicting the electron removal energies in sp -bonded molecules and clusters by comparing their predic-tions with accurate quantum chemistry calculations and ex-perimental photoelectron spectroscopy data. Modeling ex-cited states of systems containing transition metal elementsin general, and transition metal oxides in particular, withinthe GW theory have faced additional theoretical and compu-tational challenges, as enhanced electron correlations inher-ent in these systems and their propensity to having open-shellelectronic configurations necessitate the use of more sophis-ticated approaches beyond simple perturbative implementa-tions on top of density functional theory (DFT) with semi-local exchange-correlation functionals. Furthermore, conver-gence issues with respect to basis set size and other imple-mentation parameters present significant computational bot-tlenecks for maintaining a desirable level of accuracy compa-rable to what can be achieved for sp -bonded systems that istypically ∼ . GW approximation for sys-tems that contain transition metal elements. Motivated bythis observation, here we continue with our recent benchmarkstudies by focusing on the electronic structure of neg- a) Electronic mail: [email protected]. atively charged 3 d -transition metal dioxide clusters TMO − ,for TM = Sc, Ti, V, Cr, and Mn.Photoelectron spectra of early 3 d TMO − clusters have beenavailable up to photon energies of 5 − . Their structural and electronic propertieshave been investigated in various computational studiesusing methods based on DFT and quantum chemistry. Whilethese TMO − molecules are isostructural with little changes inthe bond angles and lengths upon changing the TM element,their frontier molecular orbitals display a wide range of spatiallocalization properties and have significantly varying amountsof TM 3 d and O 2 p contents. Therefore, their electronic struc-tures provide a diverse set of challenges to DFT within thegeneralized Kohn-Sham scheme and to many-body perturba-tion theory within the GW approximation, due to the need tomitigate self-interaction error (SIE) and the importance of asuitable DFT starting point or a form of self-consistency intheir GW treatment. In this work, we model the photoelectronspectra (PES) of early TMO − clusters using DFT with semi-local and hybrid exchange-correlation functionals (based onshifted eigenvalue spectra) and within the GW approximationusing one-shot perturbative schemes with different DFT start-ing points as well as eigenvalue self-consistent formalisms.This article is organized as follows. We begin with anoverview of the computational methods and parameters usedin our calculations in Sec. II. This is followed in Sec. IIIwith a brief discussion of the structural and vibrational prop-erties of the five TMO − clusters considered in this study anddetailed analyses of their PES computed with DFT and GW methods in comparison to each other and to experimental data.The overall trends in the performance of shifted DFT and vari-ants of the GW approximation are discussed in Sec. IV. Fi-nally, we summarize our findings and analyses in Sec. V. a r X i v : . [ c ond - m a t . m t r l - s c i ] J a n II. COMPUTATIONAL DETAILS
DFT computations to find the optimized structures werecarried out with the NWCHEM code, Version 6.6, usingPBE exchange-correlation functional and aug-cc-pVTZ ba-sis sets. The photoelectron spectra were simulated by broad-ening the eigenvalue spectra with a Gaussian distributionfunction of 0.1 eV smearing width, without taking into ac-count the photoionization cross sections. For DFT spectraobtained with PBE and PBE0 functionals, the Kohn-Shameigenvalues were shifted to align the first peak with the verti-cal ionization potential (IP) of the anion, which is computedas the total energy difference between the anionic and neutralclusters at the fixed anionic geometry. GW calculations were carried out using the MOLGW software within both the one-shot (perturbative) G W schemeand the eigenvalue self-consistent GW (ev GW ) scheme withtwo types of self-consistency, G n W and G n W n , which updatethe 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 E PBE α xc = α E HFx + ( − α ) E PBEx + E PBEc , (1)where E HFx , E PBEx , E PBEc and α are the Fock exact exchange,PBE exchange, PBE correlation energies, and the amount ofexact exchange, respectively. For starting point dependencyof the G W scheme, we tested three different values for α :0 (PBE starting point), 0.25 (PBE0 starting point) and 0.50.In these full-frequency G W calculations, we solved the non-linear quasiparticle equation for each state graphically usingthe secant (quasi-Newton) method. The numerical parame-ter η used to broaden the self-energy poles was chosen as η = .
001 Hartree, and the self-energy was evaluated on afrequency grid of spacing ∆ ω = .
001 Hartree. As discussedin detail in Ref. 43, full frequency G W method leads tocomplicated self-energy pole structures, typically at nonfron-tier orbitals, which results in multiple solutions of the quasi-particle equation and makes determining the correct and ac-curate quasiparticle energies difficult. This multisolution is-sue of G W is particularly prevalent in transition metal oxidesystems with a PBE starting point. As such, we checked theaccuracy of our predictions for the quasiparticle energies by(i) also using the spectral-function method, which locates thequasiparticle peak with the highest weight in the spectral func-tion, (ii) looking at the trends as a function of the basis set sizeand the amount of exact exchange, and (iii) by varying the pa-rameters η and ∆ ω , as recommended in Ref. 43.For calculations performed with the G n W and G n W n self-consistency methods, we used the iterative scheme employedin MOLGW, described in Ref. 43, where the quasiparticlerenormalization factor Z is set to Z =
1. In these calculations,we used a larger broadening parameter η = .
01 Hartree forthe self-energy poles (and accordingly, a larger frequency gridspacing of ∆ ω = .
01 Hartree) in order to avoid the oscilla-tions typically observed in the quasiparticle energies of non-frontier orbitals as a function of the iteration index. With theseparameters, self-consistency to 0.01 eV was achieved within3-9 iterations depending on the starting point, basis set size, and the particular orbital considered. In this paper, we onlypresent results from G n W and G n W n calculations with a PBEstarting point, as we observed that G n W and G n W n schemeswith PBE0 (or larger values of exact exchange) lead to worseagreement with experiment. The results for G n W with PBE0and hybrid functional with α = . GW calculations were performed using aug-cc-pVTZand aug-ccpVQZ basis sets and extrapolated to the completebasis set (CBS) limit with the following function E = a + bN BF , (2)where E is a quasiparticle energy, a and b are the fitting pa-rameters, and N BF is the number of basis functions. We testedthe accuracy of these fits by performing the G W calcula-tions with various starting points for a closed-shell (ScO − )and open-shell (TiO − ) molecule with the aug-cc-pV5Z ba-sis sets. We observed that including aug-cc-pV5Z basis setsonly had a negligible effect (few tens of meVs) on the CBSlimit. We note, however, that extrapolation to the CBS limitis a must for accurate predictions of quasiparticle energies intransition metal oxide molecular systems, as the CBS limit istypically lower than the value obtained with aug-cc-pVQZ ba-sis sets by appreciable amounts ranging from approximately0.2 eV to 0.6 eV depending on the molecule and the partic-ular orbital considered. These findings are consistent withtrends observed in Ref. 43. Finally, the Coulomb interac-tion terms were evaluated using the resolution-of-identity (RI)approximation. We checked the accuracy of the RI ap-proximation at both the DFT and G W levels, and the differ-ences were found to be negligible, typically in the 1 −
10 meVrange.
III. RESULTS
The optimized structural parameters of the molecular an-ions considered in this study are shown in Table I. Allmolecules have the C v symmetry. Our results for the TM-O bond lengths and the O-TM-O bond angles (obtained withPBE functional) are in excellent agreement with the resultsof Gutsev et al. Also shown in Table I are the harmonicvibrational frequencies, which again agree very well withthe results of Gutsev et al. , and the ground state symmetriesand spin multiplicities of the molecules. ScO − is the onlymolecule with a closed-shell ground state; all others are open-shell ranging from doublet for TiO − to quintet ground statefor MnO − . Next, we discuss the electronic structures and thephotoelectron spectra obtained at various levels of DFT and GW theory for each of the five transition metal dioxide an-ions. For this discussion, we place the molecules in the yz plane and choose z as the C symmetry axis. We use the con-vention in which the b orbital is antisymmetric with respectto reflection in the yz plane and b is symmetric. We employMulliken population analysis to quantify the 3 d characters ofthe molecular orbitals. TABLE I. Computed TM-O bond lengths, O-TM-O bond angles, vibrational frequencies, and ground state symmetries for the TMO − molecules considered in this study.Bond length (Å) Bond angle ( ◦ ) ω ( a ) ω ( a ) ω ( b ) Ground stateScO − A TiO − A VO − B CrO − B MnO − B A. ScO − Figure 1 shows the experimental photoelectron spectrum of ScO − along with the spectra calculated at various lev-els of theory. The experimental spectrum obtained witha 266 nm (4.66 eV) laser consists of three main features,starting with a peak centered at 2.32 eV, followed by amore intense and broad peak near 2.9 eV (attributed inRef. 47 to two overlapping features at 2.89 and 2.95 eV)and a higher energy broad peak centered at 3.68 eV. WithPBE (as well as PBE0) functional, we find the groundstate of ScO − to be a singlet with a valence configuration ( a ) ( b ) ( a ) ( b ) ( a ) ( b ) . We find the groundstate of the neutral molecule to be a doublet with the electronremoved from the 6 b orbital. The vertical IPs are calculatedas 2.10 and 2.17 eV with PBE and PBE0 functionals, respec-tively, which compare reasonably well with the experimentalvalue of 2.32 eV. Both the shifted PBE and PBE0 spectra agreequite well with experimental data for the first three peaks,while they underestimate the fourth peak by 0.44 eV, whichcorresponds to the removal of 1 a (or 5 b ) electron with bothfunctionals. In transition metal oxide molecules, shifted PBEspectra do not typically lead to good agreement with experi-mental data for localized orbitals or orbitals with considerable3 d character, as shown for the case of copper oxide molecu-lar anions in Ref. 42. We attribute the apparent success ofthe shifted spectra of ScO − computed with PBE for the firstthree peaks to the observation that the three highest occupiedorbitals ( b , a , b ) are delocalized and primarily due to O2 p x , 2 p y , and 2 p z states, respectively, with negligible ( < d character. The orbitals with the largest ( ∼ d con-tent for this molecule are the 1 a (with d xy character) and 5 b (with d xz character) orbitals, for which the agreement with ex-perimental data is not as satisfactory. Both the shifted PBEand shifted PBE0 spectra have the same mean absolute error(MAE) of 0.22 eV.Among the variants of the GW approximation, the pre-dictions from the one-shot GW with a PBE starting point( G W @PBE) are particularly poor. All quasiparticle ener-gies are underestimated significantly (by more than ∼ ↑ b ↑ a ↑ ↑ a ↑ ↑ ↑ ↑↑ ↑ ↑↑ FIG. 1. Experimental photoelectron spectrum of ScO − (Ref. 47) along with spectra computed with shifted PBE, shiftedPBE0, G W @PBE, G n W @PBE, G n W n @PBE, G W @PBE0 and G W @PBE α = . . Contour plots for some of the occupied molec-ular orbitals are shown on the right, with matching color codes dis-played in the spectra. this level of theory. Predictions from G W calculations withhybrid functional starting points are much better, as shown inFig. 1. In fact, both G W @PBE0 and G W @PBE α = . lev-els of theory have the lowest MAE of 0.21 eV averaged overthe four experimentally measured peaks. It is important, how-ever, to note that not all peaks are predicted equally well atthese two levels of theory: G W @PBE0 predictions for thefirst three peaks are quite accurate (within ∼ .
15 eV), whilethe fourth peak is underestimated by 0.45 eV, very similar tothe case of the shifted PBE0. G W @PBE α = . predictions,on the other hand, are not as good for the first three peaks(off by ∼ . d char-acter of the relevant orbitals: Typically, orbitals with larger3 d content require a larger amount of exact exchange in the G W starting point, as also observed for the case of copperoxide molecular anions. Another level of theory that leads togood agreement with experimental data is G n W @PBE, whichhas been argued to be a practical GW scheme for the elec-tronic structure of transition metal oxide molecular systemsas a good compromise between computational efficiency andaccuracy. The MAE for G n W @PBE is 0.24 eV. Applyingeigenvalue self-consistency also in W ( G n W n ), on the otherhand, deteriorates the agreement with experimental data, andthe MAE for this level of theory is 0.76 eV. B. TiO − Figure 2 shows the experimental photoelectron spectrumof TiO − and the spectra calculated with shifted DFT andvarious GW flavors. The most recent experimental spec-trum obtained with a 193 nm (6.42 eV) laser consists offour main peaks at 1.60, 3.90, 4.72, and 5.38 eV. Withboth PBE and PBE0 functionals, we find the ground stateof TiO − to be a doublet ( A ) with a valence configuration ( a ) ( b ) ( a ) ( b ) ( a ) , and the ground state of theneutral molecule is a singlet with the electron removed fromthe very delocalized 10 a orbital with large Ti 4 s content andsmall Ti p z and 3 d x − z admixtures. As in ScO − , the 6 b stateis entirely due to O 2 p z states, but it is considerably more lo-calized in TiO − . The higher binding energy (BE) orbitals 9 a and 3 b are also more localized in TiO − , and they have larger3 d content ( ∼ − .The vertical IPs calculated with PBE and PBE0 function-als are 1.54 and 1.66 eV, bracketing the experimental value of1.60 eV. In the following comparisons, we compare the exper-imental values for the higher BE peaks with the average of thecalculated spin-up and spin-down eigenvalues of a given dou-bly occupied orbital. While the vertical IPs calculated withboth PBE and PBE0 are in very good agreement with experi-ment, shifted PBE clearly fails to provide satisfactory predic-tions when higher BE peaks are taken into account, as shownin Fig. 2. Shifted PBE0, on the other hand, performs quitewell for the first three peaks and underestimates the energy ofthe fourth. This observation is another example of the well-known failure of semi-local functionals, such as PBE, as theyunderbind localized orbitals due to SIE, while the addition ofa fraction of exact exchange globally, as in PBE0, or in range-separated form helps to mitigate this error. The overall MAEof shifted PBE and PBE0 are found to be 0.77 and 0.21 eV, re-spectively. These observations are in good agreement with the ↑ b ↑ a ↑ ↑ a ↑↑↑↑↑ FIG. 2. Experimental photoelectron spectrum of TiO − (Ref. 51) along with spectra computed with shifted PBE, shiftedPBE0, G W @PBE, G n W @PBE, G n W n @PBE, G W @PBE0 and G W @PBE α = . . Contour plots for some of the occupied molecu-lar orbitals, in both majority ( ↑ ) and minority ( ↓ ) spin channels, areshown on the right, with matching color codes displayed in the spec-tra. For doubly occupied orbitals, the calculated energy displayedwith the higher (lower) height refers to the orbital in the majority(minority) spin channel. results of Marom et al. who reported similar findings.Among the variants of the GW approximation, G W @PBEhas the worst performance with an MAE of 1.68 eV. Similar toScO − , the best overall agreement with experiment is achievedwith the G W @PBE0 level of theory (with an MAE of 0.16eV), while adding more exact exchange worsens the agree-ment with experiment, as the MAE of G W @PBE α = . is0.51 eV, significantly larger than that of ScO − . The effect ofeigenvalue self-consistency in predicting quasiparticle ener-gies of TiO − is similar to that of ScO − . G n W provides fairlygood agreement with experiment (MAE of 0.25 eV), while it-erating eigenvalues in W worsens the agreement by overbind-ing all orbitals other than the highest occupied molecular or-bital (HOMO), leading to an MAE of 0.77 eV. C. VO − Previous studies on the electronic structure of VO − have re-vealed many possible low energy configurations as candidatesfor its ground state. Based on an earlier studyon the electronic structure of the neutral molecule, Wu andWang interpreted their experimental photoelectron spectrumof VO − by assuming that the ground state of the anion is a sin-glet of A symmetry. Later computational studies showedthat the A state is higher in energy than the triplet states B and A by 0 . − . B or A , typically by a very small energydifference. For example, the most recent DFT calculationsof Kim et al. and RCCSD(T) calculations of Hendrickx andTran have found B as the ground state by small differencesof 0.02 to 0.07 eV relative to A .In our DFT calculations with PBE functional, we find theground state of VO − to be the B triplet with the valenceconfiguration ( b ) ( a ) ( b ) ( a ) ( b ) . The groundstate of the neutral molecule is found to be a doublet withthe electron removed from the 4 b orbital, which has a large( ∼ d xz character and is significantly more localized thanthe 10 a orbital. As a result, with the addition of a frac-tion of exact exchange, which reduces SIE, the 4 b orbitalmoves down in energy much more than the 10 a orbital. Ac-cordingly in our calculations with the PBE0 functional, whilethe ground state of VO − is still found to be the B triplet,the HOMO is the 10 a orbital. The downward shifts in theKohn-Sham eigenvalues of the 4 b and 10 a orbitals in goingfrom PBE to PBE0 are 1.83 and 0.90 eV, respectively, whichleads to the observed reordering of the orbitals. We note thatour PBE calculations find the A singlet state to be 0.4 eVhigher in energy, in agreement with previous DFT and quan-tum chemistry calculations. Below, we only discuss the spec-tra obtained for the triplet configuration. As we show in thesupplementary material, none of the theoretical spectra ob-tained for the singlet configuration resembles the experimen-tal spectrum, which provides more evidence that the moleculesampled in the experiments is a triplet.Figure 3 shows the experimental photoelectron spectrum ofVO − obtained with a 193 nm laser and our spectra calcu-lated at various levels of theory. The experimental spectrumconsists of four main peaks at 2.03, 2.75, 4.10, and 4.85 eV. Inour comparisons with the third and fourth experimental peaks,we take the average of our calculated spin-up and spin-downeigenvalues of doubly occupied orbitals. The vertical IPs cal-culated with PBE and PBE0 are 1.84 and 2.04 eV, the latterbeing in nearly perfect agreement with experiment. The othermost obvious difference between shifted PBE and PBE0 spec-tra is the splitting of the first two peaks, which are calculatedas 0.26 and 0.67 eV, respectively, compared to the experimen-tal value of 0.72 eV. In fact, the whole shifted PBE0 spec-trum has excellent agreement with experimental data with anMAE of 0.06 eV. The agreement with shifted PBE, on theother hand, is not as good with a MAE of 0.34 eV. ↑ b ↑ a ↑ ↑ ↑↑↑↑ ↑↑ FIG. 3. Experimental photoelectron spectrum of VO − (Ref. 48) along with spectra computed with shifted PBE, shiftedPBE0, G W @PBE, G n W @PBE, G n W n @PBE, G W @PBE0 and G W @PBE α = . . Contour plots for some of the occupied molecu-lar orbitals, in both majority ( ↑ ) and minority ( ↓ ) spin channels, areshown on the right, with matching color codes displayed in the spec-tra. For doubly occupied orbitals, the calculated energy displayedwith the higher (lower) height refers to the orbital in the majority(minority) spin channel. The predictions from the one-shot GW approximation withPBE, PBE0, and PBE α = . starting points follow the simi-lar trends discussed earlier. G W @PBE predictions are verypoor, while G W @PBE0 provides excellent agreement withexperiment with an MAE of 0.09 eV, very similar to shiftedPBE0. More exact exchange in the starting point worsens theagreement with experiment by overbinding all orbitals signifi-cantly, leading to an MAE of 0.58 eV at the G W @PBE α = . level of theory. Eigenvalue self-consistency in G with PBEstarting point interestingly does not provide as good agree-ment with experiment compared to other molecules. Whilethe G n W @PBE prediction for the IP is very good (1.98 eV),the second peak predicted at 2.00 eV is underestimated by0.75 eV. Here, we should also note that the HOMO eigen-value calculated with G n W @PBE does not correspond to the4 b orbital (which is the HOMO at PBE and G W @PBE lev-els), but to the 10 a orbital. In particular, even though G n W predicts the quasiparticle energies of the 10 a orbital to belower than those of the 4 b orbital calculated with aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets, the eigen-value difference gets smaller as the basis set size increases,such that when the eigenvalues are extrapolated, the 10 a or-bital becomes the HOMO. This difference in the convergencebehavior of 10 a and 4 b orbitals is related to the fact thatthe 4 b orbital has a much higher 3 d content, and increasingthe basis set size for such localized orbitals typically leadsto a much larger lowering of the corresponding quasiparticleenergies. At the G n W n @PBE level, the predicted splittingbetween the 10 a (now clearly the HOMO) and 4 b (HOMO-1) orbitals significantly improves, however, the peaks withhigher BE are overbound at this level of theory, leading toa MAE of 0.64 eV. D. CrO − Figure 4 shows the experimental photoelectron spectrum ofCrO − obtained with a 193 nm laser along with the com-puted spectra. There are three main peaks in the experi-mental spectrum at energies 2.43, 3.41, and 4.25 eV. Withboth PBE and PBE0 functionals, we find the ground stateof CrO − to be the B quartet with the valence configuration ( b ) ( a ) ( b ) ( a ) . The ground state of the neutralmolecule is the B triplet with the electron removed from the11 a orbital. Unlike the HOMO of VO − , which is a local-ized orbital with a large 3 d xz character, the HOMO (11 a ) ofCrO − is a delocalized orbital of primarily Ti 4 s character withsome Ti p z and 3 d x − z admixture, similar to the 10 a orbitalsin TiO − and VO − . The frontier orbitals of large ( ∼ d character in CrO − are HOMO-1 and HOMO-2, whose DFTeigenvalues undergo a large downward shift ( ∼ ∼ a orbitalremains as the HOMO with PBE0 also, and there is no or-bital reordering as observed in VO − . We also note that in thecomputed photoelectron spectra, since the Kohn-Sham eigen-values or the GW quasiparticle energies of the 10 a and 4 b states are quite close to each other, we assign their average tothe measured peak at 3.41 eV (a similar assignment was donein Ref. 50), and the measured peak at 4.25 is compared withthe average of the spin-up and spin-down 6 b orbital energies.The vertical IPs calculated with PBE and PBE0 function-als are 2.30 and 2.48 eV, respectively, which are in reason-ably good agreement with experiment. However, as shownin Fig. 4, the shifted PBE does not have a very good agree-ment with experiment overall (MAE of 0.36 eV), since thesecond peak is underestimated significantly. This is not sur-prising, as this peak is due to 10 a and 4 b orbitals with large3 d character, which both suffer from SIE. As expected, theaddition of a fraction of exact exchange mitigates this prob-lem, and the shifted PBE0 spectrum has a much better agree-ment with experiment with a very low MAE of 0.06 eV. The GW predictions are overall similar to the general trends dis-cussed earlier. One-shot GW with PBE0 starting point has thebest overall agreement with experiment with an MAE of 0.09eV, and each predicted peak matches almost perfectly with the ↑ b ↑↑↑ ↑ FIG. 4. Experimental photoelectron spectrum of CrO − (Ref. 50)along with spectra computed with PBE, PBE0, G W @PBE, G n W @PBE, G n W n @PBE, G W @PBE0 and G W @PBE α = . .Contour plots for some of the occupied molecular orbitals, in bothmajority ( ↑ ) and minority ( ↓ ) spin channels, are shown on the right,with matching color codes displayed in the spectra. For doubly oc-cupied 6 b orbital, the calculated energy displayed with the higher(lower) height refers to the orbital in the majority (minority) spinchannel. shifted PBE eigenvalues. Different from VO − (but in line withother molecules), the G n W @PBE predictions are quite goodwith an MAE of 0.18 eV. As before, adding self-consistencyin W eigenvalues or one-shot GW with a larger exact exchangefraction worsen the agreement with experiment by overbind-ing the orbitals, and G W @PBE predicts orbitals which aretoo underbound, resulting in a very poor MAE of 1.66 eV. E. MnO − In their combined experimental and theoretical study of thephotoelectron spectra of MnO − x clusters, Gutsev et al. founda strong temperature dependence of the spectrum for MnO − .The experimental data displayed in Fig. 5 shows the spectrumobtained with a 266 nm laser under cold conditions. The smallpeak near 1.89 eV was attributed to a higher energy isomerand will not be considered in the following discussion. Thevertical IP of the MnO − was found to be 2.26 eV based onthe analysis of the vibrational progression between 2.0 to 2.7eV, followed by a broad peak at 3.09 eV. This is followed by acongested spectrum with many overlapping features between3.6 and 4.3 eV, from which we have identified peaks at 3.67,3.77, 3.85, 3.96, and 4.11 eV, giving a total of 7 peaks forcomparison with our calculated spectra.With both PBE and PBE0 functionals, we find the groundstate of MnO − to be the B quintet with a valence configu-ration ( b ) ( b ) ( a ) ( b ) ( a ) ( a ) . We find theground state of the neutral molecule to be the B quartet withthe electron removed from the 2 a orbital, which is a localizedmolecular orbital derived primarily from the Mn 3 d xy atomicorbital hybridized with O 2 p x orbitals. Similar to the case ofCrO − , 11 a is rather delocalized with significant 4 s character,while 4 b and 10 a are both localized and derived primarilyfrom Mn 3 d xz and 3 d x − z , respectively. In the following dis-cussion, we compare the energies of five majority-spin orbitalenergies and two minority-spin orbital (6 b ↓ and 3 b ↓ ) ener-gies with the seven peaks from the experimental spectrum.The vertical IPs calculated with PBE and PBE0 functionalsare 2.13 and 2.54 eV, respectively, compared to the experi-mental value of 2.26 eV. The deviation of the IP calculatedwith PBE0 from the experimental value is relatively large incomparison to other molecules considered here, however, therest of the shifted PBE0 spectrum is in very good agreementwith experiment leading to a MAE of 0.19 eV. While the IPcalculated with PBE is very close to the experimental value,overall the shifted PBE spectrum with a MAE of 0.30 eV doesnot have as good agreement with experiment as PBE0. Themost obvious difference between the two spectra is the split-ting between the first two peaks, which are found to be 1.19and 0.61 eV with PBE and PBE0, respectively, compared tothe experimental value of 0.83 eV. This can be understood interms of the difference in the spatial extents of the relevantorbitals, 2 a and 11 a , and the downward shift in energy theyexperience with the addition of a fraction of exact exchange,as discussed for several cases earlier. In the case of MnO − , theDFT eigenvalue of the more localized 2 a orbital decreases by1.87 eV, while that of the less localized 11 a orbital decreasesby 1.29 eV in going from PBE to PBE0 (keep in mind that thePBE and PBE0 spectra plotted in Fig. 5 are shifted spectra).The predictions from the variants of the GW approxima-tion are consistent with the overall trends discussed so far.Best agreement is achieved with G W @PBE0 (MAE of 0.12eV), followed by G n W @PBE with a MAE of 0.37 eV. Whilethe first two peaks are predicted well with G n W , states withhigher BEs are overbound by 0 . − . W overbinds these orbitals more,and MAE increases to 0.97 eV. PBE starting point for one-shot GW again leads to the worst agreement with experiment(MAE of 1.23 eV), and G W @PBE α = . overbinds all statesby 0 . − . ↑ a ↑ a ↑ ↑ a ↑ b ↑ ↑ FIG. 5. Experimental photoelectron spectrum of MnO − (Ref. 49)along with spectra computed with PBE, PBE0, G W @PBE, G n W @PBE, G n W n @PBE, G W @PBE0 and G W @PBE α = . .Contour plots for some of the occupied molecular orbitals, in bothmajority ( ↑ ) and minority ( ↓ ) spin channels, are shown on the right,with matching color codes displayed in the spectra. For doubly oc-cupied 6 b orbital, the calculated energy displayed with the higher(lower) height refers to the orbital in the majority (minority) spinchannel. IV. DISCUSSION
The overall trends for how the predictions from various lev-els of DFT and GW theories compare with experimental datafor all TMO anions are summarized in Fig. 6, which showsthe differences between the computed and experimental ver-tical IPs (left panel) as well as the MAEs (right panel). Wenotice that the vertical IPs computed with the PBE functionalare quite good, underestimated slightly (0 . − . − with an MAE of 0.77 eV, which high-lights the importance of comparing not just the vertical IP butall the higher BE peaks in available PES in assessing the per-formance of a particular level of theory. For shifted PBE0,on the other hand, not only are the predictions for vertical IPsin reasonably good agreement with experiment (typically lessthan 0.15 eV, but with an outlier for the case of MnO − where -1-0.5 0 0.5 1 1.5 2 S c O - T i O - V O - C r O - M n O - E ex pb i nd - E ca l b i nd ( e V ) PBEPBE0G W @PBEG n W @PBEG n W n @PBEG W @PBE0G W @PBE α =0.5 S c O - T i O - V O - C r O - M n O - M A E ( e V ) FIG. 6. Difference between experimental and our calculated vertical IPs ( E expbind − E calbind ) (left), and mean absolute error (MAE) of calculatedelectron binding energies (right) of TMO − molecules at different levels of theory. the deviation from experiment is 0.28 eV), but also the overallMAEs (less than ∼ . G W @PBE α level of theory, thestarting point makes a striking difference. As highlighted forall TMO anions in the previous section, a PBE starting pointleads to very poor agreement with experiment for both the ver-tical IPs (average underestimate compared to experiment forall cases considered being 1.36 eV) and the MAEs, which areslightly higher. PBE0 starting point, on the other hand, leadsto the best agreement overall with experimental data amongall levels of theory considered in this study, not only for verti-cal IPs (within ∼ . . − . G W @PBE α = . ,all orbitals are overbound by ∼ . − where MAE remains near 0.2 eV, close to the valueobtained with α = .
25. This is consistent with the results ob-tained for ScO − in Ref. 43, where any value of α in the DFTstarting point in the range 0 . ≤ α ≤ d character in therelevant orbitals. We also note that our finding of nearly ex-cellent agreement of G W @PBE α at the “sweet spot” valueof α = .
25 should not be taken as a general trend for the G W @PBE0 level of theory for TMO molecular systems, asthis work has focused particularly on early TM dioxide sys-tems. Typically, the amount of α needed for good agree-ment with experimental data in the DFT starting point of one-shot GW calculations for TMO molecules is dependent on theamount of the 3 d character of the orbitals as well as thecontent of TM in the molecule. For example, Shi et al. showed in their study of small copper oxide molecular anionsthat while α = . O − andCuO − , values near α = .
25 were optimal for molecules withless Cu content, such as CuO − and CuO − .For the case of ev GW , updating eigenvalues only in G ( G n W ) with a PBE starting point leads to very good agree-ment with experiment for the vertical IPs, typically within0.05 eV, except for the two slight outliers of ScO − and MnO − where the deviation is ∼ .
25 eV. When higher BE peaks areconsidered, the agreement with experiment slightly worsens,but the MAE value of 0.3 eV averaged over the five TMO anions is still reasonable. Adding self-consistency in W , onthe other hand, overbinds all orbitals, as the poles of the self-energy become more negative (for occupied orbitals) with theincrease in the neutral excitation energies. For vertical IPs ofTiO − , VO − , and CrO − , G n W n @PBE predictions are in goodagreement with experiment (overestimated by 0.11, 0.23, and0.40 eV, respectively), but the predicted IPs for ScO − andMnO − are not, as they are overestimated by ∼ . G n W n @PBE is clearly not a predic-tive level of theory, as the MAE averaged over the five TMO anions is ∼ GW theories for the TMO anionsconsidered, we now focus on the (averaged) eigenvalue shiftsfor each molecule, as shown in Fig. 7. For the GW levels oftheory, the bars show for each TMO anion the magnitudes ofthe shifts of the GW eigenvalues from the corresponding DFTstarting point eigenvalues averaged over the measured peaks.The standard deviations (SDs) of the shifts are shown witherror bars for each level of GW theory. For PBE and PBE0,the bars show the magnitudes of the shifts needed to alignthe HOMO eigenvalue with negative of the computed verticalIP of the anion discussed earlier. One immediate observationfrom Fig. 7 is that the eigenvalue shifts for the G W @PBE0(blue bars) are very close to the PBE0 shifts (cyan bars) forall molecules considered. This, combined with the observa-tion that the SDs of the shifts for G W @PBE0 are quite small(less than 0.1 eV), shows that the G W @PBE0 spectra can beobtained by almost a uniform shift of the PBE0 spectra by theamount needed to align HOMO with negative the vertical IPof the anion, irrespective of particular nature (localized, delo-calized, large or small 3 d character) of the orbitals, resulting S c O - T i O - V O - C r O - M n O - E i g e n va l u e s h i ft ( e V ) PBEG W @PBEG n W @PBEPBE0G W @PBE010a FIG. 7. Calculated eigenvalue shifts for TMO − molecules at eachlevel of theory. For PBE and PBE0, the bars show the magnitudesof shifts needed to align the HOMO eigenvalues with negative of thecomputed vertical IPs. For GW levels of theory, the bars show themagnitudes of the shifts of the GW eigenvalues from the eigenvaluesof the corresponding DFT starting points averaged over the measuredpeaks. The error bars are the standard deviations of the shifts foreach level of GW theory. The black dots show the shifts for thedelocalized 10 a and 11 a orbitals discussed in the text. This orbitalis unoccupied in ScO − . also in similar MAEs with respect to experimental data. Thevery good agreement of the G W @PBE0 level of theory withexperimental PES can then be attributed to the correct ener-getic positioning of the orbitals obtained with PBE0 and thePBE0 wave functions being good approximations to the truequasiparticle wave functions, so that a simple first-order G W perturbative correction is enough to lead to accurate quasipar-ticle energies.With a PBE starting point, on the other hand, the av-erage shifts for G W @PBE (yellow bars in Fig. 7) aremuch smaller than the PBE eigenvalue shifts required to alignHOMO with negative the vertical IP (red bars in Fig. 7).That is, unlike the very good agreement of shifted PBE0 and G W @PBE0 with each other (as well as with experiment),application of first-order G W perturbative correction to thePBE eigenvalues using PBE wave functions fails to providethe required amount of downward shift to achieve a reason-ably good level of agreement with experimental data. Toquantify this further, we compared the PBE and PBE0 eigen-values for all TMO anions and considered the amounts ofshifts that would be needed to achieve perfect agreement withexperimental data. Averaged over all experimental peaks andall TMO anions, the required shift is 3 . ± .
27 eV and1 . ± .
13 eV for PBE and PBE0, respectively, showing (i)how underbound PBE orbitals are for perturbative correctionsto be effective, and that (ii) the required shifts for PBE have amuch larger variation than those for PBE0 due to the signifi-cant difference in the SDs.These last two observations are also consistent with the trends of eigenvalue shifts for G n W @PBE shown with greenbars 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 partic-ular orbital, due to the large SD observed for most TMO molecular anions, with the exception of ScO − , where the oc-cupied frontier orbitals are composed of primarily O 2 p statesand have similar spatial extents. The large SD observed in G n W @PBE eigenvalue shifts is primarily due to a particularorbital with a character (10 a in TiO − and VO − , and 11 a in CrO − and MnO − ), which is either the HOMO or HOMO-1in the majority spin channel depending on the level of theory(see Fig. 3 for VO − ) or the molecule (in MnO − , HOMO isthe 2 a orbital). The shifts for this particular orbital also dis-played in Fig. 7 are significantly less than the averaged shift, e.g. the G n W @PBE eigenvalue for the 10 a state of TiO − is2.2 eV lower than the PBE eigenvalue, which is much smallerthan the average downward shift of 3.27 eV. As a result, theSD is particularly high at the G n W @PBE level of theory, asdifferent orbitals are shifted by varying amounts to achievea reasonably good level of agreement with experiment. Asmentioned earlier, this particular a orbital has a large TM 4 s character with some TM p z and 3 d x − z admixtures. It is im-portant to note that while the 3 d character is not very large,it is still appreciable, ranging from 15 to 30% in going fromTiO − to MnO − . The primary reason why this orbital does notundergo a large downward shift is its very delocalized naturein spite of having an appreciable 3 d content, consistent withthe observation that G n W @PBE does a particularly good jobfor the vertical IPs, and the agreement deteriorates slightly forthe more localized higher BE peaks. V. SUMMARY
In summary, we provided a detailed comparison of theDFT- and GW -computed eigenvalue spectra to experimen-tal photoelectron spectra of five 3 d -transition metal dioxidemolecular anions, ScO − , TiO − , VO − , CrO − , and MnO − . Wefocused our study on the comparison between semi-local andhybrid functional predictions (within DFT and as a startingpoint for G W calculations) as well as the effects of eigen-value self-consistency with a semi-local DFT starting point.Overall, both the shifted PBE0 and G W @PBE0 appear tostand out as the best levels of theory, as they are able to re-produce the experimental BEs with 0 . − . GW calcula-tion with a semi-local starting point like PBE is clearly a poorchoice leading to the worst agreement with experiment (aver-age error ∼ GW levels of theory. However,too much Fock exchange (e.g. α = .
5) also worsens theagreement with experiment by overbinding orbitals. We foundthat G − only eigenvalue self-consistent GW with W computedfrom PBE (G n W @PBE), while not as accurate as shiftedPBE0 or G W @PBE0, still provides a reasonably good de-scription of the quasiparticle energies for occupied frontierorbitals in these transition metal oxide molecular systems.Updating the eigenvalues in W ( G n W n @PBE), on the otherhand, leads to poor agreement with experiment overall, eventhough it may be reasonably accurate for vertical IPs of someof the molecular anions considered. We caution the reader thatour finding of very good agreement of G W @PBE0 predic-tions with experimental data for the particular case of theseearly 3 d − transition metal dioxide molecular anions shouldbe viewed in the context of our previous studies on othertransition metal oxide molecular systems, which showedthat higher values of Fock exchange might be needed in theDFT starting point of GW calculations to achieve good agree-ment with experiment in systems with more localized orbitalsand higher transition metal content. When viewed from thisperspective, the results presented in this paper still lead usto advocate using G n W @PBE as a practical GW schemethat presents a good balance between accuracy and efficiencyin predicting quasiparticle energies of transition metal oxidemolecular systems without the need for a system-dependentparameter in the starting DFT description. Future studies onthe more challenging late 3 d − transition metal dioxide molec-ular anions, such as FeO − , NiO − and CuO − , and other transi-tion metal oxide molecules with varying amounts of transitionmetal content are needed to further assess and understand thepredictive capabilities of GW methods for electronic excita-tions in molecular systems with moderate electron correlation. SUPPLEMENTARY MATERIAL
See supplementary material for more details, results, anddiscussion.
ACKNOWLEDGMENTS
This work was supported by the U.S. Department of EnergyGrant No. DE-SC0017824. This research used resources ofthe National Energy Research Scientific Computing Center,a DOE Office of Science User Facility supported by the Of-fice of Science of the U.S. Department of Energy under Con-tract No. DE-AC02-05CH11231. We would also like to thankYoung-Moo Byun for useful discussions in the earlier stagesof this work.
REFERENCES L. Reining, WIREs Comput. Mol. Sci. , e1344 (2018). P. Rinke, A. Qteish, J. Neugebauer, C. Freysoldt, and M. Scheffler, New J.Phys. , 126 (2005). F. Bruneval, N. Vast, and L. Reining, Phys. Rev. B , 045102 (2006). M. L. Tiago and J. R. Chelikowsky, Phys. Rev. B , 205334 (2006). M. Shishkin and G. Kresse, Phys. Rev. B , 235102 (2007). F. Fuchs, J. Furthmüller, F. Bechstedt, M. Shishkin, and G. Kresse, Phys.Rev. B , 115109 (2007). C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen, Phys. Rev. B ,085103 (2010). Y. Ma, M. Rohlfing, and C. Molteni, J. Chem. Theory Comput. , 257–265(2010). X. Blase and C. Attaccalite, Appl. Phys. Lett. , 171909 (2011). X. Blase, C. Attaccalite, and V. Olevano, Phys. Rev. B , 115103 (2011). C. Faber, C. Attaccalite, V. Olevano, E. Runge, and X. Blase, Phys. Rev. B , 115123 (2011). X. Qian, P. Umari, and N. Marzari, Phys. Rev. B , 075103 (2011). S.-H. Ke, Phys. Rev. B , 205415 (2011). X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko, A. Sanfilippo,K. Reuter, and M. Scheffler, New J. Phys. , 053020 (2012). N. Marom, F. Caruso, X. Ren, O. T. Hofmann, T. Körzdörfer, J. R. Che-likowsky, A. Rubio, M. Scheffler, and P. Rinke, Phys. Rev. B , 245127(2012). B. Baumeier, D. Andrienko, Y. Ma, and M. Rohlfing, J. Chem. TheoryComput. , 997–1002 (2012). M. J. van Setten, F. Weigend, and F. Evers, J. Chem. Theory Comput. ,232–246 (2013). F. Bruneval and M. A. L. Marques, J. Chem. Theory Comput. , 324–329(2013). T. A. Pham, H.-V. Nguyen, D. Rocca, and G. Galli, Phys. Rev. B ,155148 (2013). F. Caruso, P. Rinke, X. Ren, A. Rubio, and M. Scheffler, Phys. Rev. B ,075105 (2013). V. Atalla, M. Yoon, F. Caruso, P. Rinke, and M. Scheffler, Phys. Rev. B ,165122 (2013). J. Klimeš, M. Kaltak, and G. Kresse, Phys. Rev. B , 075125 (2014). C. Faber, P. Boulanger, C. Attaccalite, I. Duchemin, and X. Blase, Philos.Trans. R. Soc. A , 20130271 (2014). P. Koval, D. Foerster, and D. Sánchez-Portal, Phys. Rev. B , 155417(2014). P. Boulanger, D. Jacquemin, I. Duchemin, and X. Blase, J. Chem. TheoryComput. , 1212–1218 (2014). S. Körbel, P. Boulanger, I. Duchemin, X. Blase, M. A. L. Marques, andS. Botti, J. Chem. Theory Comput. , 3934–3943 (2014). L.-W. Wang, Phys. Rev. B , 125135 (2015). D. Hirose, Y. Noguchi, and O. Sugino, Phys. Rev. B , 205111 (2015). F. Bruneval, S. M. Hamed, and J. B. Neaton, J. Chem. Phys. , 244101(2015). D. Jacquemin, I. Duchemin, and X. Blase, J. Chem. Theory Comput. ,3290–3304 (2015). 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. , 5665–5687(2015). J. Wilhelm, M. D. Ben, and J. Hutter, J. Chem. Theory Comput. , 3623–3635 (2016). F. Kaplan, M. E. Harding, C. Seiler, F. Weigend, F. Evers, and M. J. vanSetten, J. Chem. Theory Comput. , 2528–2541 (2016). F. Caruso, M. Dauth, M. J. van Setten, and P. Rinke, J. Chem. TheoryComput. , 5076–5087 (2016). 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. ,615–626 (2016). L. Hung, F. H. da Jornada, J. Souto-Casares, J. R. Chelikowsky, S. G. Louie,and S. Ö˘güt, Phys. Rev. B , 085125 (2016). L. Hung, F. Bruneval, K. Baishya, and S. Ö˘güt, J. Chem. Theory Comput. , 2135–2146 (2017). L. Hung, F. Bruneval, K. Baishya, and S. Ö˘güt, J. Chem. Theory Comput. , 5820–5821 (2017). E. Maggio, P. Liu, M. J. van Setten, and G. Kresse, J. Chem. Theory Com-put. , 635–648 (2017). M. Govoni and G. Galli, J. Chem. Theory Comput. , 1895–1909 (2018). M. Grumet, P. Liu, M. Kaltak, J. Klimeš, G. Kresse, M. Kaltak, andG. Kresse, Phys. Rev. B , 155143 (2018). B. Shi, S. Weissman, F. Bruneval, L. Kronik, and S. Ö˘güt, J. Chem. Phys. , 064306 (2018). Y.-M. Byun and S. Ö˘güt, J. Chem. Phys. , 134305 (2019). W. Gao and J. R. Chelikowsky, J. Chem. Theory Comput. , 5299–5307(2019). D. Golze, M. Dvorak, and P. Rinke, Front. Chem. , 377 (2019). H. Wu and L.-S. Wang, J. Chem. Phys. , 8221 (1997). H. Wu and L.-S. Wang, J. Phys. Chem. A , 9129–9135 (1998). H. Wu and L.-S. Wang, J. Chem. Phys. , 5310 (1998). G. L. Gutsev, B. K. Rao, P. Jena, X. Li, and L.-S. Wang, J. Chem. Phys. , 1473 (2000). G. L. Gutsev, P. Jena, H.-J. Zhai, and L.-S. Wang, J. Chem. Phys. , 7935(2001). H.-J. Zhai and L.-S. Wang, J. Am. Chem. Soc. , 3022–3026 (2007). J. B. Kim, M. L. Weichman, and D. M. Neumark, J. Chem. Phys. ,034307 (2014). L. B. Knight, R. Babb, M. Ray, T. J. Banisaukas, L. Russon, R. S. Dailey,and E. R. Davidson, J. Chem. Phys. , 10237 (1996). G. V. Chertihin, W. D. Bare, and L. Andrews, J. Phys. Chem. A , 5090–5096 (1997). M. B. Walsh, R. A. King, and H. F. Schaefer, J. Chem. Phys. , 5224(1999). M. Zhou and L. Andrews, J. Chem. Phys. , 4230 (1999). J. M. Gonzales, R. A. King, and H. F. Schaefer, J. Chem. Phys. , 567(2000). S. F. Vyboishchikov and J. Sauer, J. Phys. Chem. A , 10913–10922(2000). G. L. Gutsev, B. K. Rao, and P. Jena, J. Phys. Chem. A , 11961–11971(2000). J. Dong, Y. Wang, and M. Zhou, Chem. Phys. Lett. , 511–516 (2002). F. Grein, J. Chem. Phys. , 034313 (2007). S. Li and D. A. Dixon, J. Phys. Chem. A , 6646–6666 (2008). E. L. Uzunova, H. Mikosch, and G. S. Nikolov, J. Chem. Phys. , 094307(2008). Y. Gong, M. Zhou, and L. Andrews, Chem. Rev. , 6765–6808 (2009). Y. Liu, Y. Yuan, Z. Wang, K. Deng, C. Xiao, and Q. Li, J. Chem. Phys. , 174308 (2009). E. P. F. Lee, D. K. W. Mok, F.-T. Chau, and J. M. Dyke, J. Comput. Chem. , 337–345 (2009). Z.-W. Qu and H. Zhu, J. Comput. Chem. , 2038–2045 (2010). N. Marom, J. E. Moussa, X. Ren, A. Tkatchenko, and J. R. Chelikowsky,Phys. Rev. B , 245115 (2011). D. J. Taylor and M. J. Paterson, J. Chem. Phys. , 1–10 (2012). J. B. Kim, M. L. Weichman, and D. M. Neumark, Phys. Chem. Chem.Phys. , 20973 (2013). M. F. Hendrickx and V. T. Tran, J. Chem. Theory Comput. , 4037–4044(2014). G. Xu and F. Yu, J. Electron Spectros. Relat. Phen. , 10–22 (2015). M. Valiev, E. Bylaska, N. Govind, K. Kowalski, T. Straatsma, H. V. Dam,D. Wang, J. Nieplocha, E. Apra, T. Windus, and W. de Jong, Comput. Phys.Commun. , 1477 (2010). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996). J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. , 9982 (1996). F. Bruneval, T. Rangel, S. M. Hamed, M. Shao, C. Yang, and J. B. Neaton,Comput. Phys. Commun. , 149 (2016). C. Hättig, W. Klopper, A. Köhn, and D. P. Tew, Chem. Rev. , 4–74(2012). J. G. Hill and J. A. Platts, J. Chem. Phys.128