Accurate and efficient computation of the Kohn-Sham orbital kinetic energy density in the full-potential linearized augmented plane wave method

2.4 MB·English

by  Lin-Hui Ye

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).