Partially deorbitalized meta-GGA

Mejia-Rodriguez and Trickey recently proposed a procedure for removing the explicit dependence of meta-GGA exchange-correlation energy functionals $E_{\rm xc}$ on the kinetic energy density $\tau$. We present a simple modification to this approach in which the exact Kohn-Sham $\tau$ is used as input for $E_{\rm xc}$ but the functional derivative of $\tau$ with respect to the density $\rho$, required to calculate the potential term $\int d^3r'\,\delta E_{\rm xc}/\delta\tau({\bf r}')|_{\rho}\cdot \delta\tau({\bf r}')/\delta\rho({\bf r})$, is evaluated using an approximate kinetic energy density functional. This ensures that the Kohn-Sham potential is a local multiplicative function as opposed to the non-local potential of a generalized Kohn-Sham approach. Electronic structure codes can be easily modified to use the new method. We validate it by quantifying the accuracy of the predicted lattice parameters, bulk moduli, magnetic moments and cohesive energies of a large set of periodic solids. An unanticipated benefit of this method is to gauge the quality of approximate kinetic energy functionals by checking if the self-consistent solution is indeed at the variational minimum.


I. INTRODUCTION
The exchange-correlation energy functional, E xc , is an essential component of the Kohn-Sham (KS) [1] formulation of density functional theory (DFT) [2].Knowledge of E xc allows for both the exact electronic density and the total energy of a system of interacting electrons in an external potential to be determined.In practice however, E xc has to be approximated, the inevitable result of which is inexact densities and energies.Consequently, much effort has been expended over that past several decades on inventing ever more sophisticated approximations.
The simplest approximation for the unknown exchange and correlation functional is the Local Density Approximation (LDA), in which the functional is supposed to depend locally on the charge density.LDA represents the zeroth order expansion of E xc in terms of electron density gradients and constitutes the first rung of the so-called 'Jacob's ladder' of density functional approximations [3].On the second rung, the Generalized Gradient Approximation (GGA), takes into account density gradients [4], and satisfies more known properties of the exact functional.However, because of the limitations of a semilocal functional form, GGAs tend to be accurate for either energies or equilibrium geometries, but not both [5].The third rung of functionals, meta-GGAs, were conceived to address this limitation by including the Kohn-Sham kinetic energy (KE) density explicitly in their formulation [6][7][8].This class of functionals is truly nonlocal because of the implicit dependence of the KE density on the density itself.This additional flexibility allows for more known exact conditions to be satisfied.In fact, the strongly constrained and appropriately normal-ized (SCAN) functional [9,10] satisfies all 17 known exact conditions that a meta-GGA can.
For the spin unpolarized case, the meta-GGA exchange-correlation energy is given by is the density and is the non-interacting, Kohn-Sham KE density, with ϕ i being the ith Kohn-Sham orbital of a state with N electrons.For its intended use, the KE density is not an independent variable but rather an implicit functional of the density, i.e. τ (r) ≡ τ [ρ](r).Thus a difficulty arises when one has to determine the exchange-correlation potential v xc as the functional derivative of E xc with respect to the density: The last term requires the functional derivative of τ with respect to ρ.This is numerically difficult to perform and requires an approach similar to that used for the optimized effective potential (OEP) method [11].Instead, codes typically calculate potentials determined from the derivative with respect to the orbitals δE xc /δϕ(r).Such an approach, however, produces a non-local potential [12] and is referred to as generalized Kohn-Sham (gKS).This presents us with a dilemma: the great effort expended to satisfy as many exact constraints as possible is undermined by violating a fundamental property of the KS potential, namely that it be a local function in r.And yet, the choice of τ by the inventors of meta-GGA forces the writers of electronic structure codes to have to deal with the difficult functional derivative, or avoid it altogether with gKS.
Mejía-Rodríguez and Trickey (MRT) neatly sidestepped this problem by replacing the τ determined from the orbitals via Eq.( 2) with one obtained from an approximate KE density functional [13][14][15].This 'deorbitalized' meta-GGA was found to produce results of accuracy which were comparable to that of the gKS method.This approach however, removes the true non-locality of E xc and in effect reduces meta-GGA to a semi-local GGA-like functional (albeit possibly with Laplacian terms [16]).
In the current work we adopt a 'half-way' strategy, in that we use the KS orbital-derived τ as input to the functional, but use an approximate KE density functional to evaluate the functional derivative in Eq. ( 3).We term this approach 'partial deorbitalization' and find that even a fairly primitive KE functional, like the Thomas-Fermivon Weizsäcker (TFvW) gradient expansion [17], yields accurate results.This deorbitalization scheme can be easily implemented in codes which already employ the generalized Kohn-Sham version of meta-GGA.Partial deorbitalization retains the 'exact' τ for the energy but utilizes both the exact and approximate τ for the potential, which favors situations in which the main error is functional-driven rather than density-driven [18].

II. APPROXIMATIONS
As will be demonstrated later, the method does not require a particularly sophisticated KE density functional for calculating the functional derivative δτ (r )/δρ(r) in Eq. ( 3).Here we choose to use the gradient expansion of τ with the TFvW terms [17] for the sake of ease of implementation: where σ(r) ≡ |∇ρ(r)| 2 .In this case, the second term in Eq. (3) becomes where The spin polarized case is a straight-forward extension to the unpolarized case: the meta-GGA functional is generalized to the collinear form where ρ ↑↓ and τ ↑↓ are the up-and down-spin density and KE density, respectively.The total Kohn-Sham KE satisfies [19] T thus we will take spin-up KE density to depend exclusively on the spin-up density [16], , and likewise for the spin-down density.This implies that δτ ↑ (r )/δρ ↓ (r The spin-up exchange-correlation potential, for example, is then given by which can be easily evaluated for the spin-polarized TFvW KE density

III. COMPUTATIONAL DETAILS
The goal of the following section is to evaluate the accuracy of our partial deorbitalization strategy.We do so by validating its ab initio predictions against a large number of experimental results.We also list the computational outcomes of the previous work by MRT [14], including both the results from their gKS and the fully deorbitalized meta-GGA simulations, in order to compare against different strategies for deorbitalizing meta-GGA.Notably, the comparison among these approximations, i.e. our work and the findings of MRT, must be performed 'cum grano salis', since the simulations have been carried out using slightly different conditions.MRT [14] used the plane-wave based code VASP [20,21], adopting PAW pseudopotentials [22], and the SCAN exchange and correlation functional [43].In what follows we will call this combination gKS-SCAN when referring to simulations performed with the gKS, while the label FD-SCAN will indicate their fully deorbitalized results.In our work we opt instead for the rSCAN functional [23] in order to overcome convergence problems with full potential (FP) simulations.The partially deorbitalized method introduced above will therefore be labelled PD-rSCAN.The label should immediately remind the reader of the two main ingredients to be considered: the deorbitalization scheme (if any) and the choice of the functional providing the exchange and correlation contribution.Finally, we also performed a number of simulations within the gKS scheme and we refer to these results as gKS-rSCAN, to distinguish them from the gKS results of MRT, labeled gKS-SCAN.
In order to compute equilibrium lattice parameters and bulk moduli with high accuracy, we opt for a FP description of the electronic wave-functions with an APW basis [24].This choice is effective for periodic systems but makes it difficult to compute the total energy of isolated atoms.To overcome this problem, we also use a planewave basis set and pseudopotentials for the estimation of the cohesive energies.
The Elk code [25]  For plane-wave simulations we use the Quantum ESPRESSO package [26] and opt for norm-conserving pseudopotentials [27] generated with the PBE [4] exchange and correlation functional.The reciprocal space sampling and the cut-off energy for the KS wavefunction expansion have been converged in order to obtain better than 1 mRy/atom accuracy in the total energy.
Estimation of atomic energies with meta-GGA requires further care.The exponentially vanishing charge produces problematic behavior in the exchange and correlation potential across the self consistent cycles.There are two options to improve the convergence.The first is to converge isolated atom simulations with GGA (we used PBE) and later reuse the converged electronic charge as the starting point of meta-GGA simulations (for both gKS-rSCAN and PD-rSCAN).The alternative method consists of introducing a cut-off for vanishing charge density that removes the ill-behaving contributions from the exchange and correlation potential.This parameter can be converged together with the remaining settings governing the basis expansion.Numerically equivalent results are obtained with both approaches, when convergence can be achieved.The details are reported in the Supplemental Information.
Finally, we point out that a more recent analysis [15] of deorbitalized meta-GGA adopts r 2 SCAN [28] and the authors report that the convergence issues discussed in their previous work [14] appear to be due to the functional itself rather than to their deorbitalization strategy.We tried r 2 SCAN but still found difficulties converging FP simulations and therefore, in order to preserve consistency among our PW and FP results, we abandoned this option.
Following MRT, the equilibrium lattice constants a 0 and bulk moduli B 0 at T = 0 K were determined by calculating the total energy per unit cell in the range V 0 ± 10% (where V0 is the equilibrium unit cell volume), followed by a twenty point fit to the stabilized jellium equation of state (SJEOS) [14,29].

IV. RESULTS AND DISCUSSION
The complete set of results obtained for the computation of the equilibrium lattice parameters of 51 solids is presented, together with statistical data analysis, in Fig. 1, where panels (a) and (b) report the absolute and relative deviations.All absolute values are reported in Table SIII of the supplemental material (SM) [44].Values for gKS-SCAN and FD-SCAN are taken from Ref. [14] while our results for PD-rSCAN are shown in the third column of the figure and in the fifth column of the table.Finally, the reference experimental data are based on zero-point corrected experimental lattice constants, detailed in Ref. [31] and references therein.
FD-SCAN and PD-rSCAN show similar trends and comparable agreement with respect to experimental lattice parameters.Alkali metals are a notable exception, where deviations by more than 2% [45] are observed.This is indeed a known issue with SCAN that originates from a poor description of the semi-core region of these elements, as explained in detail in Ref. [32].No clear indications of over-binding or under-binding can be identified when considering at the whole set.Yet we note that FD-SCAN and PD-rSCAN show similar trends across the periodic table.The mean absolute deviation (MAD) is 0.028 Å, 0.026 Å and 0.036 Å for gKS-SCAN, FD-SCAN and PD-rSCAN respectively.While the distribution of deviations for PD-rSCAN is slightly broader than the other two, it is noted that outliers are found in the upper part of the box-plot for both gKS-SCAN and PD-rSCAN meta-GGA, while the fully deorbitalized approach of MRT behaves differently, showing only outliers on the opposite side of the other two distributions.
Fig. 2 reports the same statistical analysis for the bulk moduli of various cubic and hexagonal systems.All results are also tabulated in Table SIV of the SM.The most significant discrepancies are observed also in this case for alkali metals.Additional outliers are found in transition metal elements, Cu and Au, and in wide-bandgap semiconductors GaN and BN.The distribution of PD-rSCAN and gKS-SCAN compare equally well against the experiment, while FD-SCAN is showing slightly worse performance, as it can be appreciated from the box-plots of Fig. 2. The MAD for the bulk modulus is 6.8, 9.4 and 6.7 GPa for gKS-SCAN, FD-SCAN and PD-rSCAN respectively.As expected, our approach gives bulk moduli that are similar to those gKS-SCAN, but the partially deorbitalization formulation shows slightly improved values on average with respect to FD-SCAN.
The Kohn-Sham (KS) band gaps for selected insulators and semiconductors are shown in Table I.The results obtained with PD-rSCAN are similar to the ones produced by FD-SCAN, and in both cases the values are smaller than those obtained with gKS-SCAN.This systematic difference is a well known property of the generalized KS theory and an accurate analysis of this point is presented in Ref. [33].The results of FD-SCAN and PD-rSCAN are instead obtained with a local potential and are therefore   II along with experimental atomization energies from Ref. [14,35,36].As already mentioned, these results have been obtained with PW based simulations and in this case a larger discrepancy between the ab initio predictions and the Notably, MRT results obtained with the gKS-SCAN approach are in slightly better agreement than our gKS-rSCAN (first and third columns).This may be due to the norm conserving pseudopotentials used in our simulations that have been generated with a different functional (PBE) for the core electrons and miss the KE density contribution.The comparison between PD-rSCAN and gKS simulations (third and fourth columns) is instead showing perfect agreement, with the only exception being K. Magnetic properties, reported in Table III, are the most sensitive to the choice of the deorbitalization scheme.It has indeed already extensively discussed how gKS with SCAN leads to overhestimated magnetization in transition metal elements [37][38][39].On the other hand, FD-SCAN and PD-rSCAN improve the agreement with experimental results for the elemental ferromagnets Fe, Co and Ni, and also predict the expected non-magnetic ground state for vanadium.The density of states for the four elemental solids obtained with LDA and PD-rSCAN is shown in Fig. 3. Small differences in the densities of states of the ground states can be appreciated.The plot shows that, relative to LDA, PD-rSCAN shifts the spin majority occupied states downward, while the spin minority state energies are only slightly increased in all elements but iron, where the effect is more pronounced but only in the conduction bands.The overall effect is very limited and indeed the resulting magnetic moments are very close.
For the anti-ferromagnetic, insulating magnetic oxides FeO, CoO, NiO and MnO the picture is more mixed.The atomic moment of MnO is close to the experimental value and that of FeO lies within the admittedly broad range of measured moments.However, the moments of CoO and NiO are underestimated by both FD-SCAN and PD-rSCAN.This has been attributed to strong correlation effects which are not fully described even by meta-GGA functionals.
By calculating the total energy while keeping the moment constrained to a given value (referred to as a fixed spin-moment calculation), it can be ascertained if the self-consistent solution is truly variational.The energy vs. moment for Fe, Co and Ni is plotted in Fig. 4 along with the moment from the corresponding self-consistent solutions.As can be seen, the moments are generally smaller than the location of the minima of the curves.This implies that the calculations are not perfectly variational, which in turn implies, unsurprisingly, that the approximate KE functional is inexact.Fully deorbitalized functionals will not suffer from this inconsistency because the functional derivative is determined from same KE functional as is used to evaluate the energy (inexact though it is).Nevertheless, this mismatch does serendipitously present us with a useful tool for determining the accuracy of approximate KE functionals: simply check if the solution is at the energy minimum.This is particularly useful because these calculations are for real-world atomistic systems as opposed to simplified models.Magnetic moments can be used as done here, but there are other possibilities: for example testing if the Helmann-Feynman forces on the nuclei are strictly zero at the lowest energy.

V. CONCLUSIONS
We have modified the deorbitalized approach introduced by MRT by retaining the true Kohn-Sham KE density as input to any meta-GGA exchange-correlation functional, while employing an approximate KE functional for the functional derivative of τ with respect to ρ.This represents a simple route to keeping both E xc inherently non-local and v xc as a local operator in r.DFT codes which already have the gKS version of meta-GGA implemented can be easily modified to accommodate this scheme.We find that even a relatively crude KE functional like TFvW yields results which are, on the whole, at least as accurate as competing approaches such as generalized Kohn-Sham or full deorbitaliztion.Furthermore, there is ample scope for improvement with respect to both E xc as well as the KE functional.For example, KE functionals involving the Laplacian of the density [16,19] could further enhance the accuracy of method.Showing that the self-consistent solution does indeed correspond to the energy minimum for solids and molecules is a useful measure of the quality of the KE functionals in realworld situations.Lastly, our method may be extended to the case of non-collinear exchange-correlation meta-GGA functionals, at least two of which have been developed recently [40,41].Treating this type of functional with partial deorbitalization will also require a generalization of the spin-dependent KE density functional to the non-collinear case [42].

Convergence of plane wave based meta-GGA simulations
The estimation of atomic energies with meta-GGA is not straightforward and requires the adoption of workarounds to the pathological behaviour of the rSCAN functional with vanishing charges.We found two options: 1. converge the isolated atom states with PBE and use the resulting density as an input for the meta-GGA simulation, or 2. introduce a cutoff such that We then converge the total energy against the box size, the plane wave cutoff E c , the charge density expansion (which is larger than 4E c in order to achieve numerical stability of meta-GGA) and, when using v cut xc , against ρ c .An example of convergence is shown in Fig. S1.When both method 1. and 2. achieve convergence, the results are found to be numerically equivalent.
The cutoff for charge expansion is set to 18 times E c and a total energy convergence better than 1 mRy is achieved with the following parameters for the atoms listed The estimation of the equilibrium volume of the above elements arranged in periodic solids is not problematic but we stress that unexpectedly large thresholds for the charge density expansion are required to converge the total energy of rSCAN simulations and obtain smooth E vs V curves.The values used for the elemental solids listed above is reported in Tab.S1.Monkhorst-Pack grids were used to sample reciprocal space and a Methfessel-Paxton smearing of 0.01 Ry was adopted for metals.

Tabulated results
The results obtained for all materials of our testing set are reported in Tab.S2 and Tab.S3 and compared with the outcomes of the previous investigation by Mejia-Rodriguez and Trickey [14].

FIG. 1 :FIG. 2 :
FIG.1: Visualization (violin plot and box plot) of the deviations from experimental data of the predicted lattice parameter obtained with the different implementations of meta-GGA.The upper(lower) panel gives absolute(relative) deviations.The violin plots (transparent color) represent the data distribution and are based on a Gaussian kernel density estimation implemented in seaborn[30].In the box plot, the boxes hold 50% of the data, with equal number of data points above and below the median deviation (full black line).The whiskers indicate the range of data falling within 1.5×boxlength beyond the upper and lower limits of the box (from the first quartile to the third quartile).The whiskers extend from the box by 1.5× the inter-quartile range.Outliers beyond this range are indicated with circular makers and the solids' labels are reported on the right of each point.

FIG. 3 :
FIG. 3: Spin projected density of states (DOS) in states/spin/eV for (a) Fe (b) Co (c) Ni and (d) V. Results are shown for LDA as well as PD-rSCAN.
FIG. S1: Convergence of total energy of isolated Ba atom against charge threshold ρc.The red lines highlight 1 mRy accuracy.The yellow point shows the result obtained by using converged charge density obtained with PBE as the starting point for the meta-GGA simulation.
version 8.3.15 is used to perform FP simulations.Reciprocal space sampling is performed with at least a 17 × 17 × 17 Monkhorst-Pack grid.The basis is expanded up to 2|G + k| max ≥ 8.The remaining parameters are set by the vhighq option.

TABLE I :
[34]ulated Kohn-Sham band gaps in eV for 20 insulators or semiconductors in the test set.The last column reports optimized effective potential results from Yang et al.Ref.[33]obtained with the Krieger-Li-Iafrate approximation[34].

TABLE III :
[37]ulated magnetic moment (in µB) for a selection of ferromagnets and anti-ferromagnets.Experimental results are taken from Ref.[37]and references therein.

TABLE S1 :
Parameters used to compute E vs V curves with the plane wave basis.Energies are in Ry.