Improving the applicability of the Pauli kinetic energy density based semilocal functional for solids

The Pauli kinetic energy enhancement factor $\alpha=(\tau-\tau^W)/\tau^{unif}$ is an important density ingredient, used to construct many meta-generalized gradient approximations (meta-GGA) exchange-correlation (XC) energy functionals, including the very successful strongly constrained and appropriately normed (SCAN) semilocal functional. Another meta-GGA functional, known as MGGAC [Phys. Rev. B 100, 155140 (2019)], is also proposed in recent time depending only on the $\alpha$ ingredient and based on the generalization of the Becke-Roussel approach with the cuspless hydrogen exchange hole density. The MGGAC functional is proved to be a very useful and competitive meta-GGA semilocal functional for electronic structure properties of solids and molecules. Based on the successful implication of the ingredient $\alpha$, which is also useful to construct the one-electron self-interaction free correlation energy functional, here we propose revised correlation energy for MGGAC exchange functional which is more accurate and robust, especially for the high and low-density limits of the uniform density scaling. The present XC functional, named as revised MGGAC (rMGGAC), shows an impressive improvement for the structural and energetic properties of solids compared to its previous version. Moreover, the assessment of the present constructed functional shows to be quite useful in solid-state physics in terms of addressing several current challenging solid-state problems.


I. INTRODUCTION
The Kohn-Sham (KS) density functional theory (DFT) 1 has become an indispensable computational framework for performing the electronic structure calculations of new generation electronic devices 2-4 , energy storage devices 5 , catalyst 6 and spintronics devices 7-9 . Despite the extraordinary potential and considerable progress made to date, challenges remain to bridge the gap between simulations and experimental findings. Efforts have increased in developing new strategies in understanding the theory and practical implementation of complex systems, but current progress is still far from sufficient. However, recent strategies of the development of non-empirical exchange-correlation (XC) functionals by satisfying quantum mechanical exact constraints 10 exhibit extremely prosperous to catch the interplay of dimensionality, correlation, charge, orbital character, topological, and other large reservoirs of exotic properties of condensed matter systems 11 .
The α ingredient helps to impose some exact constraints in the functional form, such as the strongly tightened bound of exchange, i.e. F x ≤ 1.174 10,122,123 . Importantly, this allows a negative slope (∂F x /∂α ≤ 0) which is related to an improvement in the prediction of the band gap 84 . Also, the use of solely α and s = |∇ρ|/[2(3π 2 ) 1/3 ρ 4/3 ] (reduced density gradient) makes the meta-GGA XC functional free from the order-of-limit anomaly, that is important for the structural phase transition pressures and energies 73,124 . Focusing on the XC functional development, a popular non-empirical meta-GGA XC functional is the Strongly Constrained and Appropriately Normed (SCAN) 10 semilocal DFA. Besides the SCAN, it was also recently proposed the MGGAC XC functional 80 , whose exchange part depends only on the α ingredient, being developed from a generalization of the Becke-Roussel approach 78,125 , and using the cuspless hydrogen model for the exchange hole density. The MGGAC functional showed its productive power over other meta-GGAs in terms of the band gap performance for bulk solids 126 , two dimensional (2D) van der Waals (vdW) materials 2,3,126 , and structural phase stability of solids 127 .
The MGGAC exchange enhancement factor is a simple monotonic function of α and satisfies (i) the strongly tightened bound of exchange F x ≤ 1.174 (F x = 1.174 for two-electron systems), (ii) the uniform electron gas limit (F x = 1) when α = 1, and (iii) the cuspless hydrogen related behavior (F x = 0.937) at α → ∞. On the other hand, the MGGAC correlation energy functional is based on the GGA PBE correlation form 128 , where the linear response parameter (β(r s ) = 0.030) was fitted to the equilibrium lattice constants of several bulk solids. Then, it is not one-electron self-interaction free, an important constraint that should be satisfied by every meta-GGA functional. Nevertheless, the MGGAC DFA has demonstrated several successes which include (i) improved thermochemical properties 80 , (ii) improved fundamental band gap for bulk and layer solids 126 , and (iii) correct energy ordering for challenging polymorphs of solids 127 . However, the functional has its own limitations such as moderate performance for lattice constants, bulk moduli, and cohesive energies for bulk solids 80 . Therefore, it is desirable to improve the accessibility of the MGGAC functional for solids by removing some deficiencies of the current version. In order to do so, we keep the MGGAC exchange unchanged, modifying only the correlation part of DFA energy such that it becomes exact for one electron densities and more robust for low-and high-density limits.
The paper is organized as follows: in the next section, we will construct the revised MGGAC correlation energy functional dependent on both α and s. Later, we will couple the proposed correlation with the MGGAC exchange to assess the functional performance for molecular and solid-state properties. Furthermore, we also investigate the low-and high-density limit of the XC functional. Lastly, we summarize our findings.
The parameters γ 1 = 0.08 and γ 2 = 0.3 are chosen from the atomization energies of the AE6 molecules 130 . Note that in the tail of the density, if the valence atomic orbital is degenerate, as in the case of Ne atom, then α → ∞, else α → 0 as in case of alkali and alkaline-earth atoms.
In Fig. 1 we show a comparison between rMGGAC, MGGAC and SCAN correlation energies per particle, for Ne atom. By construction, rMGGAC recovers ǫ 0 c at the nucleus (where α ≈ 0) and in the tail of the density (where s diverges), while it recovers ǫ 1 c in the atomic core (where the density is compact and slowly-varying). Overall, rMGGAC agrees closely with, but is smoother than SCAN, being significantly different from the MGGAC behavior.
We mention that the choice of the f 1 (α, s), f 2 (α, s), and ǫ rMGGAC c keeps the correct constraint of the correlation energy functional. Such as, in the rapidly varying density limit, recognized as t → ∞ and s → ∞, both H 0 and H 1 correctly cancel with ǫ LSDA0 c and ǫ LSDA1 c , respectively, making the correlation vanishes, an important constraint satisfied by all functionals constructed on the top of PBE form 54 . Nevertheless, as shown in Fig. 1, this vanishing process is much slower for rMGGAC and SCAN than for MGGAC. Here we also note that in the tail of the density, if the valence atomic orbital is degenerate, as in the case of Ne atom, then α → ∞, else α → 0 as in case of alkali and alkaline-earth atoms.  In the low-density limit (see APPENDIX A for details), the rMGGAC XC energy functional becomes independent of ζ for 0 ≤ ζ ≤ 0.7 (as shown in Fig. 2 for one electron Gaussian density). The ζ independence of the XC functional is an important constraint not only for Wigner crystals, but also for the atomization energies 131 .
To quantify the performance of correlation functional for high-density limit, in Table I, we consider the correlation energies for the benchmark test of several atoms and ions, also used in Ref. 133 . It is shown that the rMGGAC correlation energy considerably improves over MGGAC, being in line with the SCAN correlation.
Finally in Fig. 3 we compare the SCAN, MGGAC, and rMGGAC XC enhancement factors F xc (ζ, α, s, r s ), for the spin-unpolarized choice (ζ = 0) and α = 0 (that is describing the two-electron singlet states) and α = 1 (the case of slowly-varying density limit when s ≤ 1). In both panels, the r s = 0 represents the exchange-only (high-density) limit, and we observe that MGGAC (rMGGAC) enhancement factor is independent on s (F MGGAC x = 1.174 for α = 0, and F MGGAC x = 1 for α = 1). Overall, the rMGGAC XC enhancement factor is slightly bigger than SCAN for α = 0, and slightly smaller than SCAN α = 1. We also note the large differences between MGGAC and rMGGAC, especially for the low-density cases (r s = 10).

III. RESULTS AND DISCUSSIONS
In this section, we present the assessment of the performance of the rMGGAC XC DFA, along with SCAN and  Solids for lattice constants, bulk moduli, and cohesive energies MGGAC methods for harmonium atom, various molecular and solid-state problems including the band gap assessment and structural phase transition of challenging systems.

A. Model system: Harmonium atom
We have tested the performance of the several DFAs (including the rMGGAC) against the harmonium atom 135 , for various values of the confinement strength ω. We recall, that at small values of ω, the system is strongly correlated, whereas for large values of ω it is tightly bounded. These two regimes are very important for many condensed matter applications. Hence, the harmonium atom provides an excellent tool for testing ap-proximate density functional methods 136,137 . Similarly to our past studies 138,139 , in order to perform calculations we have utilized an even-tempered Gaussian basis set from Ref. 140 for ω ∈ [0.03, 1000]. As a reference results, we have used full configuration interaction (FCI) data obtained in the same basis set which have been proved to be close to exact values 138 .
In Fig. 4, we report the relative error (RE) on total energies calculated with respect to reference data for several XC functionals as a function of ω. All DFAs energies have been obtained in post-SCF fashion on top of self-interaction free OEPx (optimized effective potential exact exchange) densities as in Refs. 138 and 139. One can note the remarkable good performance of rMGGAC functional along with the whole range of ω with an overall MARE of 0.21%. Moreover, rMGGAC significantly outperforms the MGGAC (MARE=2.67%), being especially visible for the strongly correlated regime. This may indicate that the new correlation functional form is more compatible with MGGAC exchange. The rMGGAC behavior is also quite similar to BLOC functional 141 which also underestimate slightly the reference results with an overall MARE of 0.23%. The SCAN DFA gives similar results (MARE=0.20%), although in this case, we observe the overestimation. For comparison, we also show the TPSS and PBE data which gives an overall MARE of 0.78% and 0.28%, respectively.  134 and additionally we include the band gaps of 9 solids (Cu2O, CuBr, ScN, SnTe, MgO, NaCl, LiCl, NaF, and LiF). The MARE of the SBG31 test is also reported. For Surface energies (ǫ surf ace (in J/m 2 )) we consider the (111) surface. We also calculate the CO molecule adsorption energies (ǫ ads ) on top of the (111) surface (in eV). For both cases we consider five different transition metals taken from Ref. 80 . The SCAN values are taken from Ref. 109 . The surface energies and adsorption energies of SCAN and MGGAC functionals are taken from ref. 80 . The LC20, BM20, COH20, and SBG31 test sets results of MG-GAC are taken from ref. 80 . See ref. 132  For general assessment, we first consider the 41 bulk solids set proposed in Table II and we compute the equilibrium lattice constants (LC41), bulk moduli (BM41), and cohesive energies (COH41). This benchmark test set was also used to assess the performance of SCAN 10 and MGGAC 80 functionals. We report the error statistics (MAE and MARE) for all the solids in TABLE III. Our benchmark calculations for the LC41 test indicate that the rMGGAC functional gives the MAE ≈ 38 mÅ, which is a remarkable improvement performance compared to MGGAC, being close to the state-of-art accuracy of the SCAN meta-GGA. We also report the MAE of the LC20 test set compiled in ref. 142 for which MAEs for other methods (not included here) are also available.  From literature we found other meta-GGAs, like the revised Tao-Perdew-Staroverov-Scuseria (revTPSS) 69 and Tao-Mo (TM) 44 semilocal functionals, that give for LC20 test, the MAE ≈ 32 mÅ 77,142 .
The improvement of the lattice constants from rMG-GAC functional is also followed by its performance for the bulk moduli, showing that the rMGGAC predicts the realistic shape of the energy versus volume curve near the equilibrium point. Thus, we obtain for the BM41 test, an improved MAE ≈ 8.5 GPa from rMGGAC compared to the MGGAC (MAE ≈ 11.5). Our benchmark values are also close to that of the SCAN where the later gives MAE of 7.5 GPa. For BM20 test set we obtain MAE of 4.6 GPa from rMGGAC. The performance of the rMG-GAC for BM20 is also better or comparable than the state-of-art GGA for solids, PBEsol (MAE 6.2 GPa) 142 and meta-GGA revTPSS (MAE 8.7 GPa) 142 and TM (MAE 4.2 GPa) 77 .
The results for the COH41 test also display the improved performance from rMGGAC (MAE of 0.30 eV/atom) over the MGGAC (MAE of 0.36 eV/atom). Here, the SCAN (MAE of 0.19 eV/atom) performs better than rMGGAC. However, our reported MAE values for COH20 indicates that rMGGAC values are close to the GGA PBEsol (MAE 0.25 eV/atom) 142 and TM meta-GGA semilocal functionals (MAE 0.28 eV/atom) 77 . Note that for cohesive energies, well balanced performances for solids and atoms are required.
The different functional performances for the general purpose bulk solid properties are also shown in Fig. 5, where we observe an almost systematic improvement of rMGGAC over the MGGAC. For lattice constants, the Cs is challenging, for which we observe improvement from MGGAC to rMGGAC and to SCAN. Note that being a soft matter (i.e. small bulk moduli systems), the accurate description of the lattice constants of Cs needs a vdW interaction. Here, all methods overestimate the lattice constants indicating the lack of vdW interaction in the functional form. For cohesive energies and bulk moduli, we observe the magnetic Fe and Ni systems are bit of challenging for rMGGAC. Note that all these functionals overestimate the magnetic moments of the itinerant electron ferromagnets 109,111,[143][144][145] which is probably the source of error for bulk moduli and cohesive energies of Fe and Ni. Here, we observe that SCAN is bit better than MGGAC and rMGGAC for Ni cohesive energy.
Next, we consider the band gaps of 40 bulk semiconductors which is constructed from the 31 semiconducting band gaps of SBG31 test taken from Ref. 134 , additionally we include the band gaps of 9 solids: Cu 2 O, CuBr, ScN, SnTe, MgO, NaCl, LiCl, NaF, and LiF. All band gaps are calculated at the experimental lattice constants, where the lattice parameters are taken from ref. 146 . The rMG-GAC performance on band gaps of bulk solids is quite close to that of the MGGAC. Our test set consists of narrow to wide band gap insulators.
Interestingly, we observe that MGGAC and rMGGAC give non-zero band gap even for InSb, a difficult case for SCAN, where the generalized KS (gKS) gap of SCAN is zero. We obtain 1.34 eV gap for InSb from rMG-GAC which is larger than the experimental reported value (0.24 eV). The readers are suggested to go through refs. [147][148][149] for other different functionals performance in bandgap assessments. Nevertheless, the improvement of the MGGAC and rMGGAC band gaps over SCAN method for band gap of solids within 1.5 < E g < 6 eV are clearly visualized from Fig. 6. The improvement in the performance of MGGAC and rMGGAC is due to the use of only α ingredient in the exchange functional form, which gives more negative slope ∂F x /∂α, achieving a stronger ultra-non-locality effect 84 . However, we mention that both MGGAC and rMGGAC overestimate the band gap of low-gap semiconductors, having MARE worse than SCAN. Finally, we consider the surface energies and on top CO surface adsorption energies for (111) surfaces of transition metals Cu, Pd, Pt, Rh, and Ir. These tests are particularly interesting because they involve transition metal surface phenomena, that are important for catalytic applications. We observe that without including any van-der Waals (vdW) long-range corrections, MG-GAC and rMGGAC perform slightly better than SCAN. We recall that SCAN may perform better with a vdW non-local correction 96,150 .

C. Electronic band gap of 2D vdW Transition Metal Dichalcogenide and Monolayers
In a generalized way, 2D vdW materials offer great flexibility in terms of tuning their electronic properties 170,171 . Thus, electronic band-gap engineering can be car- ried out by changing the number of 2D layers in a given material [172][173][174] . The stacking of two or more monolayers of different vdW materials leads to the rich variety of 2D vdW quantum systems and allows researchers to explore novel and collective electronic, topological, and magnetic phenomena at the interfaces 175 . The correlation between electron and spin, phonons, and other electrons can be significantly modulated affecting the charge transport, entropy, and total energy in the quantum limit. Especially, the bandgap engineering of 2D vdW systems from a low cost semilocal methods is still very difficult to achieve where various methods based on the hybrid density functionals and Green functional (GW) level theory is proposed. Also, the GGA functionals perform moderately to estimate the bandgaps of these systems, whereas meta-GGA methods are not vastly assessed. In the present context, we show a systematically improvable performance for bandgaps of those systems from SCAN to MGGAC and/or rMGGAC functionals. To depict this we consider bandgaps of several bulk 2D vdW transition metal dichalcogenide (TMDC) and monolayers, which are presented in Table IV. These systems have electronic band gap ranging from the semiconducting monolayer 2H-MoS 2 (1.8 eV) and black phosphorene (BP) (1.2 eV) to their bulk systems 2H-MoS 2 (1.29 eV) to wide gap insulator of hexagonal boron-nitride (h-BN) (5.4 eV). The results of the aforementioned systems are presented in  TABLE IV. Inspecting TABLE IV we observe, MGGAC and rMGGAC can be considered as improved methods to address the bandgaps of the bulk TMDC to monolayer. In all the cases, the MGGAC and rMGGAC perform in a quite similar way, predicting the bandgaps better than the SCAN functional.

D. Transition pressure and structural stability
Accurate prediction of the phase transition pressure and correct ground state phase of polymorphs of solids are quite important from the point of view of the stability of the solid. Importantly, meta-GGA functionals that suffer from the order-of-limit problem, wrongly predict the phase transition pressure and structural phase transition of different polymorphs of phases 73,177 . But functionals like SCAN and MGGAC methods do not suffer from the order-of-limit problem. In Table V, we report the phase transition pressure of different phases of the solids. We notice that the MGGAC has the tendency of overestimation the phase transition pressure for Si and Ge, while rMGGAC corrects this limitation of the MG-GAC. This is because the energy differences of different phases are better described by rMGGAC than MGGAC functional. Overall, SCAN functional performs better here for all transition pressures.
The improved performance of rMGGAC is also noticeable from the structural phase stability of different polymorphs of FeS 2 , TiO 2 , and MnO 2 . These systems are known to be challenging systems for semilocal functionals as those often fail to predict the correct energy ordering of these systems 127,176,178,179 . While the SCAN functional correctly predicts all phases of MnO 2 polymorphs, it wrongly estimates the ground states of both FeS 2 and TiO 2 . The SCAN functional predicts marcasite (M) and anatase (A) as the most stable phases over the pyrite (P) and rutile (R) phases for FeS 2 and TiO 2 , respectively 127,178,179 . However, it is shown 127 , correct energy ordering of these systems can be achieved from the MG-GAC level theory. Moreover, we observe from Fig. 7 that rMGGAC further improves the performance of MGGAC in terms of the energy differences. For MnO 2 the energy differences as obtained from rMGGAC are very close or even better than those obtained using SCAN functional.

E. General thermochemical assessment
To rationalize the rMGGAC functional performance for small molecules, we consider several standard bench-  shown in Fig. 8. From Ref. 132 and Fig. 8 it is evident that all functionals correctly predict the prism isomer as the most stable water hexamer, whereas the rMGGAC significantly improves over MGGAC for the relative energies (with respect to the prism isomer) of different isomers.

IV. SUMMARY AND CONCLUDING REMARKS
We propose a revised MGGAC (rMGGAC) XC functional by combining a one-electron self-interaction free correlation functional with the original MGGAC exchange. The rMGGAC correlation functional has been constructed to satisfy important exact conditions, being exact for one-electron systems, accurate for two-electron ground-state systems, and the low-and high-density limits of the uniformly scaled density. The rMGGAC correlation energy density is smooth and without spurious structure, as shown in Fig. 1. Moreover, the rMGGAC correlation functional corrects the huge underestimation of the correlation energy of MGGAC for atoms and ions (see Table I).
The rMGGAC XC functional, as its predecessor, is an improved functional for band gap of strongly-bound solids, vdW solids, and 2D TMDC layers, but also improves considerably over MGGAC in several prospects. The performance of the rMGGAC demonstrates its successes in the broad direction of solid-state and quantum chemistry properties, which includes thermochemical, energy ordering for challenging polymorphs of solids and isomers of water hexamers, and structural properties of solids.
Lastly, we conclude that the newly designed rMGGAC functional can be further implemented to design novel materials having advance functionalities in the field of unique device architecture. Also, the rMGGAC can offer opportunities for interesting applications of advance energy storage, catalytic problems, electronics and valleytronics.

COMPUTATIONAL DETAILS
We perform molecular calculations using a development version of Q-CHEM 193 software package with def2-QZVP basis set. The XC integrals are performed with 99 points radial grid and 590 points angular Lebedev grid. Note that the present choice of the grid is adequate for the complete energy convergence of the non-bonded systems.
Solid-state calculations are performed in Vienna Ab initio Simulation Package (VASP) [194][195][196][197] . For all the simulations, we use 20 × 20 × 20 k points sampling with energy cutoff 800 eV. For cohesive energies we consider an orthorhombic box of size 23×24×25Å 3 is considered. The 3 rd order Birch-Murnaghan equation of state is used to calculate the bulk moduli. The surface energy and CO adsorption energies are calculated with 16 × 16 × 1 Γ-centered k points with an energy cutoff of 700 eV and a vacuum level of 20Å . For CO adsorption, we relax the top two layers.

SUPPORTING INFORMATION
The supporting information is attached with the manuscript and available online. Additional details are available to the corresponding author upon reasonable request. ACKNOWLEDGMENTS S.J. is grateful to the NISER for partial financial support. S.K.B acknowledges NISER for financial support. Calculations are performed in the KALINGA and NIS-ERDFT High Performance Computing Facility, NISER. Part of the simulations and/or computations are also supported by the SAMKHYA: High Performance Computing Facility provided by Institute of Physics (IOP), Bhubaneswar. S.J. and P.S. grateful to Prof. Suresh Kumar Patra for the computational facility at IOP, Bhubaneswar. S.J. and P.S. would like to thank Q-Chem, Inc. and developers for providing the source code. S.Ś. is grateful to the National Science Centre, Poland for the financial support under Grant No. 2020/37/B/ST4/02713 with χ = (β 0 /γ)c 2 e −ω/γ ≈ 0.72161, β 0 = 0.066725, γ = 0.031091, c = 1.2277, and ω = 0.046644. See ref. 54