Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correction
Puskar Bhattarai, Biswajit Santra, Kamal Wagle, Yoh Yamamoto, Rajendra R. Zope, Adrienn Ruzsinszky, Koblar Alan Jackson, John P. Perdew
EExploring and enhancing the accuracy of interior-scaled Perdew-Zungerself-interaction correction
Puskar Bhattarai, a) Biswajit Santra, Kamal Wagle, Yoh Yamamoto, Rajendra R. Zope, AdriennRuzsinszky, Koblar A. Jackson, and John P. Perdew
1, 4 Department of Physics,Temple University, Philadelphia,PA-19122 Department of Physics, University of Texas at El Paso, El Paso, Texas 79968,USA Department of Physics and Science of Advanced Materials, Central Michigan University, Mount Pleasant,Michigan 48859, USA Department of Chemistry,Temple University, Philadelphia,PA-19122 (Dated: 11 January 2021)
The Perdew-Zunger self-interaction correction (PZ-SIC) improves the performance of density functional ap-proximations (DFAs) for the properties that involve significant self-interaction error (SIE), as in stretchedbond situations, but overcorrects for equilibrium properties where SIE is insignificant. This overcorrection isoften reduced by LSIC, local scaling of the PZ-SIC to the local spin density approximation (LSDA). Here wepropose a new scaling factor to use in an LSIC-like approach that satisfies an additional important constraint:the correct coefficient of atomic number Z in the asymptotic expansion of the exchange-correlation (xc) en-ergy for atoms. LSIC and LSIC+ are scaled by functions of the iso-orbital indicator z σ , which distinguishesone-electron regions from many-electron regions. LSIC+ applied to LSDA works better for many equilibriumproperties than LSDA-LSIC and the Perdew, Burke, and Ernzerhof (PBE) generalized gradient approxima-tion (GGA), and almost as well as the strongly constrained and appropriately normed (SCAN) meta-GGA.LSDA-LSIC and LSDA-LSIC+, however, both fail to predict interaction energies involving weaker bonds,in sharp contrast to their earlier successes. It is found that more than one set of localized SIC orbitalscan yield a nearly degenerate energetic description of the same multiple covalent bond, suggesting that aconsistent chemical interpretation of the localized orbitals requires a new way to choose their Fermi orbitaldescriptors. To make a locally scaled-down SIC to functionals beyond LSDA requires a gauge transformationof the functional’s energy density. The resulting SCAN-sdSIC, evaluated on SCAN-SIC total and localizedorbital densities, leads to an acceptable description of many equilibrium properties including the dissociationenergies of weak bonds. I. Introduction
Kohn-Sham density functional theory(KS-DFT) isa computationally efficient approach to calculate theground state energy and density of many-electron sys-tems. It is widely used to predict various propertiesof atoms, molecules and solids, because of its computa-tional efficiency in comparison with many-electron wave-function calculations. The accuracy of KS-DFT forground state calculations depends upon the approxima-tion for the exchange-correlation energy ( E xc ) of the sys-tem. Jacob’s ladder is often used to classify the approxi-mate density functionals, where higher rungs incorporateadditional ingredients that increase both accuracy andcomputational cost. E xc is often approximated as E xc [ n ↑ , n ↓ , ] = (cid:90) d rn(cid:15) xc ( n ↑ , n ↓ , ∇ n ↑ , ∇ n ↓ , τ ↑ , τ ↓ ) , (1)where n ( r ) = n ↑ ( r ) + n ↓ ( r ) (2) a) Electronic mail: [email protected] is the total density and τ σ ( r ) = occup (cid:88) i |∇ ψ iσ ( r ) | (3)is the positive kinetic energy density for the occupied or-bitals ψ iσ . The first three rungs include the local spindensity approximation(LSDA), generalized gradient ap-proximation (GGA) and meta-GGA. LSDA takes the lo-cal electron density as its sole ingredient. GGAs take thegradient of electron density ∇ n σ as an additional ingre-dient, and meta-GGAs take τ σ as well. Semi-local non-empirical approximate functionals of higher rungs are of-ten constructed by satisfying more exact constraints andnorms than the lower rungs. For instance, the PBE GGA satisfies 11 known exact constraints and the SCAN meta-GGA satisfies 17 known exact constraints. Themost common norm is the uniform electron gas, forwhich almost all the functionals are exact by construc-tion.All semi-local DFAs suffer from self-interaction er-ror (SIE) due to the imperfect cancellation of theself-Hartree energy U by the self-exchange-correlationenergy E xc of a single fully occupied orbital. ThePerdew-Zunger self-interaction correction (PZ-SIC) re-moves one-electron SIE by subtracting it on an orbital-by-orbital basis. Hence, the corrected E xc can be written a r X i v : . [ phy s i c s . c h e m - ph ] J a n hattarai et al. Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correctionas E P Z − SICxc = E approxxc [ n ↑ , n ↓ ] − (cid:88) ασ δ ασ . (4)Here δ ασ is the SIE of an orbital with quantum numbers α and σ , (where σ is the spin, and α is the orbital quantumnumber). The orbital density is n ασ ( r ) = | ψ ασ ( r ) | , and δ ασ = U [ n ασ ] + E approxxc [ n ασ , . (5)Use of delocalized KS orbitals for PZ-SIC leads toa size-extensivity problem. Besides that, they can behighly noded, leading to disastrous results when ap-plied to the semi-local functionals. Fermi-L¨owdin or-bitals (FLOs) correspond to a unitary transformationof KS orbitals and are guaranteed to be localized aroundFermi-orbital descriptor (FOD) positions. FLOs are notnodeless but more weakly noded than the KS orbitals,as the overlapping real FLOs must have nodes to estab-lish orthogonality. Complex FLOs (not yet implemented)can achieve orthogonality without having noded orbitaldensities, but the orbital densities may still be lobed andpose problems for semi-local functionals. The FLOSIC method has been applied to differentdensity functional approximations (DFAs) to successfullypredict various properties of of atoms and molecules .Applying FLOSIC to the approximate semi-local func-tional greatly reduces the error for stretched bonds where SIE is large, but it introduces other errors due tothe presence of nodes and lobed structures in the orbitaldensities, which are not needed by uncorrected DFAs.Besides that, the application of PZ-SIC to the DFAs vi-olates the most common norm of exactness for the uni-form electron gas , and produces a significant error inthe exchange-correlation energy of neutral atoms in thelimit of large atomic number Z . PZ-SIC, although itimproves the barrier heights of chemical reactions thatinclude stretched bonds, often worsens equlibrium prop-erties like atomization energies , electron affinities, ion-ization potential, and bond lengths of molecules .Various scaling functions and methods havebeen introduced to improve the prediction of equilibriumproperties like atomization energy. The half-SIC approx-imation proposed in Ref. 24 does not perform better forboth atomization energy and the energy barriers. Thescaling proposed in Ref. 20 and Ref. 10 violates the cor-rect asymptotic behavior − r of the exchange-correlationpotential and produces − X HO r behavior, where X HO is the scaling factor for the highest occupied orbital.Recently, a new approach called LSIC was proposedas an interior scaling of PZ-SIC. It uses an iso-orbitalindicator( z σ ) as the scaling factor to distinguish one-electron regions from many electron regions. z σ = 0 isconsistent with a region of uniform density and z σ = 1 isconsistent with a one-electron region. z σ is defined by z σ = τ Wσ ( r ) τ σ ( r ) , (6) where, τ Wσ ( r ) is the von Weizs¨acker kinetic energy den-sity τ Wσ ( r ) = |∇ n σ | n σ . (7)LSIC uses the simplest scaling factor z σ to scale downPZ-SIC and is exact for the uniform gas limit as satisfiedby the uncorrected DFAs. LSDA-LSIC, evaluated non-self-consistently on LSDA-SIC FLO densities, improvesmany calculated properties over PZ-SIC, including cova-lent binding energies . However, it fails to recover allof the known asymptotic expansion of E xc for atoms oflarge atomic number Z .The exact large- Z asymptotic expansions of E xc is given as, E xc = − A x Z − A c Zln ( Z ) + B xc Z + C x Z + ...., (8)where LSDA, PBE and almost all other DFAs exactlyreproduce A x and A c . LSDA-sdSIC and PBE-sdSIC do not recover the coefficient B xc , but, as will be shownbelow, SCAN-sdSIC recovers the leading coefficients A x , A c , and B xc and produces a balanced result. The rele-vance of the large- Z expansion to valence-electron prop-erties is discussed in Ref. 32. SCAN-sdSIC improvesthe perfomance of SCAN where SIEs are important andrestores most of the accuracy of SCAN for equilibriumproperties, which are severely degraded by SCAN-SIC. SCAN-sdSIC violates the correct asymptotic behavior − r of the exchange-correlation potential.We present here the scaling function f ( z σ ) applied toPZ-SIC, in an approximation referred to as LSIC+, thatrecovers the leading coefficients A x , A c , and B xc and alsoretains the correct asymptotic behavior of the exchange-correlation potential. LSIC+ is exact for the uniform gaslimit and produces much less error in E xc of the neutralatoms in the limit of large Z . LSIC+ satisfies at leastone more constraint and is expected to work well for bothequilibrium and stretched-bond properties. II. Theory and computational details
The scaled down SIC using LSIC+ is given as E DF A − LSIC + xc = E DF Axc − occ (cid:88) α,σ { U LSIC + [ n ασ , E LSIC + xc [ n ασ , } , (9)where U LSIC + [ n ασ ,
0] = 12 (cid:90) dr f ( z σ ( r )) n α,σ ( r ) (cid:90) dr (cid:48) n α,σ ( r (cid:48) ) | r − r (cid:48) | , (10)2hattarai et al. Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correction E LSIC + xc [ n ασ ,
0] = 12 (cid:90) dr f ( z σ ( r )) n α,σ ( r ) (cid:15) DF Axc ([ n α,σ , , r ) . (11) (cid:15) DF Axc is the exchange-correlation energy density per elec-tron, and f ( z σ ) is the scaling factor such that f (0) = 0, f (1) = 1, i.e., PZ-SIC is scaled in such a way that there isfull correction for any one-electron orbital density and nocorrection for any uniform electron density. For LSIC+,we propose f ( z σ ) as f ( z σ ) = 12 + a ( z σ −
12 ) + b ( z σ −
12 ) (12)where, b = 4(1 − a ). The value of a is chosen to re-cover the correct coefficient B xc of the large- Z asymp-totic expansion of Eq. 8. a = 0 .
5, and b = 2 . Z asymptotic expansion very close to the exactone. Table I shows the B xc values for LSIC, LSIC+ andSCAN-sdSIC. They are based on calculations performedfor four rare gas atoms Neon( Z = 10), Argon ( Z = 18),Krypton ( Z = 36), and Xenon ( Z = 54), and extrapo-lated to Z → ∞ . Here we chose the closed shell atoms toextrapolate and predict the coefficients in order to mini-mize the effect of shell structure as discussed in Ref. 31.Fig. 1 shows the scaling functions of LSIC and LSIC+plotted as a function of z σ . LSIC+ produces more cor-rection than LSIC in the region 0f 0 . < z σ ≤ . . ≤ z σ < . E xc of rare gas atoms, werecarried out using a developmental version of the FLOSICcode based on the NRLMOL code. We performedspin unpolarized calculations for rare gas atoms and theS22 set. All the other calculations are spin polarized.We used the Sadlej basis set that includes long-rangefunctions to capture the extended nature of anion or-bitals for calculating electron affinity. Optimized all-electron Gaussian basis sets were used for all otherproperties. The FODs were optimized for SIC calcu-lations until the maximum component of force is lessthan 5 × − Hartree/Bohr with self-consistent field(SCF) convergence set to 10 − Hartree. Unless other-wise stated, all DFA and DFA-SIC calculations (whereDFA = LSDA, PBE, or SCAN) are self-consistent, whileall DFA-sdSIC, LSDA-LSIC, and LSDA-LSIC+ calcula-tions are performed as a single,non-SCF step using theoptimized DFA-SIC density and Fermi orbital descriptors(FODs).Full self-consistency requires the complicated imple-mentation of the functional derivative of f ( z σ ). There-fore, we used a quasi-SCF version of LSDA-LSIC andLSDA-LSIC+ to obtain the eigenvalue of the highestoccupied molecular orbital (HOMO) for calculating thevertical electron affinity. A more detailed descriptionof quasi-self-consistent LSDA-LSIC and LSDA-LSIC+ ispresented in Refs. 18, 37, and 38, which argue that themain effect of scaling comes from the scaled potentialterm rather than the variation of the scaling factor, whichwould have a minor effect on the results. The quasi-self-consistent version approximates and scales the PZ-SIC FIG. 1:
Scaling functions for LSIC, and LSIC+ asfunctions of z σ . FIG. 2:
Extrapolation to Z − / → B xc of the large- Z asymptotic expansion of Eq. 8 for rare gasatoms Neon( Z = 10), Argon ( Z = 18), Krypton ( Z = 36),and Xenon ( Z = 54), (in Hartree) against Z − / . The solidlines represent the extrapolated curves and the solid dotsrepresent the calculated values at Z = 10 , , , and 54 . Note that A x and A c are exact in LSDA. potential as v interior scaled SICiσ ≈ − f ( z σ ( r )) δ { U [ n iσ ] + E xc [ n iσ , } δn iσ ( r ) . (13) TABLE I:
Coefficient of the third term of the large- Z asymptotic expansion of Eq. 8 using the LSDA-LSIC+,LSDA-LSIC, and SCAN-sdSIC methods. LSIC + LSIC SCAN-sdSIC Exact B xc -0.1806 -0.0828 -0.1755 -0.1868 The LSDA E xc values of rare gas atoms calculated withLSDA-SIC density are subtracted from exact , LSDA-LSIC, and LSDA-LSIC+ E xc values, and the result isdivided by the atomic number Z . The results thus ob-tained are plotted against Z − / as shown in Fig. 2. Theintercept on the y-axis provides the B xc results. The os-cillation of the LSDA exchange curve imitates that of the3hattarai et al. Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correctionexact exchange energy curve. Subtracting the LSDA val-ues minimizes the oscillations due to the shell structure. Fig. 2 shows that the curve for LSDA-LSIC+ lies veryclose to the exact curve whereas LSDA-LSIC deviates toomuch from the exact curve, especially in the higher Z re-gion. We can observe that B xc for LSIC+ is very closeto the exact one, as presented in Table I.Fig. 3 shows the relative percentage error of E xc with LSDA-SIC, LSDA-LSIC, and LSDA-LSIC+ and itdemonstrates that LSDA-LSIC and LSDA-LSIC+ re-spect the norm of the uniform electron gas, since thepercentage error for exchange-correlation energy becomesnearly zero in the limit of large Z . These extrapolatederrors are slightly different from zero since the extrapo-lation is performed taking just four points into account.For Z → ∞ , LSDA-SIC has the highest error of 5.46%,LSDA-LSIC has an error of -0.67% and LSDA-LSIC+has an error of 0.57%. The error of these methods for E xc of rare-gas atoms is shown in table II. TABLE II:
Mean percentage error(MPE) and meanabsolute percentage error(MAPE) of E xc for four rare gasatoms Neon( Z = 10), Argon ( Z = 18), Krypton ( Z = 36),and Xenon ( Z = 54).LSIC+ LSIC LSDA-SICMPE -0.3301 -0.2721 3.6418MAPE 0.3324 0.5492 3.6418 FIG. 3:
Extrapolation as in Ref. 19 of percentage error of E xc for rare gas atoms against Z − / . The solid lines repre-sent the extrapolated curves and the solid dots represent thecalculated values at Z = 10 , , , and 54. III. Results and Discussion
LSIC+ scaling is applied to determine differentground-state properties. We compare several proper-ties of atoms and molecules calculated with LSDA,PBE, and SCAN, with and without SIC, scaled-downSIC(sdSIC) , and the interior-scaling methods LSIC andLSIC+ with LSDA. The interior-scaling methods withPBE and especially SCAN do not produce better re- sults due to gauge inconsistency, which is discussed indetail in Ref. 10. The results for PBE-LSIC and SCAN-LSIC are also presented here for comparison purposesonly. We will discuss the results for the total ener-gies of atoms, atomization energies of molecules, barrierheights of chemical reactions, ionization potentials andelectron affinities of atoms and molecules, equilibriumbond lengths of molecules, and the interaction energiesof a few organic complexes. LSIC and LSIC+ are ap-plied to LSDA, but they produce results close to SCANresults and sometimes even better. Because of the predic-tive power and success of SCAN , the LSDA-LSICand LSDA-LSIC+ results are primarily compared withSCAN, SCAN-SIC , and SCAN-sdSIC results. A. Atoms
The mean error (ME), mean absolute error (MAE) ,and mean absolute percentage error (MAPE) for the totalground state energy of atoms ( Z = 1 −
18) using differentDFAs, their SIC counterparts, and the scaling methodsare presented in Table III. Both LSDA-LSIC and LSDA-LSIC+ significantly reduce the error produced by LSDA-SIC. However, LSIC+ generates slightly more error thanLSIC. SCAN serves as the best functional in estimatingthe total ground state energy of atoms. Fig. 4 shows theperformances of LSDA-LSIC, LSDA-LSIC+, SCAN, andSCAN-SIC in predicting the total ground state energyof atoms. Both LSDA-LSIC and LSDA-LSIC+ producean improved result for Z =1-10 and generates the re-sults closer to the reference values . The error for bothLSDA-LSIC and LSDA-LSIC+ rises rapidly as we movefrom Z =11 to Z =18. FIG. 4:
Error (in Hartree) in the total ground state energyplotted against atomic number Z for atoms from H to Ar. B. Atomization energy
The atomization energy of a molecule is the energyrequired to break it into its constituent atoms. Table IVpresents the MEs, MAEs and MAPEs of the AE6 set et al. Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correction
TABLE III:
ME (in Hartree), MAE (in Hartree), andMAPEs for total ground state energies of atoms for atomicnumber Z = 1 to Z = 18 for various levels of approximation.Method MAE ME MAPELSDA 0.7261 0.7261 -1.0014LSDA-SIC 0.3808 -0.3808 -0.2625LSDA-sdSIC 0.0423 -0.0412 -0.1081LSDA-LSIC 0.0409 0.0147 -0.0753LSDA-LSIC+ 0.0661 0.0618 -0.0694PBE 0.0830 0.0830 -0.0995PBE-SIC 0.1585 0.1585 -0.1091PBE-sdSIC 0.0710 0.0710 -0.0679SCAN 0.0192 -0.0152 -0.0217SCAN-SIC 0.1474 0.1474 -0.0943SCAN-sdSIC 0.0334 0.0334 -0.0446 that comprises atomization energies of six representativemolecules. Since the reference atomization energies rangefrom 101 kcal/mol to 1149 kcal/mol, it is reasonable tocompare percentage error rather than absolute error.Percentage error = calculated value − reference valuereference value × TABLE IV:
ME in kcal/mol, MAE in kcal/mol and MAPEfor atomization energy of the AE6 set for various levels ofapproximation.Approx. ME MAE MAPELSDA 75.5 75.5 16.7LSDA-SIC 53.5 57.8 9.9LSDA-sdSIC 24.5 25.7 5.6LSDA-LSIC -0.9 9.3 3.0LSDA-LSIC+ -0.4 8.2 2.3PBE 10.6 13.8 3.8PBE-SIC -15.6 17.8 5.5PBE-sdSIC 6.8 11.7 4.1SCAN 0.3 3.0 1.4SCAN-SIC -24.4 26.1 7.0SCAN-sdSIC -3.3 5.7 2.4
LSDA-LSIC+ works better than the LSDA, PBE, theirSIC counterparts, and even LSDA-LSIC. The MAPE ofLSDA-LSIC+ is close to SCAN and SCAN-sdSIC. Fig.5 shows the comparison of percentage errors for SCAN,SCAN-SIC, LSDA-LSIC, and LSDA-LSIC+ for the indi-vidual molecules. Percentage errors are calculated us-ing Eq. 14. The reference values for calculating theerrors for the AE6 set are taken from Ref. 48. SCANoverestimates the atomization energies slightly on aver-age, SCAN-SIC over-corrects and hence underestimates
FIG. 5:
Percentage errors in atomization energies forindividual molecules of the AE6 set by SCAN, SCAN-SIC,LSDA-LSIC, and LSDA-LSIC+. them. This is because, as we proceed from separatedatoms to the molecules they form, the valence orbitals ac-quire more orbitals with which they overlap and, hence,become more noded. Thus, SIC to SCAN makes theenergy of the molecule higher than it would have beenif the orbital densities had been nodeless. SCAN-sdSICscales down the SCAN-SIC results and yields an accuracynear that of SCAN. The results for the total ground stateenergy of individual atoms and molecules using LSDA-LSIC, LSDA-LSIC+, SCAN, SCAN-sdSIC, and SCAN-SIC are presented in the supplemental information.
C. Chemical barrier heights
The barrier height of a chemical reaction is the dif-ference between the maximum energy of the transitionstate and the total energy of the reactants (forward bar-rier) or the products (reverse barrier). The transitionstate involves stretched bonds where SIE is expected tobe maximal. LSDA-LSIC and LSDA-LSIC+ both per-form almost equally well for the representative set of bar-rier heights BH6. The MAE of both the methods is bet-ter than that of SCAN-SIC and SCAN-sdSIC. Table Vpresents the MEs and MAEs of barrier heights for differ-ent approximations. The reference values are taken fromRef. 49.
D. Ionization potential
The ionization potential of a system refers to the en-ergy required to remove the most loosely bound elec-tron. The ionization potential can be calculated usingeither the vertical or the adiabatic method. Here, theadiabatic ionization potential of a system refers to thedifference in energy of the neutral system, at its moststable geometry, and of the cationic system, at its sta-ble geometry, whereas vertical ionization potential refersto the difference of the neutral and cationic system in5hattarai et al.
Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correction
TABLE V:
Mean error (ME) and mean absolute error(MAE) of barrier heights for the BH6 set for various levelsof approximation. The values are in kcal/mol.Approx. ME MAELSDA -18.1 18.1LSDA-SIC -5.1 5.1LSDA-sdSIC -4.1 4.1LSDA-LSIC 0.6 1.4LSDA-LSIC+ -0.2 2.1PBE -9.6 9.6PBE-SIC 0 4.6PBE-sdSIC -3.7 4.2SCAN -7.9 7.9SCAN-SIC -1 3SCAN-sdSIC -4.6 4.6 the stable geometry of the neutral. This means that,in the vertical method, we assume that the geometry ofthe system remains unchanged even after the removalof an electron. In the ∆
SCF method, the ionizationpotential is obtained as the difference between the to-tal ground state energy of the neutral and the cationicsystem. The vertical ionization potential can also be ac-cessed by using Janak’s theorem, i.e., the derivative oftotal energy with respect to orbital occupation is equal tothe eigenvalue of that orbital ( ∂E∂f i = ε i ), independent ofthe detailed form of the exchange-correlation functional.This implies that the exact vertical ionization potentialis simply the negative of the exact KS HOMO eigen-value of the neutral system. Table VI shows the adi-abatic ionization potential of the G2-1 set.
For theadiabatic ionization potential with the ∆
SCF method,LSDA-LSIC+ improves over LSDA-LSIC and is almostcompetitive with SCAN and SCAN-sdSIC. The verticalionization potential result of the same test set accessedusing Janak’s theorem is presented in Ref. 38. For aset of intermediate-sized organic molecules , the error inthe highest-occupied orbital energy is 0.83 eV for LSDA-LSIC+ and 1.13 eV for LSDA-LSIC. E. Electron affinity
The electron affinity of a neutral system is defined asthe change in energy when an extra electron is addedin its isolated gaseous phase. Tables VII and VIII showthe MEs and MAEs of the electron affinities of the G2-1 set using vertical and adiabatic methods, respec-tively. The adiabatic electron affinity assumes that thegeometry of the system is relaxed after the addition ofan electron. In the vertical method, we assume that thegeometry of the system remains unchanged even afteradding an extra electron. In the ∆
SCF method, elec-tron affinity is obtained as the difference between thetotal ground state energy of the neutral and anionic sys-
TABLE VI:
Mean error (ME), and mean absolute error(MAE) of the adiabatic ionization potential of the G2-1 setusing the ∆SCF method for various levels of approximation.The reference values are taken from Ref. 55. The values arein kcal/mol.Method ME MAELSDA 2.7 6.2LSDA-SIC 7.0 10.3LSDA-sdSIC 3.6 7.8LSDA-LSIC 2.7 7.1LSDA-LSIC+ 1.4 6.4PBE -0.8 4.8PBE-SIC -4.7 10.6PBE-sdSIC 0.8 6.8SCAN -0.9 5.8SCAN-SIC -5.4 8.3SCAN-sdSIC -1.2 5.5 tems. For an exact functional, the negative of the eigen-value of the HOMO of an anionic system yields the ver-tical electron affinity.
We have accessed the verticalelectron affinity of the G2-1 set with this method. Semi-local DFAs fail to capture an extra electron in the an-ionic system. Hence, LSDA, PBE, and SCAN results pre-sented here are single step non-SCF calculations with theSCAN-SIC density. LSIC and LSIC+ have equal MAEsfor the adiabatic electron affinity and perform better thanSCAN and SCAN-sdSIC. However, LSIC+ seems to besuperior to LSIC for the vertical electron affinity of theG2-1 set. As single step non-SCF calculations do notproduce eigenvalues, the quasi-SCF versions of LSIC andLSIC+ are used to calculate vertical electron affinity.
TABLE VII:
MEs and MAEs of the vertical electronaffinity of the G2-1 set, calculated by using the eigenvalue ofthe HOMO. All values are in kcal/mol.Method ME MAELSDA -75.0 75.0LSDA-SIC 30.2 30.2LSDA-LSIC 6.5 8.9LSDA-LSIC+ -1.4 6.3
F. Equilibrium bond length
Table IX shows the MEs and MAEs of the equilib-rium bondlengths of a benchmark set of 11 diatomicmolecules , evaluated by finding the minimum of thequadratic fit of the total energy over varying bond length.SCAN performs much better than any other functional.LSIC and LSIC+ values lie in between SCAN and SCAN-SIC, with LSIC having a slightly better result than theLSIC+ and closer to SCAN-sdSIC results.6hattarai et al. Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correction
TABLE VIII:
MEs and MAEs of the adiabatic electronaffinity of the G2-1 set, calculated by the ∆
SCF method.All values are in kcal/mol. The reference values are takenfrom Ref. 53.Method ME MAELSDA 5.7 5.9LSDA-SIC -1.4 5.6LSDA-sdSIC -0.3 3.4LSDA-LSIC 1.5 3.2LSDA-LSIC+ -0.2 3.2PBE 1.3 2.0PBE-SIC -13.2 13.4PBE-sdSIC -4.5 4.6SCAN -0.3 4.1SCAN-SIC -9.0 9.4SCAN-sdSIC -2.8 5.2
TABLE IX:
MEs and MAEs of the equilibrium bondlengths of a benchmark set of 11 diatomic molecules (inangstrom). Results other than those for LSIC+ are the sameas in Table VIII of Ref. 10, where the large improvementdue to interior-scaling of the SIC was first found.Method ME MAELSDA 0.0076 0.0110LSDA-SIC -0.0317 0.0392LSDA-sdSIC -0.0085 0.0189LSDA-LSIC -0.0015 0.0129LSDA-LSIC+ -0.0024 0.0147PBE 0.0123 0.0123PBE-SIC -0.0134 0.0257PBE-sdSIC -0.0019 0.0132SCAN 0.0039 0.0057SCAN-SIC -0.0190 0.0197SCAN-sdSIC -0.0110 0.0110 G. Dissociation energy of noncovalently bondedcomplexes
The S22 dataset consists of the dissociation energiesof complexes built up from pairs of small to large, mostlyorganic molecules, at reference geometries. These com-plexes are divided into three categories: hydrogen bondedcomplexes, complexes with predominant dispersion inter-actions, and mixed complexes (with both kinds of bonds).We tested the performance of different approximationson this dataset. Table X presents the MEs and MAEsin kcal/mol for these functionals. All of our calculationsin Table X are spin unpolarized and self-consistent, withthe exception of LSDA-LSIC and LSDA-LSIC+ (evalu-ated on LSDA-SIC total and localized orbital densities)and SCAN-sdSIC (evaluated on SCAN-SIC total and lo- calized orbital densities). While the other functionalsperform reasonably well, both interior-scaled function-als, LSDA-LSIC and LSDA-LSIC+, consistently under-estimate the binding of the complexes, in several casespredicting them to be unbound. Detailed results for theindividual molecules are presented in the supplementalinformation. Until we can implement LSDA-LSIC andLSDA-LSIC+ self-consistently, we cannot say whethertheir failures for weak bonds might be corrected by fullself-consistency. All we can say is that the LSDA-SICand SCAN-SIC densities and FLO densities on whichthey are evaluated are not by themselves wrong, sincethe LSDA-SIC and SCAN-SIC S22 binding energies arereasonable.It might be surprising that the original PZ SIC (la-belled SIC in Table X) is not worse than it is for the S22set, since the self-interaction corrections from the valenceorbitals are large in comparison with the weak bindingenergies of the S22 complexes. The explanation for thislies in the chemical interpretation of the FLOs as covalentdouble bonds, covalent single bonds, lone pairs, and coreorbitals. These orbitals are often changed only modestlywhen weak bonds are broken. Thus most of the self-interaction correction from the valence electrons cancelsout of the binding energies of the complexes. The can-cellation is even more perfect when the weak bonds arenot all broken but are simply re-arranged, as in the struc-tural energy differences among the water hexamers . Anexception to the partial cancellation upon weak-bondbreaking is the hydrogen bond in H O + ( H O ), which isstrong enough to induce a transition of electrons betweencovalent-bond orbitals and lone-pair orbitals.The cancellation in PZ SIC discussed in the previ-ous paragraph can be lost in LSDA-LSIC (evaluated onLSDA-SIC densities and orbitals). The FLO densities arewhat they were in LSDA-SIC, but their contributions tothe self-interaction correction are now scaled down, andthe scaling can be stronger in the overlap region thanin the corresponding regions of the same monomer orisolated molecule. This arises because of a reduction inthe scaling parameters z σ and f ( z σ ) in the region wherethe charge from the two components in a given complexoverlaps.In ( H O ) , for example, the sum of the orbital self-interaction corrections to the total energy differs by only0.2 mHa from the sum of the corresponding SIC orbitalenergies in the two isolated molecules. In the LSIC cal-culation, however, this difference is 10.0 mHa (about 6kcal/mol). This difference can be traced to contributionsfrom FLOs corresponding to the O lone-pair orbital andthe O-H bond orbital that are located in the overlap re-gion between the two water molecules. Fig. 6 showsthe comparison of the total electron density n and thescaling factor z σ of the water dimer along the bond axisof O-H in the overlapping region and the correspondingnon-overlapping region. The plot omits the very largedensity near the O atom nucleus. It can be seen that thevalue of z σ as well as n remains similar up to the position7hattarai et al. Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correctionof the H atom at about 0.96 ˚A. The density remains al-most equal until about 1.5 ˚A but z σ differs significantly.After about 1.5 ˚A, the density in the overlapping regionincreases because of the presence of the lone pair electronof another water molecule. The scaling factor z σ ( f ( z σ )for LSIC+) is seen to be smaller in the overlap regionthan at corresponding locations outside the overlap re-gion or in the isolated molecules. Because of this, theLSIC energies for these orbitals are relatively smaller inmagnitude than the LSIC energies for the correspondingorbitals in the isolated molecules, resulting in the en-ergy difference noted above. Because the LSDA-SIC andLSDA-LSIC orbital energies are negative, and the energydensities that get scaled down are primarily negative, theimplication is that the complex is relatively less stable inLSDA-LSIC than in LSDA-SIC. Similar effects are foundfor other complexes we have analyzed. FIG. 6:
Comparison of the total electron density n and thescaling factor z σ in the overlapping region (OR) and acorresponding non-overlapping region (NOR) of the waterdimer. The position of the O atom is at the origin for eachcurve. In the inset is a diagram that shows the positions ofthe O (red) and H (white) atoms, as well as the FODpositions (pink). The scaling factor and the total electrondensity are plotted along the dotted lines shown on theinset. The position co-ordinates are in ˚A and the totalelectron densities are in atomic units. Why is the scaling factor smaller in the overlap region? z σ is close to one in one-electron regions of the chargedensity. Far from the center of an isolated molecule,the density tends to be dominated by contributions froma single orbital and z σ is correspondingly close to one.In the overlapping charge region of the molecular com-plex, however, the density decreases to a minimum mov-ing away from one molecule and then increases again onnearing the second molecule. Since z σ depends on thegradient of the density (see Eq. 7), the scaling factormust be close to zero near this minimum.Performing self-consistent LSDA-LSIC and LSDA-LSIC+ calculations could improve the results for the S22 interaction energies. Shifting the electron density in theoverlap region of the complexes slightly could have a sig-nificant impact on the scaling parameters z σ and f ( z σ )without changing other quantities that affect the total en-ergy. Alternatively, employing LSDA-LSIC and LSDA-LSIC+ -type approaches with scaling functions based onimproved iso-orbital indicators α or β could also im-prove the results for the interaction energies. Unlike z σ ,these indicators can distinguish densities associated withweak (i.e., van der Waals or hydrogen) bonds from slowly-varying densities. FIG. 7:
Individual errors(in kcal/mol) for the dissociationenergy of complexes of the S22 dataset
FIG. 8:
FLO densities of the benzene molecule obtainedusing two FOD arrangements. The sym arrangement yieldsthe FLOs in panel (a) (top view) and (b) and (c) (sideviews). The asym arrangement gives the FLOs in panel (d),(e) and (f). FLOs corresponding to C-H bonds (orange),C-C single bonds (blue), and C-C double bonds (red andgreen) are shown. Isosurface values of 0.001 e/ ˚A were usedto make all FLO images. FODs corresponding to each FLOtype are shown by small spheres of the corresponding color.In the side views, only the double bond FLOs are shownalong with the corresponding FOD positions. Another conclusion from Table X and Ref. 10 isthat SCAN-sdSIC, evaluated on SCAN-SIC densities andFLOs, is almost uniformly more accurate than SCAN-LSIC or SCAN-LSIC+ evaluated on LSDA-SIC densitiesand FLOs for equilibrium properties of molecules, and8hattarai et al.
Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correction
TABLE X:
The dissociation energy of weakly boundcomplexes in the S22 dataset. All values are in kcal/mol.The calculations are spin-restricted. The reference values aretaken from Ref. 59. A negative ME indicates that anapproximation underestimates the binding, on average.ME MAELSDA 2.16 2.18LSDA-SIC 1.41 2.98LSDA-LSIC -8.95 9.77LSDA-LSIC+ -7.00 7.42PBE -2.56 2.58PBE-SIC -2.64 2.97SCAN -0.60 0.91SCAN-SIC 2.37 4.16SCAN-sdSIC -1.14 2.46 especially so for the weakly-bound complexes. Figure 7shows the individual errors in the dissociation energy ofthe complexes in the S22 database by SCAN, SCAN-SIC,and SCAN-sdSIC. SCAN-sdSIC overestimates and makesslightly more error than both SCAN and SCAN-SIC forhydrogen-bonded complexes (complexes 1-7). SCAN-sdSIC restores the capacity of SCAN to capture theshort and intermediate range van der Waals interactionsinvolved in the complexes with predominant dispersioncontribution (complexes 8-15) and the mixed complexes(complexes 16-22) which are somewhat degraded fromSCAN to SCAN-SIC. Hence, SCAN-sdSIC might be abetter option to scale down the PZ-SIC and predict theproperties of weakly bonded systems.
H. Do the FLOs have a consistent chemicalinterpretation?
The FLOs are often interpreted chemically as covalentdouble bonds, covalent single bonds, lone pairs, etc. Butthis interpretation is only useful if it can be made con-sistently. The benzene molecule is a challenging examplein this context. There are two distinct solutions for thespin unpolarized benzene monomer in LSDA-SIC. In one( sym ), the FODs are as shown in Fig. 8: for the C-Cbonds around the ring, a single FOD at the center of onebond alternates with a pair of FODs that are symmetri-cally placed above and below the bond center of the next,nominally representing single (C-C) and double (C=C)bonds, respectively. The single bond FOD gives rise tothe FLO density depicted in blue in Fig. 8, while the dou-ble bond FODs give rise to a pair of symmetric FLOs, oneabove the plane of the molecule and the below, its mirrorimage, below it (Fig. 8 (b) and (c)). In the other solu-tion ( asym ), the double bond FODs both lie on the sameside of the molecular plane. In this case, the single-bondFLO is essentially unchanged from the sym solution, but the double bond FLOs change shape dramatically to thered and green FLO densities shown in Fig. 8 (e) and (f).(Note that it is the FODs, and not the FLO densities,that can be strongly asymmetric around the plane of themolecule.) The red FLO is very similar to the single bondFLO and can be thought of as a bonding combination of sp hybrid orbitals. The green FLO can be thought of asa bonding combination of p z orbitals and has a node inthe molecular plane. In the sym solution, the symmetricdouble bond FLOs each have an SIC energy of -0.016 Ha.In the asym solution, the red and green FLOs have SICenergies of -0.048 and +0.023 Ha, respectively. Despitethis, the total energy of the asym solution is only 1.4mHa (0.9 kcal/mol) lower in energy than the sym solu-tion in LSDA-SIC. This remarkable near degeneracy isbroken in LSDA-LSIC, where the ordering switches andthe sym solution is lower by 23.2 mHa (14.6 kcal/mol). Inthe benzene-containing complexes, we only found asym -type LSDA-SIC solutions, despite using sym -type start-ing points for the FODs. The asym and sym solutionsare also nearly degenerate in PBE-SIC.There is also a third nearly degenerate solution forthe benzene monomer if spin polarized densities areconsidered. This solution ( Linnett ) has been discussedelsewhere and can be linked to the Linnett double quar-tet bonding theory, a spin-polarized generalizationof Lewis theory. In this case, sym -type FOD arrange-ments in the up and down spin channels are shifted byone bond around the ring, so that each bond has a singleFOD of one spin at the bond center and two FODs ofthe opposite spin placed symmetrically above and belowthe center. The FLOs for the Linnett solution are verysimilar to the sym
FLOs shown in Fig. 8. The totalenergy of the
Linnett solution in LSDA-SIC is slightlylower (1.3 mHa = 0.8 kcal/mol) than that of the asym solution. For the parallel-displaced benzene dimer (S22complex 11), we find that the
Linnett solution for thedimer is also slightly lower than the asym dimer, so thatthe calculated dissociation energies are nearly identical,1.6 kcal/mol and 1.1 kcal/mol for
Linnett and asym re-spectively. The data in Table X reflect spin unpolar-ized calculations only. We do not expect that significantdifferences would arise if spin-polarized
Linnett arrange-ments were used for all complexes containing benzene.The way to find the localized orbitals suggested byPerdew and Zunger 1981, and universally applied sincethen, is to choose the unitary transformation of thecanonical or generalized Kohn-Sham orbitals that min-imizes the self-interaction-corrected total energy. Butthat approach is inconsistent with the chemical interpre-tation of the FLOs. To preserve a consistent chemicalinterpretation, we need a new way to select the opti-mum unitary transformation, for example, by minimiz-ing a measure of the inhomogeneity of the FLO densi-ties. Such a measure should penalize all strongly nodedor lobed orbital densities such as the green ones in Fig.8 (d) and (e).9hattarai et al.
Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correctionThe nuclear geometry used here for benzene is un-changed under 60-degree rotations around a central axis.There must be an exact ground-state electron densitythat preserves this symmetry. The single-determinantauxiliary wave functions that we have calculated breakthis symmetry. In either the spin unpolarized sym / asym or the Linnett scheme, there are two degenerate Slaterdeterminants with different densities ( sym / asym ) or spindensities ( Linnett ), the density or spin density of onedeterminant transforming into that of the other under a60-degree rotation. A possible physical interpretation is that there are low-frequency fluctuations of the densityand spin-density that anti-correlate on alternate bondsaround the benzene ring. IV. Summary and conclusions
We introduced LSIC+ that takes the iso-orbital indi-cator z σ as its argument for interior scaling of Perdew-Zunger self-interaction correction. LSIC+ is applied withLSDA only and not with GGAs and meta-GGAs to avoidgauge-inconsistency. LSIC+ is an improved version ofLSIC that satisfies one more constraint from the large- Z asymptotic expansion of exchange-correlation energy,in addition to restoring the correct uniform density limitof the exchange energy violated by PZ-SIC. The relevanceof the uniform-density and large-Z limits to valence-electron properties was explained in Ref. 32. LSIC+ re-tains the full PZ-SIC in the region of one-electron density,giving no correction for uniform electron density, and ascaled correction in the region in between. Hence, LSDA-LSIC+ works well for both stretched bond cases and equi-librium properties that do not involve weak bonds. Wepresent here the results of LSDA-LSIC+ for ground stateenergies of atoms from Z =1-18, atomization energies ofthe AE6 set, chemical barrier heights of the BH6 set,ionization potentials of the G2-1 set, electron affinitiesof the G2-1 set, equilibrium bond lengths, and interac-tion energies of the S22 set of weakly-bonded complexes,and compare with LSDA-LSIC and other SIC-correctedand uncorrected DFAs. LSDA-LSIC+ performs betterthan LSDA-LSIC for most of the properties and compa-rably for a few. LSDA-LSIC+ is even better than SCANor SCAN-SIC for a few properties such as total groundstate energy of atoms and electron affinity. LSIC+ sat-isfies one more constraint than LSIC and produces bet-ter results than LSIC, which shows the significance ofconstraint satisfaction in the construction of functionals.The better performance of SCAN-sdSIC can also be at-tributed to the recovery of higher-order coefficients inthe asymptotic expansion of exchange-correlation energyfor atoms of large Z . However, since LSIC and LSIC+employ functions of z σ , they cannot recognize weakerbonds and (at least when evaluated on LSDA-SIC totaland localized orbital densities)produce poor results forthe van der Waals or hydrogen-bond binding energies.Surprisingly, the Fermi-L¨owdin orbitals that describemultiple covalent bonds are not unique. Two different versions give nearly degenerate LSDA-SIC energies inbenzene-like systems, though not in LSDA-LSIC. A morephysical way to determine the Fermi orbital descriptors isneeded to give a single, consistent picture of these FLOs.Unlike LSIC and LSIC+, SCAN-sdSIC is applied tothe more accurate SCAN functional, which is alreadymuch better than LSDA and hence improves its predic-tive power. It provides good results for several groundstate properties discussed here, including the dissociationenergy of weakly bonded systems. The SCAN-sdSIC, ifscaled down by functions of improved iso-orbital indi-cators α or β , may further improve the results. How-ever, since SCAN-sdSIC provides an incorrect asymptote − X HO r for the exchange-correlation potential, it does notguarantee a good description of charge transfer. The op-timal SIC that remains to be developed might be PZ SICevaluated on complex Fermi-L¨owdin orbitals (withnodeless orbital densities) and with Fermi orbital de-scriptors chosen to minimize a measure of the inhomo-geneity of the orbital densities. Acknowledgments
This work was supported by the U.S. Department ofEnergy, Office of Science, Office of Basic Energy Sciences,as part of the Computational Chemical Sciences Programunder Award No. DE-SC0018331. The work of P.B.and KW. was supported by the U.S. National ScienceFoundation under Grant No. DMR-1939528. This re-search includes calculations carried out on HPC resourcessupported in part by the National Science Foundationthrough major research instrumentation grant number1625061 and by the US Army Research Laboratory un-der contract number W911NF-16-2-0189. A few calcu-lations were carried out on the EFRC cluster supportedby the Department of Energy (DOE), Office of Science(OS), Basic Energy Sciences (BES) through Grant No.DE-SC0012575 to the Energy Frontier Research Center:Center for Complex Materials from First Principles. Theplots are generated by MATPLOTLIB. V. Supplemental material
The detailed results of ground state energy of atoms( Z =1-18), the energy of individual atoms, and moleculesof the AE6 set and the S22 dataset can be accessedthrough supplementary material. Other data in detailthat support the findings of the study are available fromthe corresponding author upon reasonable request. W. Kohn and L. J. Sham, Physical Review , A1133 (1965). J. P. Perdew and K. Schmidt, AIP Conference Proceedings ,1 (2001). J. Sun, A. Ruzsinszky, and J. P. Perdew, Physical Review Letters , 036402 (2015). et al.
Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correction J. P. Perdew, K. Burke, and M. Ernzerhof, Physical ReviewLetters , 3865 (1996). J. P. Perdew and Y. Wang, Physical Review B , 13244 (1992). J. P. Perdew and A. Zunger, Physical Review B , 5048 (1981). C. Shahi, P. Bhattarai, K. Wagle, B. Santra, S. Schwalbe,T. Hahn, J. Kortus, K. A. Jackson, J. E. Peralta, K. Trepte,S. Lehtola, N. K. Nepal, H. Myneni, B. Neupane, S. Adhikari,A. Ruzsinszky, Y. Yamamoto, T. Baruah, R. R. Zope, and J. P.Perdew, The Journal of Chemical Physics , 174102 (2019). M. R. Pederson, A. Ruzsinszky, and J. P. Perdew, The Journalof Chemical Physics , 121103 (2014). D. Hofmann, S. Kl¨upfel, P. Kl¨upfel, and S. K¨ummel, PhysicalReview A , 062514 (2012). P. Bhattarai, K. Wagle, C. Shahi, Y. Yamamoto, S. Romero,B. Santra, R. R. Zope, J. E. Peralta, K. A. Jackson, and J. P.Perdew, The Journal of Chemical Physics , 214109 (2020). Z.-h. Yang, M. R. Pederson, and J. P. Perdew, Physical ReviewA , 052505 (2017). M. R. Pederson, The Journal of Chemical Physics , 064112(2015). K. Sharkas, L. Li, K. Trepte, K. P. Withanage, R. P. Joshi, R. R.Zope, T. Baruah, J. K. Johnson, K. A. Jackson, and J. E. Per-alta, The Journal of Physical Chemistry A , 9307 (2018). R. P. Joshi, K. Trepte, K. P. Withanage, K. Sharkas, Y. Ya-mamoto, L. Basurto, R. R. Zope, T. Baruah, K. A. Jackson,and J. E. Peralta, The Journal of Chemical Physics , 164101(2018). K. P. K. Withanage, S. Akter, C. Shahi, R. P. Joshi, C. Diaz,Y. Yamamoto, R. Zope, T. Baruah, J. P. Perdew, J. E. Peralta,and K. A. Jackson, Physical Review A , 012505 (2019). A. I. Johnson, K. P. Withanage, K. Sharkas, Y. Yamamoto,T. Baruah, R. R. Zope, J. E. Peralta, and K. A. Jackson, TheJournal of Chemical Physics , 174106 (2019). K. Sharkas, K. Wagle, B. Santra, S. Akter, R. R. Zope,T. Baruah, K. A. Jackson, J. P. Perdew, and J. E. Peralta,Proceedings of the National Academy of Sciences , 11283(2020). K. P. K. Withanage, P. Bhattarai, J. E. Peralta, R. R. Zope,T. Baruah, J. P. Perdew, and K. A. Jackson, accepted in TheJournal of Chemical Physics. B. Santra and J. P. Perdew, The Journal of Chemical Physics , 174106 (2019). O. A. Vydrov, G. E. Scuseria, J. P. Perdew, A. Ruzsinszky,and G. I. Csonka, The Journal of Chemical Physics , 094108(2006). O. A. Vydrov and G. E. Scuseria, The Journal of ChemicalPhysics , 191101 (2006). A. Ruzsinszky, J. P. Perdew, G. I. Csonka, O. A. Vydrov, andG. E. Scuseria, The Journal of Chemical Physics , 104102(2007). A. Ruzsinszky, J. P. Perdew, G. I. Csonka, O. A. Vydrov, andG. E. Scuseria, The Journal of Chemical Physics , 194112(2006). S. Kl¨upfel, P. Kl¨upfel, and H. J´onsson, The Journal of ChemicalPhysics , 124102 (2012). R. R. Zope, Y. Yamamoto, C. M. Diaz, T. Baruah, J. E. Peralta,K. A. Jackson, B. Santra, and J. P. Perdew, The Journal ofChemical Physics , 214108 (2019). Y. Yamamoto, S. Romero, T. Baruah, and R. R. Zope, TheJournal of Chemical Physics , 174112 (2020). L. Li, K. Trepte, K. A. Jackson, and J. K. Johnson, The Journalof Physical Chemistry A , 8223 (2020). K. Burke, A. Cancio, T. Gould, and S. Pittalis, The Journal ofChemical Physics , 054112 (2016). J. Schwinger, Physical Review A , 2353 (1981). P. Elliott and K. Burke, Canadian Journal of Chemistry , 1485(2009). A. Cancio, G. P. Chen, B. T. Krull, and K. Burke, The Journalof Chemical Physics , 084116 (2018). A. D. Kaplan, B. Santra, P. Bhattarai, K. Wagle, S. T. u. R.Chowdhury, P. Bhetwal, J. Yu, H. Tang, K. Burke, M. Levy,and J. P. Perdew, The Journal of Chemical Physics , 074114(2020). FLOSIC 0.2 developed by R. R. Zope, T. Baruah, J. E. Peralta,and K. A. Jackson, . D. Porezag and M. R. Pederson, Physical Review A , 2840(1999). P. Jureˇcka, J. ˇSponer, J. ˇCern`y, and P. Hobza, Physical Chem-istry Chemical Physics , 1985 (2006). A. J. Sadlej, Theoretica Chimica Acta , 45 (1991). S. Akter, Y. Yamamoto, C. M. Diaz, K. A. Jackson, R. R. Zope,and T. Baruah, The Journal of Chemical Physics , 164304(2020). S. Adhikari, B. Santra, S. Ruan, P. Bhattarai, N. K. Nepal, K. A.Jackson, and A. Ruzsinszky, The Journal of Chemical Physics , 184303 (2020). A. D. Becke, Physical Review A , 3098 (1988). E. B. Isaacs and C. Wolverton, Physical Review Materials ,063801 (2018). J. H. Yang, D. A. Kitchaev, and G. Ceder, Physical Review B , 035132 (2019). Y. Zhang, J. W. Furness, B. Xiao, and J. Sun, The Journal ofChemical Physics , 014105 (2019). N. K. Nepal, L. Yu, Q. Yan, and A. Ruzsinszky, Physical ReviewMaterials , 073601 (2019). S. Adhikari, H. Tang, B. Neupane, A. Ruzsinszky, and G. I.Csonka, Physical Review Materials , 025005 (2020). N. K. Nepal, A. Ruzsinszky, and J. E. Bates, Physical ReviewB , 115140 (2018). Y. Yamamoto, C. M. Diaz, L. Basurto, K. A. Jackson, T. Baruah,and R. R. Zope, The Journal of Chemical Physics , 154105(2019). S. J. Chakravorty, S. R. Gwaltney, E. R. Davidson, F. A. Parpia,and C. F. Fischer, Physical Review A , 3649 (1993). B. J. Lynch and D. G. Truhlar, The Journal of Physical Chem-istry A , 8996 (2003). L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi, andS. Grimme, Physical Chemistry Chemical Physics , 32184(2017). J. Janak, Physical Review B , 7165 (1978). J. P. Perdew and M. Levy, Physical Review B , 16021 (1997). J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, PhysicalReview Letters , 1691 (1982). L. A. Curtiss, K. Raghavachari, G. W. Trucks, and J. A. Pople,The Journal of Chemical Physics , 7221 (1991). L. A. Curtiss, K. Raghavachari, P. C. Redfern, and J. A. Pople,The Journal of Chemical Physics , 1063 (1997). L. A. Curtiss, P. C. Redfern, and K. Raghavachari, The Journalof Chemical Physics , 084108 (2007). K. Wagle, B. Santra, P. Bhattarai, C. Shahi, K. A. Jackson, andJ. P. Perdew, in preparation. J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria,Physical Review Letters , 146401 (2003). J. W. Furness and J. Sun, Physical Review B , 041119 (2019). M. S. Marshall, L. A. Burns, and C. D. Sherrill, The Journal ofChemical Physics , 194102 (2011). S. Schwalbe, K. Trepte, L. Fiedler, A. I. Johnson, J. Kraus,T. Hahn, J. E. Peralta, K. A. Jackson, and J. Kortus, Jour-nal of Computational Chemistry , 2843 (2019). J. Linnett, Nature , 859 (1960). J. Linnett, Journal of the American Chemical Society , 2641(1961). J. P. Perdew, A. Ruzsinszky, N. K. Nepal, and A. D. Kaplan,Proceedings of the National Academy of Sciences (to appear). J. Sun, B. Xiao, Y. Fang, R. Haunschild, P. Hao, A. Ruzsinszky,G. I. Csonka, G. E. Scuseria, and J. P. Perdew, Physical reviewletters , 106401 (2013). S. Kl¨upfel, P. Kl¨upfel, and H. J´onsson, Physical Review A ,050501 (2011). et al. Exploring and enhancing the accuracy of interior-scaled Perdew-Zunger self-interaction correction M.R. Pederson, Carlos M. Diaz, work in progress. J. D. Hunter, Computing in Science Engineering , 90 (2007)., 90 (2007).