Functional designed to include surface effects in self-consistent density functional theory

We design a density-functional-theory DFT exchange-correlation functional that enables an accurate treatment of systems with electronic surfaces. Surface-specific approximations for both exchange and correlation energies are developed. A subsystem functional approach is then used: an interpolation index combines the surface functional with a functional for interior regions. When the local density approximation is used in the interior, the result is a straightforward functional for use in self-consistent DFT. The functional is validated for two metals Al, Pt and one semiconductor Si by calculations of i established bulk properties lattice constants and bulk moduli and ii a property where surface effects exist the vacancy formation energy . Good and coherent results indicate that this functional may serve well as a universal first choice for solid-state systems and that yet improved functionals can be constructed by this approach.


I. INTRODUCTION
Kohn-Sham ͑KS͒ density functional theory 1 ͑DFT͒ is a method for electronic structure calculations of unparalleled versatility throughout physics, chemistry, and biology.In principle, it accounts for all many-body effects of the Schrödinger equation, limited in practice only by the approximation to the universal exchange-correlation ͑XC͒ functional.In this paper we present an improved XC functional, created with a methodology entirely from first principles, that incorporates a sophisticated treatment of electronic surfaces-i.e., strongly inhomogeneous electron densities.3][4] The result is a systematic improvement of bulk properties of solid state systems and a qualitative improvement for systems with strong surface effects.
The XC functional suggested in the early works on the theoretical foundation of DFT, 1 the local density approximation ͑LDA͒, was derived from the properties of a uniform electron gas, but has shown surprisingly wide applicability for real systems.For solid-state calculations the LDA is still often the method of choice.The next level in functional development, the generalized gradient approximations ͑GGA's͒, in many cases significantly improves upon the LDA.The GGA functionals popular for solid-state applications 5,6 are constructed to fulfill constraints that have been derived for the true XC functional.However, the resulting functionals improve results in an inconsistent way ͑see, e.g., Ref. 4͒.Even worse, these functionals often are less accurate than the LDA for properties involving strong surfaces effects, such as the generalized surfaces of metal monovacancies.Recent work has explained this as a systematic underestimation of the surface-intrinsic energy contribution that, for simple surface geometries, can be estimated by a posteriori procedure. 2,3A recently developed meta-GGA functional by Tao, Perdew, Staroverov, and Scuseria 7 ͑TPSS͒ is able to fulfill yet more constraints of the exact XC func-tional by allowing for a more complicated electron density dependence ͑i.e., through the kinetic energy density of the KS quasiparticle wave functions͒ than the present work does.However, it appears that TPSS does not fully rectify the surface energy problems found for the GGA's.We repeated the post-correction scheme in Ref. 3 for TPSS, using published TPSS jellium XC surface energies, 7 and from this a remaining surface error is predicted.
The present work follows an alternate route to functional development from the traditional path described above.The LDA's use of the uniform electron gas model system leads to physically consistent approximations ͑e.g., compatible exchange and correlation that compose the XC functional͒.Our subsystem functional approach, 8 aims to preserve this propitious property of the LDA through the use of region-specific functionals derived from other model systems.A first effort in this direction was made with the local airy gas 9 ͑LAG͒.It extends the LDA by an exchange surface treatment derived from the edge electron gas model system, 10 but keeps the LDA correlation.This first step is completed with the optimized, compatible, correlation introduced here.It is in this sequence of functional development, the LDA, LAG, and then our functional, that the contribution of the present work is most clear.
The XC energy functional E xc ͓n͔ operates on the groundstate electron density n͑r͒.It is usually decomposed into the XC energy per particle ⑀ xc , E xc ͓n͑r͔͒ = ͵ n͑r͒⑀ xc ͑r;͓n͔͒dr.

͑1͒
Exchange and correlation parts are treated separately, with ⑀ xc = ⑀ x + ⑀ c .We put special emphasis on the conventional, local, inverse radius of the exchange hole 10 definition of the exchange energy per particle, ⑀ ˆx.This is in contrast to expressions based on transformations of Eq. ͑1͒ that arbitrarily delocalize ⑀ x and therefore cannot directly be combined with each other within the same system. 8The LDA is local in this sense, while common GGA functionals 5,6 are not.The LDA exchange term is, in Rydberg atomic units, ⑀ ˆx LDA "n͑r͒… = − 3/͑2͓͒3 2 n͑r͔͒ 1/3 .͑2͒

II. FUNCTIONAL CONSTRUCTION
Kohn and Mattsson 10 put forward the Airy electron gas as a suitable model for electronic surfaces.The Airy gas is a model of electrons in a linear potential, v eff ͑r͒ = Lz.L sets an overall length scale and ⑀ ˆx and n͑r͒ can be rescaled by ⑀ ˆx,0 Airy = L −1/3 ⑀ ˆx͑r ; ͓n͔͒ and n 0 = L −1 n͑r͒.Parametrizations are constructed from the exact ⑀ ˆx,0 Airy and n 0 expressed 11 in Airy functions Ai,

͑6͒
The LAG functional of Vitos et al. 9 uses the ⑀ c of the Perdew-Wang ͑PW͒ LDA ͑Ref.12͒ combined with ⑀ ˆx = ⑀ ˆx LDA F x LAG from an Airy gas corresponding to a generic system's density n͑r͒ and scaled gradient s = ͉ ٌ n͑r͉͒ / ͓2͑3 2 ͒ 1/3 n 4/3 ͑r͔͒.The refinement factor is where a ␣ = 2.626 712, a ␤ = 0.041 106, a ␥ = 0.092 070, and a ␦ = 0.657 946.F x depends only on s since n͑r͒ just sets a global scale of the model via L.However, far outside the electronic surface, F x LAG does not reproduce the right limiting behavior.We have derived an improved parametrization by using ͑i͒ the leading behavior of the exchange energy far outside the surface, 10 ⑀ ˆx,0 Airy → −1 / ͑2͒, ͑ii͒ asymptotic expansions of the Airy functions in Eqs.͑4͒ and ͑5͒, and ͑iii͒ an interpolation that ensures the expression approaches the LDA appropriately in the slowly varying limit, , n ˜0͑s͒ = ˜͑s͒ 3/2 3 2 s 3 , ͑8͒ using a superscript LAA for the local Airy approximation, the Lambert W function, 13 and where c = 0.7168 is from a least-squares fit to the true Airy gas exchange.Figure 1 shows that the improvement of the LAA over the LAG is small in the intermediate region, but pronounced outside the surface.
The Airy exchange parametrizations are designed to accurately model the electron gas at a surface.Hence, they cannot be assumed to successfully work for interior regions.The subsystem functional approach 8 uses an interpolation index for the purpose of categorizing parts of the system as surface or interior regions.We use a simple expression where ␣ is determined below.
In the present work the LDA is used in the interior.In the limit of low s, the LAG and LAA already approach the LDA exchange.The end result for the interpolated exchange functional is therefore only slightly different from using the LAG or LAA in the whole system.However, interpolation is needed for the correlation and to enable future use of other interior exchange functionals.
No "exact" correlation has been worked out for electrons in a linear potential.To obtain a correlation functional, we combine the LAA or LAG exchange with a correlation based on the LDA, but with a multiplicative factor ␥. The numerical value of ␥ is given by a fit to jellium surface energies xc .For a functional ⑀ xc ͑r ; ͓n͔͒, xc = ͐n͑z͓͒⑀ xc ͑r ; ͓n͔͒ − ⑀ xc LDA ͑n ¯͔͒dz, where n͑r͒ is from a self-consistent LDA calculation on a system with uniform background of positive charge n ¯for z ഛ 0 and 0 for z Ͼ 0 ͑Ref.14͒.The value of n ¯is commonly expressed in terms of r s = ͓3/͑4n ¯͔͒ 1/3 .The most accurate XC jellium surface energies are given by the improved random-phase approximation scheme presented by Yan et al. 15 RPA+.We minimize a least-squares sum ͚ r s ͉ xc approx − xc RPA+ ͉ 2 , using values for r s = 2.0, 2.07, 2.3, 2.66, 3.0, 3.28, and 4.0.The surface placement ␣ and the LDA correlation factor ␥ are fitted simultaneously 16 : ␣ LAA = 2.804, ␥ LAA = 0.8098.

͑11͒
The resulting fit reproduces the jellium XC surface energies with a mean absolute relative error ͑MARE͒ less than half a percent; cf.Fig. 2 and Table I.
The final form of the functional is where F x ͑s͒ is either from Eq. ͑7͒ or from Eq. ͑8͒, and ⑀ c LDA is the PW LDA correlation. 12

III. TESTS
Numerical tests were performed with the plane-wave code SOCORRO.17,18 Pseudopotentials ͑PP's͒ were generated with the FHI98PP code, 19 modified to obtain the XC potential from a numerical functional derivative.We use settings provided by the included element library. 18The PP's and code modi-fications have been extensively tested.In addition to the functionals presented by this paper, PP's were generated for the LDA, the GGA of Perdew and Wang ͑PW91͒ 5 , and the GGA of Perdew, Becke, and Ernzerhof ͑PBE͒ 6 .For the latter, bulk calculations with PP's constructed with our numerical functional derivatives agree with the results of PPs based on analytical functional derivatives within 0.001%. 6We also obtain reasonable agreement with the all-electron bulk results in Ref. 4. As the tools for PP analysis could not easily be made to use numerical derivatives, an analysis was done for PP's of the above functionals with analytical derivatives using identical settings.These PP's were found to have satisfactory logarithmic derivatives and pass the built-in ghoststate tests. 18he tests presented here have been chosen from a condensed-matter point of view: three elements for which the LDA and PBE give similar as well as different results.
The tests include materials where the GGA ͑Al͒ and LDA ͑Si͒ are considered to work well.Furthermore, we include a transition metal, Pt, as a more complex material.Established bulk properties are examined to make sure the new functionals do not significantly worsen established results.Then vacancy formation energies are studied, a property known to include strong surface effects and which none of the presently established functionals describe correctly.No other functional has been initially tested on this intricate property.
Bulk properties only include weak surface effects.The equilibrium lattice constant a 0 and bulk modulus B 0 = ͉ − V‫ץ‬ 2 E / ‫ץ‬V 2 ͉ V 0 are obtained from the energy minimum given by a fit of seven points in a range about ±10% of the cell volume at equilibrium V = V 0 to the Murnaghan equation of state. 20As seen in Table II our functionals improve on the results of other functionals.A convincing sign of general improvement is the tendency for values to stay between the LDA and PBE, as they are known to overbind and underbind, respectively.As a measure of overall performance, the table shows the mean absolute relative error x ¯and its standard deviation = ͓͚͑x i − x ¯͒2 / N͔ 1/2 for N absolute relative errors x i .The value of gives the spread of the errors independently of their overall magnitude.If further testing confirms the LDA-LAG͑LAA͒'s robustness to be universal for solid-state systems, they should be considered as a "first  choice" for such applications.Furthermore, an explicit trend is seen in the sequence LDA, LAG, and LDA-LAG͑LAA͒.Throughout the table LAG shifts LDA values towards the PW91/PBE values, while LDA-LAG͑LAA͒ corrects them back towards ͑and occasionally even beyond͒ the LDA.This behavior illustrates the importance of compatible correlation.We now turn to tests of the strong surface effects manifest in calculations of the monovacancy formation enthalpy H V F = E V − ͑N −1͒E / N, where E V and E are total energies for the system with and without a vacancy, and N is the number of atoms in the fully populated supercell.Monovacancy energies are calculated using 64-atom cells.The vacancy cell is geometrically relaxed, and both vacancy and bulk cells are volume relaxed.The number of k points used is 4 3 for Pt, 6 3 for Al, and 3 3 for Si.The Si calculations are for the T d structure. 18For Pt and Si the supercells are too small for the results to be directly compared to experiment but are sufficient to allow for comparison between functionals.
Strong surface effects are seen for Al and Pt, but not in Si.This is seen by the widely different results between functionals for the metals.Similar to the bulk properties, our surface correlation corrects LAG results in the right direction, but it is apparent that it is still too crude to give truly quantitative results.The surprisingly good LDA result for Al might draw some attention, but as has been pointed out before, 2 it is not reflected in any other property of Al and is thus coincidental.The unexpected discrepancy between PW91 and PBE monovacancy energies will be addressed in another publication. 21e examine only solid-state systems; we do not assess performance for atoms and molecules.However, a hint is provided by the atomic XC energies given from the allelectron calculations used for constructing PP's ͑cf.Table II͒.The present functionals give results close to the LDA, with a slight adjustment towards the PBE.For atoms, the PBE is expected to be more accurate than the LDA. 4

IV. CONCLUSIONS
In conclusion, we have presented two promising functionals for use in DFT calculations.The method of their construction is generic and could potentially be used with any local approximation to ⑀ ˆxc in the interior region.The locality criteria precludes using, e.g., the PBE for this region, 8 and the effect of a localized equivalent cannot be inferred from GGA results.We are working on a gradient-corrected interior functional and an improved surface correlation.The two varieties of edge treatment, LAG and LAA, behave similarly but we recommend the LAA based on its better behavior far outside the edge.

FIG. 1 .
FIG.1.Parametrizations of Airy exchange ⑀ ˆx Airy vs scaled spatial coordinate .The solid black line is the true Airy exchange from Eq. ͑3͒.The inset shows the difference between the parametrizations and the true exchange.Far outside the edge, the LAA is more accurate than the LAG due to the former's proper limiting behavior.

FIG. 2 .
FIG. 2. Local surface XC energy for the r s = 2.66 jellium surface.The main figure shows the quantity that integrates to the surface energy xc in ergs/ cm 2 .The upper inset shows the difference between the functionals and LDA.The lower inset shows the interpolation indices X. Integration gives in ergs/ cm 2 for LDA 1188, for LAG 1121, and for LDA-LAG͑LAA͒ the "exact" RPA+ value of 1214.

TABLE I .
2ellium XC surface energies in erg/ cm2.RPA+ values are from Ref.15and are taken as exact.The LDA-LAG and LDA-LAA functionals are created using a two-parameter fit to values for r s up to 4.00.

TABLE II .
Results of electronic structure calculations for materials exhibiting widely different properties; Al, a free-electron metal; Pt, a transition metal; and Si, a semiconductor.The LDA-LAG and LDA-LAA functionals are from this paper, Eq. ͑12͒.Values given as percent are relative errors as compared to experimental values.Values in boldface are mean absolute relative errors.The standard deviation of absolute relative errors is defined in the text.LDA-LAG͑LAA͒ are not fitted to any values shown in this table, but to jellium surface energies.1.35± 0.05 eV from Ref. 22. e 0.68± 0.03 eV from Ref. 2. f 3.6± 0.2 eV from Ref. 23.
a Reference 24.b Reference 2. c Reference 25. d