Enhancing the efficiency of density functionals with a novel iso-orbital indicator
EEnhancing the efficiency of density functionalswith a novel iso-orbital indicator.
James W. Furness ∗ and Jianwei Sun ∗ Department of Physics and Engineering Physics, Tulane University, New Orleans, LA70118, USA
E-mail: [email protected]; [email protected] a r X i v : . [ phy s i c s . c h e m - ph ] A ug bstract The accuracy and efficiency of a density functional is dependent on the basic in-gredients it uses, and how the ingredients are built into the functional as a whole. Aniso-orbital indicator based on the electron density, its gradients, and the kinetic energydensity, has proven an essential dimensionless variable that allows density functionalsto recognise and correctly treat various types of chemical bonding, both strong andweak. Density functionals constructed around the iso-orbital indicator usually requiredense real-space grids for numerical implementation that deteriorate computationalefficiency, with poor grid convergence compromising the improved accuracy. Here, anovel iso-orbital indicator is proposed based on the same ingredients that retains thecapability to identify the same chemical bonds while significantly relieving the require-ment of dense grids. Furthermore, the novel iso-orbital indicator gives an improvedrecognition for tail regions of electron densities and is constraint-free for the exchange-correlation potential. The novel iso-orbital indicator is therefore expected to be theprime choice for further density functional development.
Graphical TOC Entry τ - τ vW τ + τ UEG β = τ - τ vW τ UEG α = Radius (a ) ∂α∂ n ∂β∂ n Overlap IndicatorWith SmootherFunctional DerivativesImproved NumericalEfficiency ntroduction Density functional theory (DFT) is, in principle, exact for the ground state energy and elec-tron density of a system of electrons under a scalar external potential, conventionally solvedthrough a set of Kohn–Sham (KS) auxiliary single-particle Schr¨odinger-like equations.
Inpractice however, the exchange-correlation energy, an essential but usually small portionof the total energy, must be approximated as a functional of the electron density. Thecomputational efficiency of this scheme and the accuracy of modern exchange-correlationapproximations has resulted in DFT becoming one of the most widely used electronic struc-ture theories and arguably the only practical method for high-throughput discovery of novelmaterials currently available.Exchange-correlation approximations can be broadly categorised by the ingredients usedinto 5 levels of increasing non-locality. The accuracy of an approximation usually increaseswhen more ingredients are included through increased flexibility though this enhancement isoften accompanied by a deterioration of efficiency, especially when non-local information isincluded. It is therefore critical to understand the ingredients of common density functionalapproximations and how they can be utilised to increase the accuracy of a functional whilstmaximising the possible computational efficiency. As functionals of higher levels are usuallydeveloped based on functionals of lower levels, knowledge obtained for the ingredients andtheir combinations in lower-level functionals can be transferred to the development of morecomplex functionals at higher levels.The lowest three levels of efficient semi-local functionals include the local spin densityapproximation (LSDA), the generalised gradient approximation (GGA), and meta-GGA.
LSDA uses only the electron density and recovers the uniform electron gas (UEG)limit. GGAs add the electron density gradient from which two standard dimensionlessvariables are constructed; s = |∇ n | / (2 k F n ) with k F = (3 π n ) / relevant for exchange, and t c = |∇ n | / (2 k s n ) with k s = q k F /π for correlation measure the inhomogeneity of electrondensities at the length scales of local Fermi wavelength 2 π/k F and Thomas screening length3 /k s , respectively. Commonly used meta-GGAs include one more ingredient, the kineticenergy density τ ( r ) = 1 / P occ .i |∇ ϕ i ( r ) | where { ϕ i } are KS orbitals.The inclusion of τ in meta-GGAs naturally arises from the Taylor expansion of the exactspherically averaged exchange hole and provides a simple and straightforward way to makea correlation functional exactly one-electron self-interaction free. The flexibility due to theinclusion of τ improves the accuracy of meta-GGAs over GGAs, by either better fitting toexperimental data empirically or satisfying more exact constraints and appropriatenorms non-empirically. Early attempts and the recent Tao-Mo functional construct non-empirical meta-GGAs through an iso-orbital indicator defined as, z = τ vW τ , (1)where τ vW ( r ) = |∇ n ( r ) | / n ( r ) is the von-Weisz¨acker kinetic energy density that recovers τ in the single orbital limit. Whilst z can identify single orbital densities, z = 0, andslowly varying densities, z ≈
1, it is unable to distinguish slowly varying densities from thenon-covalent closed shell overlap densities important in intermediate-range van-der-Waalsbonding. A different indicator widely used in empirical constructions is, t − = ττ UEG , (2)where τ UEG = (3 / π ) / n / is the kinetic energy density of uniform electron gas. Whilstthe t − indicator has been shown to differentiate covalent from non-covalent bonding itcannot uniquely identify single-orbital regions for which t − = 5 s / . This is likely oneof the major reasons for overfitting in the M06L meta-GGA and the resulting numericalstability problem. t − is semi-infinite, [0 , ∞ ], and usually mapped to a finite domain via w = (1 − t − ) / (1 + t − ) = ( τ UEG − τ ) / ( τ UEG + τ ).4 further iso-orbital indicator has been constructed, α = τ − τ vW τ UEG , (3)which is able to uniquely identify single orbital, slowly varying and non-covalent overlap den-sities as α = 0, ≈ (cid:29) z in earliermeta-GGA functionals to enforce the correct fourth-order gradient expansion, thoughits importance for characterising bonding densities was not fully recognised until later. The α variable is directly related to the electron localization function (ELF) of Refs. as f ELF = 1 / (1 + α ) which has been used to give a rigorous topological classification ofchemical bonding. In addition, recent α dependent functionals have been found toeffectively handle properties that have traditionally been challenging for semi-local func-tionals, including the intermediate-range van-der-Waals bonding and metal-insulatortransitions. Despite the general success enjoyed by recent meta-GGA functionals for a wide rangeof systems, it has been observed that many meta-GGAs suffer numerical in-stabilities in self-consistent field (SCF) calculations, and have an unacceptably slowconvergence with respect to the density of the numerical integration points. The increasedcomputational cost of dense numerical grids severely limits the usefulness of many meta-GGAfunctionals and restricts the complexity of systems to which they can be applied. Addition-ally, the uncertainty in overall grid convergence necessitates a time-consuming validationthat calculated properties are properly converged in grid density, placing an undesirableburden of expertise on the user.Within α -based functionals, it is understood that this numerical sensitivity originatesfrom sharp oscillations in the exchange-correlation potential that are only properly capturedby very fine grids, particularly in inter-shell regions where the local orbital overlap charactercan vary rapidly. Here we show that the rapid variations in the derivatives of α are largely5esponsible for these undesirable oscillations. Further, we propose a modified iso-orbitalindicator quantity, termed β , from which new functionals can be constructed that do notsuffer this problem, β ( r ) = τ ( r ) − τ vW ( r ) τ ( r ) + τ UEG ( r ) (4)= α ( r ) τ UEG ( r ) τ ( r ) + τ UEG ( r ) ! . (5)The β variable contains similar information about local orbital overlap environments to α whilst having smoother derivatives more easily amenable to evaluation on numerical grids.This allows local orbital overlap information to be used in functionals at the meta-GGA leveland higher without suffering the numerical problems associated with analogous functions ofthe α variable. Table 1: Values of common iso-orbital indicators for important chemical envi-ronments.
Region t − z α β Single orbital 5 s / ≈ ≈ ≈ ≈ Non-covalent bonding (cid:29) ≈ (cid:29) (cid:28) β < α and β variables can be shown by examining theirlimits in the three important orbital overlap regions. Firstly, in single orbital regions τ ( r ) = τ vW ( r ) and α ( r ) = β ( r ) = 0. This bound is important in exchange-correlation functionaldevelopment, allowing the strong lower bound on exchange energy as well as correlationenergy to be enforced for all spin-unpolarised single-orbital systems. Secondly, in slowlyvarying densities τ vW ( r ) → τ ( r ) ≈ τ UEG ( r ), so α ( r ) ≈ β ( r ) ≈ /
2. Finally,in non-covalent density overlap regions where n ( r ) → τ UEG ( r ) → n ( r ) / while τ ( r ) − τ vW ( r ) → n ( r ), and thus the denominator, τ UEG ( r ), decays faster than thenumerator, τ ( r ) − τ vW ( r ), leading to α ( r ) → ∞ . Here β ( r ) approaches 1 however, sinceboth its numerator and denominator decay as n ( r ). The relations between α , β and other6so-orbital indicators are summarised in Table 1. r ( a )0123 a) Li Maj. Min. Maj. Min. r ( a )4002000200400 c) Li d / dnd / dn d / dnd / dn d / dn Maj. d / dn Min. d / dn Maj. d / dn Min. r ( a )0123 b) C Maj. Min. Maj. Min. r ( a )10050050100 d) C d / dnd / dn d / dnd / dn d / dn Maj. d / dn Min. d / dn Maj. d / dn Min.
Figure 1:
Radial plots of α ( r ) and β ( r ) functions and their density derivativesfor the lithium and carbon atoms. The α (solid red) and β (dashed blue) functions anddensity derivatives are plotted for the majority (Maj., (cid:78) ) and minority(Min., (cid:72) ) spin fromself-consistent PBE densities. Radial distances are in units of Bohr.This similarity between α and β in highlighting local orbital overlap environment isshown graphically in Figures 1 a) and b) for the lithium and carbon atoms respectively.The lithium and carbon atoms were chosen as typically difficult systems with pronouncednumerical sensitivity for density functionals. The similar inter-shell peak behaviour is clearlyvisible for both variables as a rise away from single orbital like values of the 1 s core. Thedifference between α and β is seen in the tail region of the carbon atom for the majorityspin electrons, where α diverges upwards whilst β decays slowly towards 0. For the minorityspin, both α and β approach 0 in the tail region. In general, β consistently indicates the tailregions of all densities with 0 as τ vW ( r ) is the leading order of τ ( r ) there, and identifies tail7egions completely when combined with the dimensionless reduced density gradient, s , thatis divergent. This is impossible for α which can have different values of α ( r ) = 0 for the tailregions of single orbital systems and α → ∞ otherwise.The benefit of β over α for functional development is most clearly seen in the derivativeswith respect to density, exemplified in Figure 1 c) and d) for the lithium and carbon atomsrespectively. The derivative of β with respect to density does not show the same rapidvariation as that of α in valence regions where density is significant. Similar behaviouris observed for |∇ n | and τ derivatives, as noted for α in Ref. and in the more complexmolecular systems presented in supplementary material Section 1.These much better behaved derivatives of β over those of α are expected to remedy thenumerical problems observed in α -dependent meta-GGA functionals. To validate this, wemodify the existing meta-GGA made simple 2 (MS2) functional by substituting 2 β inplace of α within the exchange enhancement factor F x ( s, α ), defined by E meta − GGAx [ n ] = Z e UEGx ( n ) F meta − GGAx ( s, α ) d r , (6)where e UEGx ( n ) = − π ) / n / / π is the exchange energy density of uniform electron gas.MS2 has a simple construction of F MS2x ( s, α ) = F x ( s ) + f MS2 x ( α )[ F x ( s ) − F x ( s )] with f MS2x ( α ) = (1 − α ) α + bα (7)that interpolates between F x ( s ), a GGA for single-orbital systems ( α = 0), and F x ( s ), aGGA for slowly-varying densities ( α ≈ α (cid:29)
1) for noncovalentbonds.This functional was chosen for the relatively simple construction and well reported numer-ical instabilities. Modification of the simple exchange enhancement factor for β dependencewas trivial by replacing α with 2 β to guarantee the interpolation between F x ( s ) for singleorbital systems (2 β = α = 0) and F x ( s ) for slowly varying densities (2 β ≈ α ≈ able 2: Parameters for the MS2 and MS2 β functionals. The notation of Ref. is followed, with constants k = 0 . and µ GE = 10 / . † The original publicationof MS2 gives c = 0 . . We find this gives a small error (on the order of µE h )in the exchange energy of the hydrogen atom that is corrected using c = 0 . .We note however, that this change has a completely negligible impact on MS2predicted atomisation energies and barrier heights. MS2 κ = 0 . b = 4 . c = 0 . † MS2 β κ = 0 . b = (27 b MS2 − / c = 0 . β -modified meta-GGA functional is termed MS2 β . The parameter b is determinedsuch that f MS2 β x ( β = 1) = f MS2x ( α → ∞ ) for noncovalent bonds. These minor adjustmentsto the balances between internal parameters to preserve exact constraints obeyed by theparent functional are summarised in Table 2 and detailed in supplementary material. MS2 β was implemented into the Turbomole package using the XCFun library to automat-ically evaluate functional derivatives. We present MS2 β simply as a convenient means forpreliminary investigation of the β variable rather than as a viable general functional.
50 100 150 200 250 300Total Radial Grid Points9876543 l o g ( | EE c o n v . | ) a) Li 50 100 150 200 250 300Total Radial Grid Points9876543 l o g ( | EE c o n v . | ) b) C MS2 MS2
Figure 2:
Convergence behaviour with respect to numerical grid density of con-ventional and β -modified meta-GGA made simple 2 (MS2) functional for (a)lithium and (b) carbon atoms. Difference in self-consistent atomic total energy relativeto the total energy from a converged grid of 520 points ( E conv . ) for the MS2 (red, × ) andMS2 β (blue, dashed, +) functionals using the aug-cc-pVQZ basis set. An acceptable con-vergence is assumed when difference remains below the SCF convergence tolerance of 1 µE h (dotted line). 9he grid convergence of the MS2 and MS2 β functionals is shown in Figure 2 as a plotof the difference in self-consistent electronic energy relative to the same calculation using avery fine benchmark grid, as a function of grid point density, for the carbon and lithiumatoms. In both cases the β -modified functional shows a convergence in total energy at lowergrid density than the parent α functional, indicated by the difference to the benchmark gridenergy remaining below the SCF convergence threshold of 1 µE h .As previously noted, the sensitivity of α -dependent functionals to numerical integrationgrid point density is understood as an effect of rapid oscillations in the exchange-correlationpotential, expressed in the generalised Kohn–Sham scheme as, Z ϕ p ˆ v xc ϕ q d r = Z ϕ p ∂e xc ∂n ϕ q d r + Z ∂e xc ∂ ∇ n · ∇ ( ϕ p ϕ q ) d r + 12 Z ∇ ϕ p · ∂e xc ∂τ ! ∇ ϕ q d r (8)Here, e xc is the exchange-correlation energy density of a meta-GGA functional. Hence,functionals showing sharp oscillations in functional derivatives will show sharp variationsin exchange-correlation potential that are challenging to capture in numerical integrationschemes. By modifying meta-GGA functionals to use the β indicator in place of α , oscilla-tions in functional derivatives are minimised and the resulting exchange-correlation potentialis smoothed, as exemplified by MS2 and MS2 β which have the same exchange enhancementfactor form in terms of α or 2 β dependence. We can therefore understand the reduced gridsensitivity of the β -modified functional as a consequence of smoother functional derivativesin the inter-shell and valence regions producing smoother exchange-correlation potentialsthat can be more easily captured by numerical integration grids.The divergences of derivatives at tails seen in Figures 1 c) and d) are not problematic forthe exchange-correlation potential of any analytic β -dependent meta-GGAs because e UEGx ( n )decays faster than the divergence of β derivatives there. However, the divergences of α derivatives at tails are much stronger than those of β , as also seen in Figures 1 c) andd), and can cause problems for the exchange-correlation potential of α -dependent meta-10GAs. At the tail regions, τ decays as n , whilst τ UEG decays more quickly as n / , so thederivative of α with respect to, for example, τ diverges as n − / , faster than the decay of e UEGx ( n ) proportional to n / . This can potentially lead to divergences of exchange-correlationpotential for α -dependent meta-GGAs as long as the derivative of F meta − GGAx ( s, α ) withrespect to α is non-zero at the tail. Such non-zero values are encountered by the meta-GGAmade very simple (MVS) and SCAN meta-GGAs at the tail regions of single-orbitalsystems. This behaviour is avoided in β as the decay of both the numerator and denominatoris determined by τ as decaying with n .The potential divergence of exchange-correlation potential at tail regions resulting fromthe strong divergence of α derivatives is undesirable and problematic especially for the con-struction of pseudo-potentials from isolated atoms, acting as a hitherto unrecognised ad-ditional constraint on α dependent functional design. This constraint was not previouslynoticed and is not obeyed by the MVS and SCAN functionals, though it is fortuitouslyobeyed by MS2. This constraint is not necessary in β dependent functionals however, asthe asymptotic behaviour of β is properly controlled and thus β offers simplicity and greaterflexibility in functional design compared to the α indicator. Table 3: Mean error (ME) and mean absolute error (MAE) in kcal mol − for thesmall test sets of atomisation energy (AE6) and reaction barrier height (BH6)for the MS2 and MS2 β functionals. All calculations were performed fully self-consistently with the 6-311++G(3df,3pd) basis set on a dense numerical grid(Turbomole level 7). MS2 MS2 β AE6 (ME) -0.78 3.95AE6 (MAE) 4.40 6.10BH6 (ME) -5.79 -6.32BH6 (MAE) 5.79 6.32To show that the MS2 β functional retains the ability to effectively distinguish differentchemical environments we calculate the mean error (ME) and mean absolute error (MAE)of small data sets of atomisation energies (AE6) and barrier heights (BH6), the resultsfor which are shown in Table 3. For these small datasets the accuracy of the β -modified11 M A E A E M A E c o n v . A E ( k c a l m o l ) M A E B H M A E c o n v . B H ( k c a l m o l ) MS2M 2
Figure 3: Convergence of mean absolute error for the AE6 and BH6 small molecule datasets, with respect to increasing numerical grid density for the MS2 (red, × ) and MS2 β (blue, dashed, +) functionals. The presence of both angular and radial grid components andthe variety of molecules in the sets prevents a simple quantitative measure of grid density. Assuch, grid density is reported as the chosen Turbomole grid level, such that a higher gridlevel (larger x axis value) directly corresponds to more dense integration grids, though stepsbetween levels may not be uniform. All calculations were performed fully self-consistentlywith the 6-311++G(3df,3pd) basis set. β . We conclude therefore that whilst it is possible to use 2 β directly inplace of α in existing functionals, some adjustment of internal parameters and interpolation-extrapolation functions is likely to be necessary to construct useful general β dependentfunctionals.A comparison of the data set grid convergence for the MS2 and MS2 β functionals issummarised in Figure 3 as a comparison of increasing grid density with respect to the densestgrid. A similar trend to the convergence of atomic total energies, Figure 2, is seen with MS2 β reaching a converged error at coarser numerical grids than the parent MS2 functional. Wenote that whilst convergence for the BH6 set is relatively rapid for both functionals muchgreater variation is seen for the atomisation energies in the AE6 set, consistent with thelarger energy scale of the test set. Å )0.40.30.20.10.00.1 I n t e r a c t i o n E n e r g y ( k c a l m o l ) MS2MS2
Figure 4:
Dissociation curve of the argon dimer.
Calculated using MS2 (red, × ) andMS2 β (blue, +) functionals. Calculations performed self-consistently using the aug-cc-pV5Zbasis set and a very fine numerical grid (grid 7 of the Turbomole program). A benchmarkcurve (black) is included from Ref. The different behaviour of the α and β variables is most pronounced in regions of non-bonded density overlap where α ( r ) (cid:29) . β is in the region 0 . < β ( r ) < .
0. This13ifference is most clearly examined in the Ar dissociation curve, for which it has been shownmeta-GGA functionals can be accurate around equilibrium. The dissociation curvesfor Ar are shown in Figure 4 for both the conventional and β -modified MS2 functionalswith benchmark data included from Ref. for comparison. In contrast to the small testsets, MS2 β closely matches the performance of the original MS2 functional across the wholeof the binding curve, suggesting that the MS2 interpolation-extrapolation function may beas appropriate for 2 β as it is for α in non-covalently bound systems.In conclusion, we have identified the source of the numerical sensitivities commonly suf-fered in calculations employing meta-GGA functionals as resulting from sharp oscillationsin the functional derivatives of the commonly employed dimensionless variable α . We haveaddressed these sensitivities by constructing a related dimensionless variable, β , which im-parts similar information about local orbital overlap environment whilst having smootherfunctional derivatives. The enhanced numerical stability, improved recognition of tail re-gions, and freedom from constraints on exchange-correlation potential presented by β arean appealing opportunity for enhancing the numerical performance of future functionals forboth the meta-GGA level and higher level fully non-local functionals.We have used the simple MS2 meta-GGA functional to show a simple proof of concept bysubstituting 2 β for α in the functional with minimal adjustment of internal parameters. Wefind improved numerical performance in all cases, with the new functional giving convergedproperties from much coarser integration grids than the original MS2 functional. Whilst theMS2 β functional preserves much of the good accuracy of its progenitor, a slight degradationof accuracy is observed for small molecule test sets though performance is preserved for thenon-covalently bound Ar diatomic. This suggests a need for the underlying functional to beadjusted for β rather than naive substitution if good performance for appropriate norms isto be preserved alongside an adherence to exact constraints. We are optimistic that whollynovel functionals utilising the β iso-orbital indicator can provide ever greater accuracy andutility for the wider DFT community. 14 cknowledgements We thank Jefferson Bates for his invaluable technical advice in running and modifying the
Turbomole package. We also thank John Perdew, Adrienn Ruzsinszky, Mark Pederson,Weitao Yang, Arun Bansil, Samuel Trickey, Susi Lehtola and Albert Bartok-Partay for theirfruitful discussions around the ideas presented here. The work at Tulane University wassupported by the start-up funding from Tulane University, and by the U.S. DOE, Office ofScience, Basic Energy Sciences grant number de-sc0235021(core research). Calculations werecarried out on the Cypress cluster at Tulane University, and we thank Dr. Hoang Tran, Dr.Carl Baribault, and Dr. Hideki Fujioka for their computational support
References (1) Hohenberg, P.; Kohn, W. Inhomogeneous electron gas.
Physical Review , , 864–871.(2) Kohn, W.; Sham, L. J. Self-consistent equations including exchange and correlationeffects. Physical Review , , 1133–1139.(3) Perdew, J. P.; Schmidt, K. Jacob’s ladder of density functional approximations for theexchange-correlation energy. AIP Conference Proceedings , , 1–20.(4) Ma, S.-K.; Brueckner, K. A. Correlation Energy of an Electron Gas with a Slowly VaryingHigh Density. Physical Review , , 18–31.(5) Von Barth, U.; Hedin, L. A local exchange-corremation potentiel for the spin polarizedcase: 1. Journal of Physics C: Solid State Physics , , 1629–1642.(6) Perdew, J. P.; Zunger, A. Self-interaction correction to density-functional approximationsfor many-electron systems. Physical Review B , , 5048–5079.157) Perdew, J. P.; Wang, Y. Accurate and Simple Analytic Representation of the Electron-Gas Correlation-Energy. Physical Review B , , 13244–13249.(8) Sun, J.; Perdew, J. P.; Seidl, M. Correlation energy of the uniform electron gas from aninterpolation between high- and low-density limits. Physical Review B , , 085123.(9) Becke, A. D. Density-functional exchange-energy approximation with correct asymptoticbehavior. Physical Review A , , 3098–3100.(10) Perdew, J. P.; Chevary, J.; Vosko, S.; Jackson, K.; Pederson, M.; Singh, D.; Fiolhais, C.Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approx-imation for exchange and correlation. Physical Review B , , 6671–6687.(11) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation MadeSimple. Physical Review Letters , , 3865–3868.(12) Keal, T. W.; Tozer, D. J. A semiempirical generalized gradient approximation exchange-correlation functional. Journal of Chemical Physics , , 5654–5660.(13) Armiento, R.; Mattsson, A. E. Functional designed to include surface effects in self-consistent density functional theory. Physical Review B , , 085108.(14) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Con-stantin, L. A.; Zhou, X.; Burke, K. Generalized gradient approximation for solids andtheir surfaces. Physical Review Letters , , 136406.(15) Constantin, L. A.; Fabiano, E.; Laricchia, S.; Della Sala, F. Semiclassical neutral atomas a reference system in density functional theory. Physical Review Letters , ,186406.(16) Vela, A.; Pacheco-kato, J. C.; G´azquez, J. L.; Campo, J. M.; Trickey, S. B. Improvedconstraint satisfaction in a simple generalized gradient approximation exchange func-16ional Improved constraint satisfaction in a simple generalized gradient. The Journal ofChemical Physics , , 144115.(17) Becke, A. D.; Roussel, M. R. Exchange holes in inhomogeneous systems: A coordinate-space model. Physical Review A , , 3761–3767.(18) Van Voorhis, T.; Scuseria, G. E. A novel form for the exchange-correlation energyfunctional. Journal of Chemical Physics , , 400–410.(19) Becke, A. D. A new inhomogeneity parameter in density-functional theory. Journal ofChemical Physics , , 2092–2098.(20) Perdew, J. P.; Kurth, S.; Zupan, A.; Blaha, P. Accurate Density Functional with CorrectFormal Properties: A Step Beyond the Generalized Gradient Approximation. PhysicalReview Letters , , 2544–2547.(21) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Func-tional Ladder: Non-Empirical Meta-Generalized Gradient Approximation Designed forMolecules and Solids. Physical Review Letters , , 146401.(22) Becke, A. D. A real-space model of nondynamical correlation. Journal of ChemicalPhysics , , 2972–2977.(23) Zhao, Y.; Truhlar, D. G. A new local density functional for main-group thermochem-istry, transition metal bonding, thermochemical kinetics, and noncovalent interactions. Journal of Chemical Physics , , 194101.(24) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Constantin, L. A.; Sun, J. WorkhorseSemilocal Density Functional for Condensed Matter Physics and Quantum Chemistry. Physical Review Letters , , 026403.(25) Peverati, R.; Truhlar, D. G. An improved and broadly accurate local approximationto the exchange-correlation density functional: The MN12-L functional for electronic17tructure calculations in chemistry and physics. Physical Chemistry Chemical Physics , , 13171–13174.(26) Yu, H. S.; He, X.; Truhlar, D. G. MN15-L: A New Local Exchange-Correlation Func-tional for Kohn-Sham Density Functional Theory with Broad Accuracy for Atoms,Molecules, and Solids. Journal of Chemical Theory and Computation , , 1280–1293.(27) Sun, J.; Xiao, B.; Ruzsinszky, A. Communication : Effect of the orbital-overlap de-pendence in the meta generalized gradient approximation. Journal of Chemical Physics , , 051101.(28) Del Campo, J. M.; G´azquez, J. L.; Trickey, S. B.; Vela, A. A new meta-GGA exchangefunctional based on an improved constraint-based GGA. Chemical Physics Letters , , 179–183.(29) Mardirossian, N.; Head-Gordon, M. Mapping the genome of meta-generalized gradientapproximation density functionals: The search for B97M-V. Journal of Chemical Physics , , 074111.(30) Sun, J.; Ruzsinszky, A.; Perdew, J. P. Strongly Constrained and Appropriately NormedSemilocal Density Functional. Physical Review Letters , , 036402.(31) Perdew, J. P.; Constantin, L. A.; Sagvolden, E.; Burke, K. Relevance of the slowlyvarying electron gas to atoms, molecules, and solids. Physical Review Letters , ,1–4.(32) Becke, A. D. Hartree-Fock exchange energy of an inhomogeneous electron gas. Inter-national Journal of Quantum Chemistry , , 1915–1922.(33) Wang, Y.; Jin, X.; Yu, H. S.; Truhlar, D. G.; He, X. Revised M06-L functional forimproved accuracy on chemical reaction barrier heights, noncovalent interactions, and18olid-state physics. Proceedings of the National Academy of Sciences , , 8487–8492.(34) Tao, J.; Mo, Y. Accurate Semilocal Density Functional for Condensed-Matter Physicsand Quantum Chemistry. Physical Review Letters , , 073001.(35) Sun, J.; Xiao, B.; Fang, Y.; Haunschild, R.; Hao, P.; Ruzsinszky, A.; Csonka, G. I.;Scuseria, G. E.; Perdew, J. P. Density functionals that recognize covalent, metallic, andweak bonds. Physical Review Letters , , 106401.(36) Svendsen, P.; von Barth, U. Gradient expansion of the exchange energy from second-order density response theory. Physical Review B , , 17402–17413.(37) Becke, A. D.; Edgecombe, K. E. A simple measure of electron localization in atomicand molecular systems. The Journal of Chemical Physics , , 5397–5403.(38) Silvi, B.; Savin, A. Classification of Chemical-Bonds Based on Topological Analysis ofElectron Localization Functions. Nature , , 683–686.(39) Savin, A.; Silvi, B.; Colonna, F. Topological analysis of the electron localization functionapplied to delocalized bonds. Canadian Journal of Chemistry , , 1088–1096.(40) Savin, A.; Nesper, R.; Wengert, S.; F¨assler, T. F. ELF: The Electron LocalizationFunction. Angewandte Chemie International Edition in English , , 1808–1832.(41) Noury, S.; Colonna, F.; Savin, A.; Silvi, B. Analysis of the delocalization in the topo-logical theory of chemical bond. Journal of Molecular Structure , , 59–68.(42) Paul, A.; Sun, J.; Perdew, J. P.; Waghmare, U. V. Accuracy of first-principles inter-atomic interactions and predictions of ferroelectric phase transitions in perovskite oxides:Energy functional and effective Hamiltonian. Physical Review B , , 054111.1943) Zhang, Y.; Sun, J.; Perdew, J. P.; Wu, X. Comparative first-principles studies of pro-totypical ferroelectric materials by LDA, GGA, and SCAN meta-GGA. Physical ReviewB , , 035143.(44) Zhang, Y.; Kitchaev, D. A.; Yang, J.; Chen, T.; Dacek, S. T.; Sarmiento-P´erez, R. A.;Marques, M. A. L.; Peng, H.; Ceder, G.; Perdew, J. P. et al. Efficient first-principlesprediction of solid stability: Towards chemical accuracy. npj Computational Materials , .(45) Sun, J.; Remsing, R. C.; Zhang, Y.; Sun, Z.; Ruzsinszky, A.; Peng, H.; Yang, Z.;Paul, A.; Waghmare, U.; Wu, X. et al. Accurate first-principles structures and energiesof diversely bonded systems from an efficient density functional. Nature Chemistry , , 831–836.(46) Furness, J. W.; Zhang, Y.; Lane, C.; Buda, I. G.; Barbiellini, B.; Markiewicz, R. S.;Bansil, A.; Sun, J. An accurate first-principles treatment of doping-dependent electronicstructure of high-temperature cuprate superconductors. Communications Physics , , 1–6.(47) Yang, Z.-h.; Peng, H.; Sun, J.; Perdew, J. P. More realistic band gaps from meta-generalized gradient approximations: Only in a generalized Kohn-Sham scheme. PhysicalReview B , , 205205.(48) Johnson, E. R.; Becke, A. D.; Sherrill, C. D.; Dilabio, G. A. Oscillations in meta-generalized-gradient approximation potential energy surfaces for dispersion-bound com-plexes. The Journal of Chemical Physics , , 034111.(49) Mardirossian, N.; Head-Gordon, M. Thirty years of density functional theory in com-putational chemistry: An overview and extensive assessment of 200 density functionals. Molecular Physics , , 2315–2372.2050) Perdew, J. P.; Ruzsinszky, A.; Sun, J.; Burke, K. Gedanken densities and exact con-straints in density functional theory. The Journal of Chemical Physics , ,18A533.(51) Sun, J.; Perdew, J. P.; Ruzsinszky, A. Semilocal density functional obeying a stronglytightened bound for exchange. Proceedings of the National Academy of Sciences , , 685–689.(52) Sun, J.; Haunschild, R.; Xiao, B.; Bulik, I. W.; Scuseria, G. E.; Perdew, J. P. Semilocaland hybrid meta-generalized gradient approximations based on the understanding of thekinetic-energy-density dependence Semilocal and hybrid meta-generalized gradient ap-proximations based on the understanding of the kinetic-energy-density depend. Journalof Chemical Physics , Journal of Chemical Theory and Computation , , 1971–1980.(55) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. Theatoms boron through neon and hydrogen. The Journal of Chemical Physics , ,1007–1023.(56) Neumann, R.; Nobes, R. H.; Handy, N. C. Exchange functionals and potentials. Molec-ular Physics , , 1–36.(57) Seidl, A.; G¨orling, A.; Vogl, P.; Majewski, J.; Levy, M. Generalized Kohn-Sham schemesand the band-gap problem. Physical Review B , , 3764–3774.2158) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. R. Efficient diffusefunction-augmented basis sets for anion calculations. III. The 3-21+G basis set for first-row elements, Li-F. Journal of Computational Chemistry , , 294–301.(59) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-consistent molecular orbital methods25. Supplementary functions for Gaussian basis sets. The Journal of Chemical Physics , , 3265–3269.(60) Lynch, B. J.; Truhlar, D. G. Small Representative Benchmarks for ThermochemicalCalculations. Journal of Physical Chemistry A , , 8996–8999.(61) Haunschild, R.; Klopper, W. Theoretical reference values for the AE6 and BH6 testsets from explicitly correlated coupled-cluster theory. Theoretical Chemistry Accounts , , 1112.(62) Patkowski, K.; Murdachaew, G.; Fou, C. M.; Szalewicz, K. Accurate ab initio potentialfor argon dimer including highly repulsive region. Molecular Physics , , 2031–2045. 22 upplementary Material Behaviour of α and β and their derivatives in a range ofchemical environments The α and β quantities and their d/dn , d/d |∇ n | and d/dτ derivatives are plotted for atomicand molecular systems, chosen to illustrate a range of chemically diverse systems. All eval-uated from self consistent Perdew-Burke-Ernzerhof (PBE) densities. A common scheme of α (blue) and β (red, dotted) for majority ( (cid:78) ) and minority ( (cid:72) )spin channels is used for all plots. r ( a )0123 Maj. Min. Maj. Min. 0 1 2 3 4 5 r ( a )4002000200400 d /dn Maj.d /dn Min.d /dn Maj.d /dn Min.0 1 2 3 4 5 r ( a )100806040200 d /d| n| Maj.d /d| n| Min.d /d| n| Maj.d /d| n| Min. 0 1 2 3 4 5 r ( a )020406080100 d /d Maj.d /d Min.d /d Maj.d /d Min. Figure 1:
Radial plots of the α and β iso-orbital indicators for the lithium atom. r ( a )0123 Maj. Min. Maj. Min. 0 1 2 3 4 5 r ( a )10050050100 d /dn Maj.d /dn Min.d /dn Maj.d /dn Min.0 1 2 3 4 5 r ( a )100806040200 d /d| n| Maj.d /d| n| Min.d /d| n| Maj.d /d| n| Min. 0 1 2 3 4 5 r ( a )020406080100 d /d Maj.d /d Min.d /d Maj.d /d Min. Figure 2:
Radial plots of the α and β iso-orbital indicators for the carbon atom. r ( a )0123 0 1 2 3 4 5 r ( a )201001020 d /dnd /dn0 1 2 3 4 5 r ( a )2010010 d /d| n| d /d| n| r ( a )10010203040 d /dd /d Figure 3:
Radial plots of the α and β iso-orbital indicators for the argon atom. r ( a )0123 Ar 0 2 4 6 8 10 r ( a )402002040 Ar d /dnd /dn0 2 4 6 8 10 r ( a )2010010 Ar d /d| n| d /d| n| r ( a )10010203040 Ar d /dd /d Figure 4:
The α and β iso-orbital indicators for the Ar vand-der-Waals diatomic,plotted along the internuclear axis evaluated for self consistent PBE density. The origin is positioned at the bond center and argon atom position is marked with a dottedvertical line. 26 r ( a )0123 Na 0 2 4 6 8 10 r ( a )3002001000100200300 d f ( r ) / d n Na d /dnd /dn0 2 4 6 8 10 r ( a )2010010 Na d /d| n| d /d| n| r ( a )10010203040 Na d /dd /d Figure 5:
The α and β iso-orbital indicators for the Na diatomic showing metallicbonding, plotted along the internuclear axis. The origin is positioned at the bondcenter and sodium atom position is marked with a dotted vertical line.27 r ( a )0123 NaCl 10 5 0 5 10 r ( a )402002040 NaCld /dnd /dn10 5 0 5 10 r ( a )2010010 NaCl d /d| n| d /d| n|
10 5 0 5 10 r ( a )10010203040 NaCld /dd /d Figure 6:
The α and β iso-orbital indicators for NaCl showing ionic bonding,plotted along the internuclear axis. r ( a )0123 f ( r ) Si H 10 5 0 5 10 r ( a )2001000100200 d f ( r ) / d n Si Hd /dnd /dn10 5 0 5 10 r ( a )10080604020020 d f ( r ) / d | n | Si Hd /d| n| d /d| n|
10 5 0 5 10 r ( a )10010203040 d f ( r ) / d Si H d /dd /d
Figure 7:
The α and β iso-orbital indicators for SiH showing covalent bonding,plotted along one Si-H internuclear axis. i g u r e : C o n t o u r p l o t s o f α ( r ) ( l e f t) a nd β ( r ) ( r i g h t) f o r t h e h y d r og e nb o nd e d ( H O ) w a t e r d i m e r e v a l u a t e d f o r s e l f c o n s i s t e n t P B E d e n s i t y . A t o m i c p o s i t i o n s a r e l a b e ll e d , t h e l a b e l s f o r o u t o f p l a n e h y d r og e n a t o m s h a v e b ee np r o j ec t e d o n t o t h e p l o tt i n g p l a n e . A s t h e r a n g e o f t h e f un c t i o n s i s d i ff e r e n t , c o l o u r m a pp i n g i s n o t s h a r e db e t w ee np l o t s . A dd i t i o n a ll y , a s α ( r ) i s unb o und e d f r o m a b o v e r e g i o n s f o r w h i c h α ( r ) > . h a v e b ee n c o l o u r e d r e d . L i n e a r s li ce s a l o n g t h e i n t e r - m o l ec u l a r a x i s ( d o tt e db l a c k li n e o n c o n t o u r p l o t) a r e i n c l ud e d ,i n w h i c h t h e s l o w l yv a r y i n g li m i t s o f . nd . f o r α a nd β r e s p ec t i v e l y a r e m a r k e d w i t h a r e dd o t e d li n e . G e o m e tr y i s t a k e n f r o m t h e S t e s t s e t , o p t i m i s e d a tt h e B L Y P l e v e l. erivation of the MS2 β functional The semilocal exchange functional can be written as a construction around an exchangeenhancement factor modulating the exchange energy density of the uniform electron gas, E x [ n ] = Z e UEGx ( n ) F x ( s, α ) d r , (1)where e UEGx ( n ) = − (cid:16) π (cid:17) / n ( r ) / and F x is the exchange enhancement factor.The MS2 functional is designed to recover the second order gradient expansion of ex-change and correlation energies of slowly varying densities. For slowly varying densities,lim |∇ n |→ α = τ UEG + |∇ n | n + ∇ n − |∇ n | n τ UEG (2)= 1 − A + 16 B + O ( ∇ ) , (3)where, A = |∇ n | e UEGx , (4) B = ∇ n/τ UEG , (5)Here, a similar derivation is made for the β quantity,lim |∇ n |→ β = τ UEG + |∇ n | n + ∇ n − |∇ n | n τ UEG + |∇ n | n + ∇ n + τ UEG , (6)= 12 (cid:20) − A + 112 B (cid:21) + O ( ∇ ) . (7)The MS2 exchange enhancement function, F MS2x ( p, α ) is constructed as an interpolation, f ( α ), of two enhancement factor functions of the squared reduced density gradient, p = s = 14(3 / π / ) |∇ n | n / , (8)31uilt for the α = 0 and α = 1 limits as F ( p ) and F ( p ) respectively, F ( p ) = 1 + κ − κ µ GE p + cκ (9) F ( p ) = 1 + κ − κ µ GE pκ , (10)where µ GE = 10 /
81, and κ = 0 .
504 and c = 0 . F and F are independent of α , and hence β , though the κ and c parameters should beadjusted accordingly.The enhancement factor interpolation is a direct function of α , f ( α ) = (1 − α ) α + bα , (11)where b = 4 . α > κ in Ref. by fitting to atomisation energy and barrier height training sets, AE6 andBH6 respectively. Here, we construct a similar interpolation function, exhibiting the samelimiting behaviours, in terms of β . We write,lim |∇ n |→ f ( α ) = lim |∇ n |→ (1 − α ) α + bα (12)= [ − − A + 16 B )] + O ( ∇ ) (13)and this term only contributes to the order of ∇ and higher. We can therefore replace α by 2 β in Eq. 11 without affecting the second order gradient expansion of the MS2 exchangeenhancement factor F MS2x of Ref. We propose, f m ( β ) = [1 − (2 β ) ] β ) + b m (2 β ) , (14)where the m subscript indicates a modified term. Firstly, we confirm that Eq. 14 exhibitsthe same behaviours as Eq. 11 with respect to the equivalent points highlighted in the main32ext, f m ( β = 0) = f ( α = 0) = 1 , (15) f m ( β = 0 .
5) = f ( α = 1) = 0 , (16)the α → ∞ limit corresponds to β → f m ( β = 1) = (1 − b m . (17)The f ( α ) function gives, f m ( β = 1) = lim α →∞ f ( α ) , (18)= − b , (19)(20)hence for f m ( β →
1) to coincide f ( α → ∞ ) then the b m parameter is defined relative to theoriginal b as, b m = 27 b − . (21) Atomic Total Energies
Table 1: The total energies of the neon and argon atoms calculated in the aug-cc-pVQZ basis using a very fine (Turbomole Level 7) integration grid convergedto − E h in the SCF cycle. MS2 MS2 β Ne -128.96625003 -129.02147785Ar -527.53776112 -527.5377611233 agnetic Fields
We note that the required modification for kinetic energy density dependent functionals instrong magnetic fields of substituting the orbital kinetic energy density with the phys-ical kinetic energy density, as derived by Dobson in Ref., indicator quantity is equallyapplicable to the β quantity as α to construct current sensitive meta-GGA functionals. References (1) Jureˇcka, P.; ˇSponer, J.; ˇCern´y, J.; Hobza, P. Benchmark database of accurate (MP2and CCSD(T) complete basis set limit) interaction energies of small model complexes,DNA base pairs, and amino acid pairs.
Physical Chemistry Chemical Physics , ,1985–1993.(2) Becke, A. D. A new mixing of HartreeâĂŞFock and local density-functional theories. The Journal of Chemical Physics , , 1372.(3) Bates, J. E.; Furche, F. Harnessing the meta-generalized gradient approximation fortime-dependent density functional theory. Journal of Chemical Physics , ,164105.(4) Furness, J. W.; Verbeke, J.; Tellgren, E. I.; Stopkowicz, S.; Ekstr¨om, U.; Helgaker, T.;Teale, A. M. Current Density Functional Theory Using Meta-Generalized GradientExchange-Correlation Functionals. Journal of Chemical Theory and Computation , , 4169–4181.(5) Dobson, J. F. Spin-density functionals for the electron correlation energy with automaticfreedom from orbital self-interaction. Journal of Physics: Condensed Matter , ,7877–7890. 346) Dobson, J. F. Alternative expressions for the Fermi hole curvature. The Journal ofChemical Physics ,98