Accurate and efficient computation of the Kohn-Sham orbital kinetic energy density in the full-potential linearized augmented plane wave method Lin-Hui Ye Key Laboratory for the Physics and Chemistry of Nanodevices, Department of Electronics, Peking University, Beijing 100871, P.R. China (Dated: January 8, 2015) (cid:80) (cid:12) (cid:12)2 The Kohn-Sham orbital kinetic energy density τσ(r) = iwiσ(cid:12)∇ψiσ(r)(cid:12) is one fundamental quantityforconstructingmeta-generalizedgradientapproximations(meta-GGA)forusebydensity functional theory. We present a computational scheme of τ (r) for full-potential linearized aug- σ mented plane wave method. Our scheme is highly accurate and efficient and easy to implement to 5 existing computer code. To illustrate its performance, we construct the Becke-Johnson meta-GGA 1 exchange potentials for Be, Ne, Mg, Ar, Ca, Zn, Kr, Cd atoms which are in very good agreement 0 withtheoriginalresults. Forbulksolids,weconstructtheTran-BlahamodifiedBecke-Johnsonpo- 2 tential(mBJ)andconfirmitscapabilitytocalculatebandgaps,withthereportedbadconvergence of the mBJ potential being substantially improved. The present computational scheme of τ (r) n σ should also be valuable for developing other meta-GGA’s in FLAPW as well as in similar methods a J utilizing atom centered basis functions. 7 PACSnumbers: 71.15.Mb,77.22.Ej ] i c I. INTRODUCTION which involves only the total spin density ρ (r), the s σ - Kohn-Shameffectivepotentialvσ(r)andtheKohn-Sham rl Nowadays, density functional theory1 has become the orbital eigenvalues εiσ. Eq.(3) is the form being imple- mt dominant method for electronic structure calculations. mented. It is more efficient than Eq.(1) since the depen- Accompanying its 50 year development, the quality of dence on individual orbital has been removed. . t the approximated energy functionals are also constantly Nevertheless, it is often unnoticed that the conver- a m improved.2 In order of higher accuracy, Perdew has cat- sion from (1) to (3) is not always valid, since typically egorized various functionals in terms of the “Jacob’s not all electrons states are solved by Kohn-Sham equa- - d ladder”:3 The first rung of the ladder is the local den- tion. In fact, core states are usually solved with rel- n sity approximation (LDA) which depends only on the ativistic corrections which is the case in the standard o spindensities(ρ ,ρ ). Thesecondrungconsistsofgener- FLAPW implementation. Then, expression (3) is even ↑ ↓ c alized gradient approximations (GGA) which depend on not non-negative-definite. Invalid numeric values are [ (ρ↑,ρ↓,∇ρ↑,∇ρ↓). Thethirdrungconsistsof“metagen- mostlyfoundaroundnucleiwheretheDirac-Kohn-Sham 1 eralized gradient approximations (meta-GGA)” which solutiondeviatesthemostfromtheKohn-Shamsolution. v further include the Kohn-Sham orbital kinetic energy Technically,theseinvaliddatamaybezero’edwithoutaf- 0 density τ (r), so that the functionals now depend on fectingthebandstructureverymuch,forwhichthedom- σ 9 (ρ ,ρ ,∇ρ ,∇ρ ,τ ,τ ). inant effects come from potential in the chemical bond- 3 ↑ ↓ ↑ ↓ ↑ ↓ As a fundamental component of meta-GGA, τ (r) ing regions. It is unavoidable, however, that the conver- 1 σ 0 should be implemented in a way which is highly accu- genceofself-consistencyisoftenruined,andsomeauthor . rateandefficientatthesametime. Forthefull-potential questionswhetherthebadconvergenceisintrinsictothe 1 linearized augmented planewave method4–6 a previous meta-GGA exchange potential itself.8 0 5 implementation of τσ(r) exists7. By definition τσ(r) re- Ofcourse, onecanintentionallyuseKohn-Shamequa- 1 quires the gradients of all occupied Kohn-Sham orbital tiontosolvenotonlythevalencestatesbutalsothecore : (w is the occupation number): states. With the sacrifice of the relativistic corrections, v Eq.(3) becomes equivalent to Eq.(1). Even so, however, i (cid:88) (cid:12) (cid:12)2 X τ (r)= w (cid:12)∇ψ (r)(cid:12) (1) the equivalence is purely mathematical rather than nu- σ iσ(cid:12) iσ (cid:12) r i=1 merical. This is because in almost all cases Kohn-Sham a Since direct computation of the gradients on dense equation is solved only approximately by use of finite meshes in the unit cell would be too costly, the authors basis sets. Through the conversion from (1) to (3) the make use of the Kohn-Sham equation incomplete basis set error is brought into the values of τ (r). (cid:26) ∇2 (cid:27) σ − 2 +vσ(r) ψiσ(r)=εiσψiσ(r) (2) theItoirsigainlwaalydsefidnesitiiroanblEeqt.o(1c)alscinuclaetethτisσ(erq)udatirioecntliys vfraolmid to convert τ (r) to the following form regardless of the details of how the orbital are solved. σ The purpose of this work is to present a computational τσ(r)= 12∇2ρσ(r)−vσ(r)ρσ(r)−(cid:88)εiσρiσ(r) (3) spclhanemeweavoef τmσe(trh)ofdor(FfuLlAl-pPoWte)n.tOiaulrlisncehaermizeedisahuiggmhleyntaecd- i=1 2 curate and efficient and requires minimum change to the our work in Section V. Within the FLAPW framework, existing computer code. Moreover, the increased accu- theunitcellispartitionedintoatomicspheres(calledthe racy of τ (r) substantially improves the convergence of muffin-tin region) and the space in between (called the σ self-consistency. interstitialregion). Implementationofτ (r)intheinter- σ This paper is organized as the following: In section II stitial region is simple and straightforward. Therefore, the FLAPW method is briefly reviewed. As this method most discussions are focused on the muffin-tin region. has become very popular in recent years and there have Throughout this work atomic units are used. To sim- been numerous literatures about it,9 our introduction is plify the notations, spin index and some other symbols limited to the minimum content which is relevant to the are suppressed wherever confusion is avoided. present work. In section III, we derive the formulae of our computational scheme of τ (r) for both the valence σ and the core states. The performance of this scheme is II. BRIEF INTRODUCTION TO THE FLAPW illustrated in section IV. We first construct the Becke- METHOD Johnson meta-GGA exchange potential10 for atoms, and then the Tran-Blaha modified version of the potential11 The FLAPW method solves the valence states and for bulk materials, by which we discuss its major feature core states by different techniques. Valence states are forcalculatingsemiconductorbandgaps. Wesummarize expanded by the following basis functions:  (cid:20) (cid:21) (cid:80) a (g)u (r)+b (g)u˙ (r) Y(cid:63) (gˆ)Y (ˆr) (muffin-tins) φ (r)= l l l l lm lm (4) g lm  √1 eig·r (interstitial) Ω Eachbasisfunctionisofahybridform: Insidethemuffin- muffin-tinregion. Inparticular,withinanatomicsphere, tins it takes the linear combination of atomic orbital, a lattice harmonics is constructed as while within the interstitial it is simply plane wave. In (cid:88) K (ˆr)= c (m )Y (ˆr) (6) thisfashion,thebasisfunctionbestaccountsforthegen- ν ν ν lνmν eral feature of the wave functions in the whole space mν which behave like atomic wave functions close to nuclei so that it is invariant under any local point group sym- and become free-electron-like in-between. The optimal metry R : i shape of the basis function allows the expansion of the valence state RiKν(ˆr)=Kν(ˆr) (7) (cid:88) Expand the density and the potential by the lattice har- ψ (r)= z (G)φ (r), (g=k+G) (5) nk nk g monics G (cid:88) ρ(r)= ρ (r)K (ˆr) (8) to use small cutoff for the reciprocal lattice vector G. ν ν In(4),Y issphericalharmonicswhichiswidelyused ν lm (cid:88) in atom centered basis functions. The u (r) is solved by v(r)= v (r)K (ˆr) (9) l ν ν the radial differential equation with a pre-chosen energy ν parameter and keeping only the spherical component of symmetry requirements are automatically fulfilled. the potential. The combination u (r)Y (ˆr) is called a l lm By applying R to Eq.(3), we realize that τ(r) has the muffin-tin function. The other radial function u˙ (r) is i l same symmetry as ρ(r) and v(r), and therefore can also the energy derivative of u (r). It serves as the first order l be expanded by the same set of lattice harmonics correction to u (r), so that the basis function becomes l applicable to a wide energy range. The use of the two (cid:88) τ(r)= τ (r)K (ˆr) (10) ν ν muffin-tinfunctionsinsideoneatomicsphereisproposed ν byAnderson5asthelinearisedapproximationtotheorig- inalaugmentedplanewavemethodofSlater4. Thelinear The expansion coefficients are calculated by: coefficients a (g) and b (g) are so chosen that the ba- (cid:90) l l sis function and its first derivative are both continuous τν(r)= τ(r)Kν(cid:63)(ˆr)dΩ across the sphere boundary. Toaccountforspacegroupsymmetry, densityandpo- Core states in the standard FLAPW implementation tentialareexpandedbysymmetrizedplanewaves(called are solved by Dirac-Kohn-Sham equation: starfunctions)intheinterstitialregionandsymmetrized (cid:110) (cid:111) spherical harmonics (called lattice harmonics) in the cα·p+c2β+v(0)(r) ψκµ(r)=εκµψκµ(r) (11) 3 to better account for relativistic effects. In (11), α and III. IMPLEMENTATION SCHEME OF τ(r) FOR β are the standard Dirac matrices. v(0)(r) only contains THE FLAPW METHOD the spherical component of v(r). The four component solution to Eq.(11) is A. τv(r) of the muffin-tin region: contribution from the valence states (cid:32) g (r)χ (ˆr) (cid:33) iκ κµ ψκcµ(r)= (12) Within the muffin-tin region, the expansion of the va- if (r)σˆ χ (ˆr) iκ r κµ lence state can be written as: (cid:20) (cid:21) withg (r)andf (r)beingtheradialpartofthemajor (cid:88) nκ nκ ψ (r)= A (nk)u (r)+B (nk)u˙ (r) Y (ˆr) (14) nk lm l lm l lm and minor components, respectively. Note that the core lm statesareassumedtobecompletelyrestrictedwithinthe muffin-tins so that their wave functions have no intersti- in which we have defined the coefficients tial part. Besides, g (r) and f (r) are not expanded (cid:88) nκ nκ A (nk)= z (g)a (g)Y(cid:63) (gˆ) bybasisset, butaredirectlyintegratedontherealspace lm nk l lm grid. G (cid:88) The total kinetic energy density is contributed by the Blm(nk)= znk(g)bl(g)Yl(cid:63)m(gˆ) core states and all occupied valence states: G The total kinetic energy density of the valence states in- τ(r)=τc(r)+τv(r) (13) side the muffin-tin region can be written as: (cid:88) (cid:12) (cid:12)2 (cid:88) (cid:26) (cid:27) τv(r)= w (cid:12)∇ψ (r)(cid:12) =∇2ρv(r)− w ψ(cid:63) (r)∇2ψ (r)−ψ (r)∇2ψ(cid:63) (r) (15) nk(cid:12) nk (cid:12) nk nk nk nk nk nk nk Note that to reach at (15) we did not make use of any On the right hand side of Eq.(15), the first term is physics equation such as Kohn-Sham equation. There- already calculated by GGA. The remaining two terms fore, (15) is equivalent to (1) regardless of the details of can be derived by use of the following relation: how the orbital are solved. (cid:110) (cid:111) (cid:18) 2W(cid:48)(r) l(l+1)W(r)(cid:19) ∇2 W(r)Y (ˆr) = W(cid:48)(cid:48)(r)+ − Y (ˆr) (16) lm r r2 lm Eq.(16)isafavorablepropertyofallmuffin-tinfunctions, for (15), we get τv(r) in the muffin-tin region. Projected by which the laplacian operation only affects the radial to the lattice harmonics representation, the final expres- part while leaves the angular part unaltered. Using (16) sion of the expansion coefficients is: τv(r)=(cid:88)(cid:40)(cid:20)u(cid:48)u(cid:48) + l(l+1)+l(cid:48)(l(cid:48)+1)−lν(lν +1)u u (cid:21)(cid:88)w (cid:88)A(cid:63) A (cid:68)Y (cid:12)(cid:12)Y (cid:12)(cid:12)Y (cid:69) ν l l(cid:48) 2r2 l l(cid:48) nk l(cid:48)m(cid:48) lm lm(cid:12) lνmν(cid:12) l(cid:48)m(cid:48) ll(cid:48) nk mm(cid:48) +(cid:20)u˙(cid:48)u˙(cid:48) + l(l+1)+l(cid:48)(l(cid:48)+1)−lν(lν +1)u˙ u˙ (cid:21)(cid:88)w (cid:88)B(cid:63) B (cid:68)Y (cid:12)(cid:12)Y (cid:12)(cid:12)Y (cid:69) l l(cid:48) 2r2 l l(cid:48) nk l(cid:48)m(cid:48) lm lm(cid:12) lνmν(cid:12) l(cid:48)m(cid:48) nk mm(cid:48) +(cid:20)u(cid:48)u˙(cid:48) + l(l+1)+l(cid:48)(l(cid:48)+1)−lν(lν +1)u u˙ (cid:21)(cid:88)w (cid:88)A(cid:63) B (cid:68)Y (cid:12)(cid:12)Y (cid:12)(cid:12)Y (cid:69) l l(cid:48) 2r2 l l(cid:48) nk l(cid:48)m(cid:48) lm lm(cid:12) lνmν(cid:12) l(cid:48)m(cid:48) nk mm(cid:48) +(cid:20)u˙(cid:48)u(cid:48) + l(l+1)+l(cid:48)(l(cid:48)+1)−lν(lν +1)u˙ u (cid:21)(cid:88)w (cid:88)B(cid:63) A (cid:68)Y (cid:12)(cid:12)Y (cid:12)(cid:12)Y (cid:69)(cid:41) (17) l l(cid:48) 2r2 l l(cid:48) nk l(cid:48)m(cid:48) lm lm(cid:12) lνmν(cid:12) l(cid:48)m(cid:48) nk mm(cid:48) 4 For comparison, we also write out the expansion coefficients of the density: (cid:88)(cid:40)(cid:20) (cid:21)(cid:88) (cid:88) (cid:68) (cid:12) (cid:12) (cid:69) ρv(r)= u u w A(cid:63) A Y (cid:12)Y (cid:12)Y ν l l(cid:48) nk l(cid:48)m(cid:48) lm lm(cid:12) lνmν(cid:12) l(cid:48)m(cid:48) ll(cid:48) nk mm(cid:48) (cid:20) (cid:21)(cid:88) (cid:88) (cid:68) (cid:12) (cid:12) (cid:69) + u˙ u˙ w B(cid:63) B Y (cid:12)Y (cid:12)Y l l(cid:48) nk l(cid:48)m(cid:48) lm lm(cid:12) lνmν(cid:12) l(cid:48)m(cid:48) nk mm(cid:48) (cid:20) (cid:21)(cid:88) (cid:88) (cid:68) (cid:12) (cid:12) (cid:69) + u u˙ w A(cid:63) B Y (cid:12)Y (cid:12)Y l l(cid:48) nk l(cid:48)m(cid:48) lm lm(cid:12) lνmν(cid:12) l(cid:48)m(cid:48) nk mm(cid:48) (cid:20) (cid:21)(cid:88) (cid:88) (cid:68) (cid:12) (cid:12) (cid:69)(cid:41) + u˙ u w B(cid:63) A Y (cid:12)Y (cid:12)Y (18) l l(cid:48) nk l(cid:48)m(cid:48) lm lm(cid:12) lνmν(cid:12) l(cid:48)m(cid:48) nk mm(cid:48) Eqs.(17) and (18) are similar. Their differences are lim- labor. ited to the radial parts while their angular parts are the same. This is not a coincidence, but a consequence of the property (16) of the muffin-tin functions. Eq.(17) is B. τc(r) of the muffin-tin region: contribution from the expression being implemented in the present work. the core states It is highly accurate because it is equivalent to Eq.(1) and no extra error is introduced throughout the deriva- Corestatesalsocontributetotheτ(r)ofthemuffin-tin tion. Moreover, the expressoin of the lattice harmonics region. Although core states are in the four-component expansion also ensures the correct symmetry of τ(r). form, each row is a separate muffin-tin function so that Eq.(17) is equally efficient because τv(r) can be cal- ν Eq.(16) is equally applicable. At ground state, all core culated together with ρv(r). During the construction of ν sub-shells are fully filled. Consequently, core density is thevalencedensity,almostallcomputationalcostof(18) spherical and τ(r) also is because they must have the is spent on the angular part, while the radial part takes same angular parts: no more than a few percent. By calculating τv(r) to- ν gether with ρv(r), the costly angular part needs to be τc(r)=τc(r) (19) ν 0 evaluated only once. The extra cost for the radial part of (17) is very small. Note that the radial part of (17) The expansion coefficient of the first lattice harmonics has no orbital dependence since it only requires the ra- K (ˆr) is contributed by all occupied sub-shells: 0 dial basis functions u (r) and u˙ (r) which are universal l l (cid:88) for all valence states. Technically, the implementation of τc(r)= τc (r) (20) 0 0,l the radial part of (17) to the construction of the valence l density is also straightforward and requires little human For the l=κ sub-shell: (cid:40) (cid:41) 2l l(l+1) l(l−1) τc (r)= √ g(cid:48) 2+ g2 +f(cid:48) 2+ f2 (21) 0,l 4π iκ r2 iκ iκ r2 iκ For the l=−κ−1 sub-shell: (cid:40) (cid:41) 2(l+1) l(l+1) (l+1)(l+2) τc (r)= √ g(cid:48) 2+ g2 +f(cid:48) 2+ f2 (22) 0,l 4π iκ r2 iκ iκ r2 iκ Since the radial solution of Dirac-Kohn-Sham equation thevalenceexpression(17)inwhichoneonlyneedstoset is slightly irregular, τ(r) diverges at the nucleus which l =0andl=l(cid:48) andremovethelinearizationtermu˙ (r). ν l is similar to that of all GGA potentials. Alternatively, C. τv(r) of the interstitial region core states can also be solved by Kohn-Sham equation. Thentheexpressionofτ0c(r)canbedirectlyderivedfrom Intheinterstitialregion,onlyvalencestatescontribute to τ(r): 5 1 (cid:88) (cid:88) (cid:16) (cid:17)(cid:16) (cid:17) τv(r)= w k+G k+G z(cid:63) (G )z (G )ei(Gj−Gi)·r (23) Ω nk i j nk i nk j nk GiGj Eqs.(17), (21), (22) and (23) have been implemented totheFLAPWcodeofNorthwesternUniversity.6 Totest 0 0 -2 the performance of our computational scheme, we next -1 cfoornsattroumctsthaendBebcuklek-Jsoehmniscoonndtyupcetomrse.ta-FGoGr Aatopmotse,ntcioarles σVx(au)--32 σVx(au)---864 states are solved by Kohn-Sham equation for compar- -4 Be -10 Ne ison with the earlier work. For semiconductors, both 0.01 0.1 1 10 0.01 0.1 1 10 r(au) r(au) Kohn-Sham equation and Dirac-Kohn-Sham equation are tested for solving the core states, and we find that 0 0 -2 the results as well as the speed of convergence are nearly -5 tahreessaomlveed. TbyheKreofhorne-,Sthharmouegqhuoauttiotnh.isworkallcorestates σVx(au)---864 σVx(au)-10 -10 Mg -15 Ar -12 -20 0.01 0.1 1 10 0.01 0.1 1 10 IV. THE BECKE-JOHNSON TYPE META-GGA r(au) r(au) POTENTIALS 0 0 -5 -5 -10 A. The atomic potentials σVx(au)--1150 σVx(au)---221550 TheBecke-Johnson(BJ06)potential10 isameta-GGA -20 Ca -30 Zn exchangepotentialattemptingtoapproachtheoptimized 0.01 0.1 1 10 -35 0.01 0.1 1 10 effective potential (OEP)12 of atoms. BJ06 consists of r(au) r(au) 5 two terms: 0 0 -5 -10 vxB,Jσ06(r)=vxB,Rσ89(r)+ π1(cid:114)152(cid:115)ρτσσ((rr)) (24) σVx(au)----22115050 σVx(au)--2300 -30 -40 -35 Kr -50 Cd The first term is the Becke-Roussel exchange potential13 -40 0.01 0.1 1 10 0.01 0.1 1 10 which is derived by assuming hydrogen-like exchange r(au) r(au) hole. vBJ06(r) is determined by a nonlinear equation x,σ involving ρ (r),∇ρ (r),∇2ρ (r) and τ (r), and closely FIG. 1. Becke-Johnson exchange potential for Be, Ne, Mg, σ σ σ σ Ar,Ca,Zn,Kr,CdatomscalculatedbytheFLAPWmethod resembles the Slater averaged exchange.14 Compared to with the implementation of τ (r) of this work. OEP, vBR89(r) is too deep and lacks the characteristic σ x,σ shell-structures. Therefore, the second term is added for correction. This term is repulsive which reduces the strength of the exchange. Besides, it strongly varies bottom also agrees well with the original results. Over- across atomic shells which restores the OEP shell struc- all, all atomic potential are in excellent agreement with tures. Summing up the two terms, the BJ06 potential is Ref.10. The remaining slight differences may be ascribed close to the accuracy of the OEP of atoms, yet is much to the use of different boundary conditions in the two simpler since it is directly constructed from local quan- groupsofcalculations: Whiletheoriginalcalculationsare tities rather than being solved by the complicated OEP performed with the molecular code, the present calcula- equation. Nevertheless,BJ06isamodelexchangepoten- tions use the FLAPW code with supercells and periodic tial without corresponding exchange energy. Therefore, boundary condition. it cannot be used for calculating total energy. For Be, Ne, Mg, Ar, Ca, Zn, Kr, Cd atoms, the BJ06 potentials have been published by Becke and Johnson.10 We have generated the same potentials by our FLAPW B. The bulk potentials code and the results are plotted in Figure 1. As can be seen,boththelocationandtheheightofthe“bumps”be- tweenatomicshellsarewellreproduced. Besides,closeto To test our implementation of τ(r) for bulk materi- the nuclei all potentials have the correct, asymptotically als, we construct the Tran-Blaha modified version of the flat shape, indicating the adequacy of the treatment of BJ06 potential (mBJ)11 for solids. The TB09 exchange τ(r) in the muffin-tin region. The depth of the potential is essentially the same as BJ06, with only the weights of 6 the BR89 term and the shell term being changed: (cid:114) (cid:115) 1 5 τ (r) vTB09(r)=cvBR89(r)+(3c−2) σ (25) x,σ xσ π 12 ρ (r) σ In Eq.(25) the weights are determined by the c parame- ter, which is essentially adjustable, but can also be cal- culated “pseudo-ab initio-ly” by empirical formulae pro- posed by the authors15. For c = 1 TB09 reverts to the original BJ06. The total mBJ potential is formed either byTB09exchangealone,orbycombiningitwithacorre- lationpotentialsuchasLDA.Itissuchslightadjustment oftheweightswhichleadstothesurprisingdiscoverythat the mBJ potential is capable of calculating semiconduc- tor band gaps. Albeit being a local potential so that its computational cost is essentially maintained at the LDA level, in many cases the mBJ potential can achieve band gapaccuracywhichisevencomparabletothemuchmore FIG. 2. (Color online) Rubber-sheet plot of the up-spin ki- sophisticated GW approximation. netic energy density log10τ↑ of the NiO (001) plane. Atomic Applications of the mBJ potential has been quickly positions are illustrated at the left bottom of the plot. growing in recent years which also motivates the present work. However, the existing FLAPW implementation uses Eq.(3) to calculate τ(r), and therefore suffers from the problems mentioned in the Introduction. It is thus desirable to recheck the mBJ potentials by our imple- mentation of τ(r). In all our calculations, the c parame- ter of (25) is automatically determined by the suggested formula:15 c=A+Bg¯ (26) with A=0.488, B =0.5, and 1 (cid:90) 1(cid:18)|∇ρ↑(r)| |∇ρ↓(r)|(cid:19) g¯= + d3r (27) V 2 ρ↑(r) ρ↓(r) cell cell For the correlation potential we use LDA.16 One differ- encewiththeearlierworkisthatinourcalculationsspin- orbit coupling (SOC) is included perturbatively through the second variational approach. Typically, the effect of FIG.3. (Coloronline)Rubber-sheetplotoftheup-spinmBJ SOConthevalueofbandgapissmallerthan0.1eV.But potential v↑mBJ of the NiO (001) plane. Atomic positions are illustrated at the left bottom of the plot. for materials (such as ZnTe) with heavy elements SOC can change the band gap by as much as 0.3 eV. In Figures 2 and 3 we plot17 the spin-up kinetic en- ergy density and the mBJ potential of the NiO (001) first derivative while the meta-GGA require the second plane for comparison with the similar plot of Ref.18. At derivative of the density. In fact, the basis function is thenucleiτ(r)showssharpspikeswhilevmBJ(r)achieves already discontinuous even at the zero’th order because the minimum, although neither quantities diverges. The in Eq.(4) the summation is cutoff at finite (lm). By the present implementation of τ(r) allows us to reveal even same reason, the GGA potential19 is also discontinuous thefinestdetailsofthemBJpotentialinFigure3. Espe- across sphere boundary, although we notice that usually cially, all shell structures are clearly seen. Another fea- the discontinuity of the meta-GGA is more serious than ture is that, although τ(r) and vmBJ(r) are very smooth the GGA. within both the muffin-tin and the interstitial regions, With our highly accurate τ(r), we have calculated the acrossthesphereboundarytheyareobviouslydiscontin- mBJ band gaps for all semiconductors in the SC40 test uous. Suchdiscontinuityisalsoobservedinthepotential set.20 The results are plotted in Figure 4 together with plot of NiO in Ref.18 and is irrelevant to the potential theLDAbandgapstocomparewithexperiment. Indeed, itself or to the implementation of τ(r). Rather, it is for all semiconductors mBJ systematically improves the caused by the discontinuity of the basis functions. By bandgaps. Especially,theimprovementisnearlyasgood (4), the FLAPW basis function is continuous only to the forthesmallestbandgap(ofInSbwithE =0.23eV)as g 7 8 potential 50% more iterations are needed to achieve self- 7 LDA(withSOC) consistencythanLDAorGGA.Foranexample,forMoS2 mBJ(withSOC,autoc) LDAtakes22iterationstoconverge,whilemBJwithour 6 implementation of τ(r) requires 32 iterations. Neverthe- (eV)5 lweshsi,chwewehafvoeunadlsothtartietdhiemsppleeemdeonfticnognτve(rrg)ebnyceEiqs.(d3r)a,sbtiy- p ga4 cally ruined. For MoS2 the mBJ calculation with Eq.(3) al requires 79 iterations to converge, while calculation of c eti3 NiO does not even converge at all. In fact, in all our heor2 testings, the implementation of τ(r) by (1) always con- T verges better than by (3), and the speed of convergence 1 is not sensitive to whether the core states are solved by Kohn-ShamequationorDirac-Kohn-Shamequation. Be- 0 sides,formorethan70solidswehavecalculatedwith(1), 0 1 2 3 4 5 6 7 8 all converge well except for CoO and FeO which fail be- Experimentalgap(eV) cause mBJ is incapable of lifting the near degeneracy of the 3d states. Therefore, it is safe to conclude that at FIG. 4. (Color online) Comparison between the mBJ band least a large part of the reported convergence problem is gaps,theLDAbandgapsandtheexperimentalbandgapsfor not intrinsic to the potential but is caused by technical the 40 semiconductors of the SC40 test set. reasons. for the largest (of MgO with E =7.31 eV). Besides, we g V. SUMMARY also confirm that in most cases the size for the band gap increases monotonically with increasing c. Nevertheless, seriouserrorsstillexistforthethreeBacompounds,BaS, Inthisworkwehavepresentedanhighlyaccurateand BaSe, BaTe, suggesting that there might be problems in efficient computational scheme for the Kohn-Sham or- the Ba potential or the experimental data. bital kinetic energy density τσ(r) to the full-potential Wehavethusfoundthattheimprovementofτ(r)does linearised augmented plane wave method. To test its not qualitatively change the most important feature of performance,wehaveconstructedtheBecke-Johnsonex- the mBJ potential. This is not surprising since the im- changepotentialforatomswhichareinverygoodagree- provement of τ(r) is most significant around the nuclei ment with the original results. For bulk solids we have while the band gap correction mainly depends on po- constructed the Tran-Blaha modified Becke-Johnson po- tential far away.18 In more details, the dominant contri- tential for semiconductors and confirmed its capability bution to the band gap correction comes from the shell to calculate semiconductor band gaps. As to the conver- term which is repulsive so that it pushes all band states genceproblemreportedbefore,wefoundthatalargepart upward. Away from the nuclei the shell term increases is due to the impropriate implementation of τ(r). With andforopensystemsitapproachesapositiveconstantat respect to accuracy, efficiency, easiness of implementa- r →∞. Because the conduction states are usually more tionandspeedofconvergence, ourschemeallsupersedes delocalizedthanthevalencestates,theupshiftofthecon- the existing implementation. We expect this work to be ductionstatesbytheshelltermislargerthanthevalence valuablefordevelopingothermeta-GGA’sinFLAPWas states, and therefore the band gap is enlarged. Through well as in similar methods utilizing atom centered basis this mechanism, band gap correction is not substantially functions. affected by the improvement of τ(r) around the nuclei. For the same reason, it is neither affected very much by the divergence of τ(r) if core states are solved by Dirac- VI. ACKNOWLEDGEMENT Kohn-Sham equation. Ithasbeenreportedbeforethatself-consistencybythe This work is supported by the Ministry of Science and mBJ potential is harder to achieve than LDA or GGA.8 Technology of China (Grant Nos. 2011CB933001 and This is confirmed in our work. Typically, using the mBJ 2011CB933002). 1 P.HohenbergandW.Kohn,Phys.Rev.136,B864(1964); 5 O.K. Anderson, Phys. Rev. B 12, 3060 (1975). W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965). 6 E. Wimmer, H. Krakauer, M. Weinert, A.J. Freeman, 2 A.D. Becke, J. Chem. Phys. 140, 18A301 (2014). Phys. Rev. B 24, 864 (1981). 3 J.P. Perdew and K. Schmidt, AIP Conf. Proc. 577, 1 7 F. Tran, P. Blaha and K. Schwarz, J. Phys.: Condens. (2001). Matter 19, 196208 (2007). 4 J.C. Slater, Phys. Rev. 92, 603 (1953). 8 H. Jiang, J. Chem. Phys. 138, 134115 (2013). 8 9 See, for example, D.J. Singh and L. Nordstro¨m, 16 S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, Planewaves, Pseudopotentials and the LAPW Method,2nd 1200 (1980). Edition, Springer. 17 TheplotsaremadebytheOpenDXvisualizationpackage. 10 A.D.BeckeandE.R.Johnson,J.Chem.Phys.124,221101 See http://www.opendx.org/. (2006). 18 D.Koller,F.Tran,andP.Blaha,Phys.Rev.B83,195134 11 F.TranandP.Blaha,Phys.Rev.Lett.102,226401(2009). (2011). 12 J.D. Talman and W.F. Shadwick, Phys. Rev. A 14, 36 19 J.P.Perdew,K.Burke,andM.Ernzerhof,Phys.Rev.Lett. (1976). 77, 3865 (1996). 13 A.D. Becke, M.R. Roussel, Phys. Rev. A 39 3761 (1989). 20 J.Heyd,J.Peralta,G.E.Scuseria,andR.Martin,J.Chem. 14 J.C. Slater, Phys. Rev. 81, 385 (1951). Phys. 123, 174101 (2005). 15 D.Koller,F.Tran,andP.Blaha,Phys.Rev.B85,155109 (2012).