Implementation of an all-electron GW approximation based on the PAW method without plasmon pole approximation: application to Si, SiC, AlAs, InAs, NaH and KH

11 Pages·0.47 MB·English

by  S. Lebègue

Implementation of an all-electron GW approximation based on the PAW method
without plasmon pole approximation: application to Si, SiC, AlAs, InAs, NaH and KH
S. Leb`egue,1,2 B. Arnaud,3 M. Alouani,1,4 and P. E. Bloechl4,5
1Institut de Physique et de Chimie des Mat´eriaux de Strasbourg (IPCMS),
UMR 7504 du CNRS, 23 rue du Loess, 67037 Strasbourg, France, EU
3 2Department of Physics, University of California, Davis, California 95616
0 3Groupe Mati`ere condens´ee et Mat´eriaux (GMCM),
0 Campus de Beaulieu - Bat 11A 35042 Rennes cedex, France, EU
2 4Kavli Institute of Theoretical Physics, University of California, Santa Barbara, Ca 93111
n 5Institute of Theoretical Physics, Clausthal University of Technology,
a Leibnizstr. 10 D-38678 Clausthal Zellerfeld, Germany, EU
J (Dated: February 2, 2008)

7
A new implementation of the GW approximation (GWA) based on the all-electron Projector-
1
Augmented-Wavemethod(PAW)ispresented,wherethescreenedCoulombinteractioniscomputed

withintheRandomPhaseApproximation(RPA)insteadoftheplasmon-polemodel. Twodifferent
]
i ways of computing the self-energy are reported. The method is used successfully to determine the
c
quasiparticle energies of six semiconducting or insulating materials: Si, SiC, AlAs, InAs, NaH and
s
- KH.Toillustratethenoveltyofthemethodtherealandimaginarypartofthefrequency-dependent
l self-energytogetherwiththespectralfunctionofsiliconarecomputed. Finally,theGWAresultsare
r
t compared with other calculations, highlighting that all-electron GWA results can differ markedly
m
from those based on pseudopotential approaches.
.
t
a PACSnumbers: 71.15.-m,71.15.Mb,71.20.Nr
m
-
d I. INTRODUCTION tion the validity of the former approaches.
n However,someattemptshavebeenmadetogobeyond
o the plasmon-pole approximation.7,10,17,18,19,20,21 In par-
c For many weakly correlated materials, the density- ticular, Aryasetiawanhas approximately determined the
[ functionaltheory1(DFT)inthelocal-densityapproxima-
screening within the RPA using a linear muffin-tin or-
1 tion (LDA) provides a good description of their ground- bital method within the atomic sphere approximation7 (
v state properties. However, DFT is not able to describe LMTO-ASA).Thismethod,althoughfast,approximates
0 correctly their excited states. Thus, for example, the the space by atomic centered overlapping spheres, thus
2 band gaps in the LDA are typically much smaller than
completely neglecting the interstitial region, and hence
3 the experimental values. Quasiparticle (QP) electronic-
making the reliability of the GW method uncertain.
1
structure calculations beyond the DFT are therefore
0 Kotani and van Schilfgaarde based their full-potential
3 highly desirable. LMTOGWcalculation17ontheworkofAryasetiawanby
The GW approximation (GWA) of Hedin2,3, which
0 takingintoaccountcorrectlytheinterstitialregion. Nev-
/ produces a good approximation for the electron’s self- erthelesstheirmethodisnotquiteaccuratesinceintheir
t
a energy Σ, enables us to make first-principle QP calcula- implementation they didn’t take into account the multi-
m tions for realistic materials. Thus, the GWA has been
plicityofthesameangularmomentaforagivenprincipal
successfully applied to the calculation of QP electronic
- quantumnumberinthebasisset(likesimultaneouslyus-
d structures of many kinds of materials4,5,6,7. In particu-
ing the 3d and 4d states). Finally, Ku and Eguiluz pro-
n lar, recent success has been achieved on predicting the
ducedselfconsistentandnon-selfconsistentQPbandgaps
co metal-insulator transition of bcc hydrogen,8 electronic based on an approximate Luttinger-Ward functional,18
excitationsofyttriumtrihydride,9aswellastheQPelec-
: the non-selfconsistent results are much smaller than all
v tronic structure of copper.10 Unfortunately, most of the
existing GW calculations. Since these results are based
i
X GWA implementations are based on the pseudopoten- onadifferentschemewehavechosennottodiscusstheir
tialtypeofapproachestogetherwithplasmon-pole(PlP)
r methodfurther. Ontheotherhand,severalpseudopoten-
a models.11,12,13,14,15,16 Theweaknessofthesetypesofcal-
tial have produced GW results without resorting to the
culations is that the imaginary part of the self-energy is
plasmon-pole approximation. These methods, although
not accessible, making it impossible to determine spec-
interesting,usepseudowavefunctionsandhencecanonly
tral functions and hence to interpret photoemission ex-
determine pseudo-matrix elements of operators, making
periments. In addition the PlP approximation is ex-
themdifficulttojustifyasquantitativeandreliablemeth-
pected not to hold for systems with localized electrons.
ods for computing QP properties.
Moreover, it has been noticed recently15,17 that GWA
The major purpose of this paper is then to present
implementations basedonpseudopotentialmethods lead
a new implementation of the GWA method using
to larger and more k-dependent shifts than calculations
theall-electronfull-potentialProjectorAugmented-Wave
based on all-electron DFT methods, bringing into ques-
2
method(PAW)completebasisset,andwithoutusingany interaction W yields the successful GW approximation
PlP model for the determination of the dielectric func- for Σ. In this approximation, both the non-locality and
tion. The screening of the Coulomb interaction is thus the dynamical correlations are included. Assuming that
described in the Random Phase Approximation (RPA), the difference Σˆ −Vˆ between the self-energy and the
xc
avoiding further approximations. Kohn-Sham exchange and correlation potential is small,
The paper is organizedas follows. In section II we de- we can use a perturbation theory approach to solve the
scribe our implementationofthe GW approximation. In effective QP Hamiltonian Hˆqp
Section III we present our QP calculations for Si, SiC,
AlAs, and InAs and also for the alkali hydrides com- Hˆqp =Hˆ +(Σˆ −Vˆ ) (2)
KS xc
pounds NaH and KH (to our knowledge this is the first
GWA study of these alkali hydrides). At the end of this and determine the QP energies by expanding the real
section we compare and discuss our results with other partofselfenergytofirstorderaroundǫ (k)thusmaking
n
calculations and experiments. the comprison with the PlP models possible
ReE (k))=ǫ (k)+Z ×
n n nk
II. FORMALISM [hΨ |ReΣ(r,r′,ǫ (k))|Ψ i−hΨ |VLDA(r)|Ψ i](3)
kn n kn kn xc kn
A. The PAW method
where the QP renormalizationfactor Z is given by
nk
The PAW formalism has been well-described ∂
elsewhere,22,23,24,25 so we will not discuss it in this Znk =[1−hΨkn|∂ωReΣ(r,r′,ǫn(k))|Ψkni]−1. (4)
paper. The PAW method22 is a very powerful all-
electron method for performing electronic structure This assumptionis valid for simple sp bonded materials,
calculations within the framework of the LDA. It takes since it was shown that the QP wave function ψ and
kn
advantage of the simplicity of pseudopotential methods, Kohn-Shamwavefunction Ψ are almostidentical, i.e.,
kn
but describes correctly the nodal behavior in the aug- the QP Hamiltonian Hˆqp is diagonal in the Ψ basis
kn
mentation regions. The selfconsistent calculation of the for simple sp bonded semiconductors.11,12 We therefore
electronic structure is performed using the Car-Parinello assume this behavior for the materials studied in this
method over the occupied states. To determine the paper. According to this equation, the LDA eigenval-
eigenvaluesandeigenvectorsofallunoccupiedstates (up ues ǫ (k) are then corrected by the GW approximation.
n
to 200 eV above the top of the valence states) needed The numerical work is therefore considerably reduced,
for the GW calculations, we have extracted the selfcon- but still computationally demanding.
sistent full-potential, constructed and diagonalized the Inourimplementation,wehavecalculatedtheGreen’s
PAW Hamiltonian for every irreducible k point in the functiononlyforthevalenceandconductionstates. One
Brillouin zone. has then to subtract out only the valence exchange and
correlationpotentialinEq. (3). Tochecktheaccuracyof
this procedure, we have also used the so-called Hartree-
B. The GW approximation Fockdecoupling,15,26 andhavefoundthattheaverageer-
rorintheQPenergiesofSiwithrespecttothetopofthe
1. Quasiparticle energies valencestatesis0.05eV.Theapproximationmadehereis
theonecurrentlyusedinallpseudopotential-basedGWA
In general, the QP energies E (k) and wave function calculations, making our method compatible with exist-
n
ψ (r) are determined from the solution of the QP ing GW implementations.
kn
equation
(T +V +V )ψ (r)+ d3r′Σ(r,r′,E (k))ψ (r′)
ext h kn n kn 2. Screened Coulomb interaction
R
=E (k)ψ (r) (1)
n kn
Forthecalculationoftheself-energy,oneneedstoeval-
whereT is thefree-electronkinetic energyoperator,V uate the dynamically screened interaction W(r,r′,ω),
ext
theexternalpotentialduetotheioncores,V theaverage which can be rewritten in reciprocal space as:
h
electrostatic(Hartree)potential,andΣ the electronself-
eEnqe.rgy(1)opiserafitnodri.ngTahne amdaeqjouratdeiffiacpuplrtoyxicmonantieocntefdorwtihthe WG,G′(q,ω)=4π|q+1G|ǫ˜−G1,G′(q,ω)|q+1G′|. (5)
self-energy operator Σ(r,r′,E (k)). Nonetheless, it was
n
shownbyHedin2,3 thatwritingtheself-energyasaprod- The symmetrized dielectric matrix ˜ǫGG′(q,ω) is defined
uct of the Green’s function and the screened Coulomb in the random phase approximation (RPA) by27
3
8π
ǫ˜GG′(q,ω) = δGG′ − Ω|q+G||q+G′| MGvc(k,q)[MGvc′(k,q)]∗
vX,c,k
1 1
× − , (6)
(cid:18)ω+ǫ (k−q)−ǫ (k)−iδ ω−ǫ (k−q)+ǫ (k)+iδ(cid:19)
v c v c
with the following notation: integraloftheselfenergynumerically. Inourimplementa-
tion,weneedtocomputeǫ˜GG′(q,ω)alongtheimaginary
MGnm(k,q)=hΨk−qn|e−i(q+G)r|Ψkmi, (7) axis and for some real frequencies. This technical point
will become clearer in the next subsection.
where v and c denote, respectively, the valence and con- To reduce the computationalcostofthe GWA, we use
ductionstates,andδ apositiveinfinitesimal. Thematrix
symmetryproperties. Detailsabouttheutilizationofthe
elements given by Eq. (7) are evaluated using the PAW
symmetry for the static dielectric matrix has been al-
basis set as described in Ref. 15. readygivenelsewhere5,15,28,29,sowejustdescribebriefly
Most of GWA calculations use a kind of PlP approx-
how to use the symmetry in the case of the dynamical
imation. This is computationally efficient since one ob-
dielectric function. For the case of pure imaginary fre-
tains an analytic expression for the integral in the self-
quencies,wecouldsafelyignorethebroadeningfactoriδ;
energy. Itisnotclear,however,thatthiskindofapprox- inthiscase˜ǫGG′(q,iω)ishermitianandwecouldusethe
imation is valid for describing the QP of different kind
symmetry just as in the static case. We can then write
of materials. It is for this reason that we have chosen to
the symmetrized dielectric matrix as
avoidthe PlPmodelaltogether,tocompute the dynami-
caldielectricfunctionintheRPA(6),andtoperformthe
8π
ǫ˜GG′(q,iω) = δGG′ − Ω|q+G||q+G′| MRvcG(k,q)[MRvcG′(k,q)]∗
k∈XBZqXv,c RX∈Gq
1 1
× − , (8)
(cid:18)iω+ǫ (k−q)−ǫ (k) iω−ǫ (k−q)+ǫ (k)(cid:19)
v c v c
where Gq is the little group of the point group G such WG,G′(q,ω) only for irreducible points of the first Bril-
that Rq=q; R being a symmetry operation. The com- louin zone (BZ). We then determine easily the screened
putational cost is further reduced by noticing that interaction for all k-points in the Brillouin zone using
symmetry properties.
˜ǫGG′(Rq,iω)=ǫ˜R−1GR−1G′(q,iω). (9)
for real ω, although the dielectric matrix is not hermi-
tian, we could use the symmetry by making a decom-
3. Self Energy
position into hermitian and anti-hermitian parts of the
polarizability PGG′(q,ω). if we define AGG′(q,ω) and
BGG′(q,ω) by The selfenergyΣ is the keyquantity ofany GWAcal-
culation. As previously noticed, we have chosento avoid
AGG′(q,ω)= PGG′(q,ω)+2 PG†G′(q,ω) (10) pdleanscmeoonf-tphoelescmreoedneelsdainntdercaocmtiponutWe ΣwwitihthintthheeωRPdeAp.en-
First,wesplittheintegralofthe selfenergyintoabare
and
exchange or Hartree-Fock contribution ΣX and an en-
BGG′(q,ω)= PGG′(q,ω)−PG†G′(q,ω), (11) eenrgeyrgdyepeffenecdtesntbecyoonntrdibΣutXio.nTΣhCe(ωm)awtrhixicheledmesecnrtibseosfstehlfe-
2i
self-energy are now given by the sum of
then equations (8) and (9) still hold, allowing us
to perform the same computational tasks as for the 4π occ |Mmn(k,q)|2
symmetrized dielectric matrix with imaginary frequen- hΨ |ΣX|Ψ i=− G (12)
kn kn Ω |q+G|2
cies. This procedure makes it possible to first compute Xq Xm XG
4
where the summation is over occupied states, and
TABLEI:Calculatedquasiparticleenergiesofsiliconforsome
pointsintheBrillouin zonewithourtwodifferentimplemen-
hΨ |ΣC(ω)|Ψ i
kn kn tations. The results are in good agreement with each other.
In thelast line, theminimum band gap Eg is presented.
1
= Ω [MGmn(k,q)]⋆MGm′n(k,q) First methoda Second methodb
Xq GXG′Xm Γ1v −11.85 −11.87
× CGmG′(k,q,ω) (13) Γ25′v 0.00 0.00
Γ15c 3.09 3.09
with Γ2′c 4.05 4.06
Cm (k,q,ω) X1v −7.74 −7.68
GG′ X4v −2.90 −2.91
i WC (q,ω′) X1c 1.01 1.03
= dω′ GG′ (14X)4c 10.64 10.59
2π Z ω+ω′−ǫm(k−q)+iδsgn(ǫm(k−q)−µ)
L2′v −9.57 −9.50
where WC is defined as WC = W − v, with v being L1v −6.97 −6.90
the bare Coulomb potential. To evaluate this integral L3′v −1.16 −1.17
directly on the real axis one should compute WC for L1c 2.05 2.03
many points ω′ since the shape of WC along the real L3c 3.83 3.83
axis is rather ragged. Even though this has been done
by some authors,10 we choose to avoid this difficulty Eg 0.92 0.90
by using the fact that WC is well behaved along the aContour deformationmethod
imaginary axis. In the present work, we have performed bAnalyticcontinuation method
this integral using two different methods:
In the first one, the contour of the frequency integral
(14) is deformedina wayto obtainanintegralalongthe
imaginary axis plus contributions from the poles of the
analyticallycontinuedto the realaxis by fitting it to the
Green’s function. In this case, we obtain the following
following Pad´e form
expression:
Cn (k,q,ω)
GG′ a +a z+a z2+..+a zN
0 1 2 N
P(z)= (15)
= −π1 Z0∞dω′′WGCG′(q,iω′′)(ω−ωǫn−(kǫn−(kq)−)2q+) ω′′2. b0+b1z+b2z2+..+bMzM
±WGCG′(q,±(ω−ǫn(k−q)))θ(±(ω−ǫn(k−q))) where ai and bi are complex parameters that are deter-
×θ(±(ω−µ))θ(±(ǫ (k−q)−µ)) mined during the fit along the imaginary axis. Values
n
of 5 for N and of 6 for M provided us with an accu-
Thefirsttermrepresentsthecontributionalongtheimag- rate and stable fit. The same kind of continuation has
inary axis and is evaluated by Gaussian quadrature. alsobeenappliedwithsuccesstocomputethedynamical
The second is from the poles of the Green’s function response function,30,31 so we are confident of its relia-
and its computation is done by fitting values of WC bility. The main difference between the work presented
at ±(ω − ǫ (k − q)) from values on a given mesh of hereandthatofRef. (19)isthatourcodestartsfroman
n
frequencies54. HereµdenotestheFermilevelintheLDA all-electron basis, so we are not using fast Fourier trans-
andω′′isdefinedtobereal. Thismethodissimilartothe forms to switch between real and reciprocal spaces and
one used by Aryasetiawanfor the implementation of the between time and frequency domains. Our expression
GWA based on the LMTO method in the atomic sphere (15) is also different, but we believe that this is of minor
approximation7 (ASA), and within the GWA of Kotani importance55.
and coworkers based on the full-potential linear muffin- In both cases, the integration over the first Brillouin
tin orbital (FP-LMTO) method.17 The reader can find zone is done by the special-point technique.32 The num-
more details about this integration procedure in Refs. ber of bands as well as the number of G vectors in (13)
7,17. Similar work has been also carried out by Bechst- is increased until the QP energies are converged. Simi-
edt and coworkers20 as well as by Fleszar and Hanke21 larly, the number of frequency points ω′ for which WC
starting from a pseudopotential approach. is computed is increased until Cn (k,q,ω) is well con-
GG′
Inoursecondimplementation,whichissimilartothatof verged. The two different implementations allow us to
Ref. 19, we evaluate the matrix elements of the correla- check carefully our results, and as can be seen in Table I
tive part of the self-energy hΨ |ΣC(ω)|Ψ i for a set of for the case of silicon, the QP energies are insensitive to
mk lk
imaginary frequencies iω, the resulting quantity is then the method used to compute the self-energy.
5
4. Treatment of the Coulomb divergence
TABLEII: Lattice constants a(in atomic units)and energy
cutoffsEcut(inRydberg)usedforourPAWcalculations. The
The last point we wish to discuss is an additional dif- lattice constants are from Ref. 35, unless stated otherwise.
ficulty which occurs when evaluating the self-energy by
a summation of the q points over the full BZ. We can- a Ecut
Si 10.26 20
not apply the special-point technique directly since the
SiC 8.24 25
integrands have a 1/q2 singularity for q → 0 as can be
AlAs 10.67 20
seen for example by putting G = 0 in the expression InAs 11.41 20
of the exchange term given by (12). The difficulty can NaH 9.28a 40
be removedby adding and subtracting a term which has KH 10.83a 40
thesamesingularitiesastheinitialexpression,andwhich
canbe evaluatednumericallyandanalytically. As acon- aReference36
sequence, the integrals over the BZ are rewritten
G(q)= [G(q)−A F(q)]+A F(q), (16) toprovideus withthe coordinatesofthe new 6points in
Xq Xq Xq the BZ.56 The computational cost is further reduced by
finding equivalent points among those 6 points, and we
whereF(q)isanauxiliaryperiodicfunctionthatdiverges end up with wellbehavedand easily evaluatedBZ sums.
like 1/q2 as q vanishes. The term is regular and can be
evaluatedbythespecialpointtechniquewhereasthelast
sum is evaluated analytically. For the exchange term, it III. NUMERICAL RESULTS AND DISCUSSION
isnotdifficulttoevaluateAin(16),butitgetsmorecom-
plicated for the correlative part of the self-energy (13).
In this section we present our theoretical quasiparti-
ThepurposeoftheoffsettedΓ-pointmethod17 istoavoid
cle energies for the six materials studied in this paper,
the evaluation of the quantity A, but still to be able to
together with the available theoretical and experimen-
dealwith the divergence. The main idea is to find a new
tal results. In Sec. IIIA we report our results for four
mesh of points such that
semiconductors(Si,SiC,AlAs, InAs)ofzinc-blendetype
structure, while Sec. IIIB is devoted to studying the al-
F(q)= F(q′) (17)
kali hydrides NaH and KH in the rocksalt phase. Table
Xq Xq′ II presents the experimental lattice parameters and the
energy cut offs E used for the final convergedcalcula-
wheretheΓ-pointisincludedintheoldmeshqbutnotin cut
tions.
the new one q′: the Γ-point is replaced by other points
(different from Γ) to construct the q′ grid in order to
fulfill Eq. (17). Eq. (16) is therefore rewritten as:
A. Results for Si, SiC, AlAs, and InAs
G(q)= [G(q)−A F(q)]+A F(q′), (18)
Xq Xq Xq′ Siliconisprobablythemostcarefullystudiedsemicon-
ductor, and several GWA results are available. Using
Then we show by inspection that the term silicon as a prototype will allow us to test our method
[G(q)−A F(q)] is equal to [G(q′)−A F(q′)] by making careful comparisons with previous GWA cal-
q q′
wPith a controlled error,Eq. (18) tPransforms to culations. As mentioned earlier, our code presents two
different ways for calculating the self-energy. We there-
G(q) = G(q′) (19) fore test their accuracy for silicon in Table I. We find
Xq Xq′ that the results of the two methods are almost identical,
showing that they are equally reliable for computing the
because the two terms which contain the function F(q) self-energy. In particular, our implementation with the
cancel out since they are evaluated on the same q′ grid. extrapolation procedure makes it possible to represent
We have therefore avoided the evaluation of the compli- the full-frequency dependence of the self-energy with a
cated A term in Eq. (16). small additional computational cost. Fig. 1 shows the
Theremainingpointstobeaddressedarethechoiceof real and imaginary parts of Σ along the real axis for sil-
the function F(q) and the number of additional points icon at the Γ point for a wide range of frequencies. The
forthenewmeshusedtosolve(17). Inourcase,wewrite agreement with previous work is excellent.19
A special feature of our work is the possibility of
exp(−|q+G|2) obtaining the imaginary part of the self-energy (see
F(q)=
|q+G|2 Fig. 1), a task virtually impossible within the PlP
XG
approximation57. The spectral function which can be
and choose to add 6 points to the original q mesh in obtained directly from the self-energy
order to get the new mesh. Equation (17) is then solved
6
10 0.2
00
band1
bands 2-4
V) -1V) bands 5-7
> (e-10 ω) (e0.1 band 8
ω) (A
(
Σ
<-20 band 1
e
R bands 2-4
bands 5-7
-30 0
band 8 -60 -40 -20 0 20 40 60
ω (eV)
-4-0100 -50 0 50 100 FIG. 2: Spectral function hΨkm|A(ω)|Ψkmi for the first 8
40 bandsforsiliconatk=0. Thezeroofenergyisatthecenter
ofthebandgap. ThesharppeaksareQPpoles,theirweights
band 1 correspond to thefactor Z definedin Eq. 4.
30 bands 2-4
bands 5-7
V)20 band 8
this energy step we determine an energy grid which we
e
> ( use to produce an accurate fit to WGCG′(q,±(ω−ǫn(k−
ω) 10 q))) for the different points ±(ω−ǫ (k−q)). All these
n
(
Σ high values of the parameters ensure the convergence of
<
m 0 the QP energies to within 0.05 eV.
I
Table III shows the excellent agreement of our results
-10 with twoother all-electronGWAimplementations ofthe
QP energies of silicon. From this table it seems that, at
least for Si, the overall difference between the RPA and
-20
-100 -50 0 50 100 the PlP results is small.58. Nevertheless, a discrepancy
ω (eV) of as much as 0.18 eV for the energy of L is obtained.
1v
It seems then, at least for Si, the PlP model overesti-
FIG.1: RehΨkm|Σ(ω)|Ψkmi andImhΨkm|Σ(ω)|Ψkmishown mates only slightly the differences between the energy
for the first 8 bands for silicon at k=0. The zero of energy
levels within the GWA.
is at thecenter of the band gap.
Table IV compares the calculated QP energies for
3C-SiC (also known as β-SiC), AlAs and InAs with
experimental data as well as with pseudopotential-GWA
hΨ |A(ω)|Ψ i=
km km
(PP-GWA) calculations. The band gaps are given at Γ,
|ImhΨ |Σ(ω)|Ψ i| X,andLandareunderlinedinthis table. Thesestudies
km km
[ω−ǫ (k)−RehΨ |Σ(ω)|Ψ i]2+[ImhΨ |Σ(ω)|Ψ i]2are motivated by the fact that SiC is a material of high
m km km km km
current technological interest and that InAs is predicted
is of major interest since it can be used for the in- to be metallic in the LDA, whereas AlAs is used as a
terpretation of experimental photoemission and inverse- test case. A general trend of our implementation is
photoemissionspectra. Asanexample,thespectralfunc- that the agreement with experiment as well as with
tion of silicon at the Γ point is shown in Fig. 2. The PP-GWA results is not perfect, as also found by other
sharp peaks correspond to QP excitations, while the in- implementations based on all-electron methods;15,17
coherentpartofthefunction,thespectralbackground,is reporting differences up to 0.4 eV. In particular, the
muchcomplicatedandcouldcorrespondtoplasmontype largestdifferenceoccursfortheΓ′ →X transition.59
25v 1c
excitations. We showed in a previous study15 for the case of silicon
The QPcalculations havebeen performedusing 256k that these differences are mainly traced back to differ-
points in the full BZ. The size of the dielectric matrix ences between the exchange-correlationmatrix elements
defined in Eq. (6) is 137×137 for siliconand SiC, 169× obtained by the two methods. We believe that this can
169 for AlAs, 181×181 for InAs. 200 bands were used be extended to other materials as well, since it seems to
forthesumoverconductionstatesinEq. (6)andforthe be a general feature60 of all-electron GWA calculations.
sum over m in Eq. (13). Due to the smoothness of the In fact, the difference between all-electron and PP
integrand along the imaginary axis, 11 points are found based GW calculations is not surprising since the use
sufficient to obtain well convergedquantities. An energy of pseudo-wave functions for evaluating matrix elements
step of 1.5 eV is used for the part of Eq. (15) which of a general operator produce results that may not be
corresponds to the poles of the Green’s function. Using sufficiently precise, because the pseudo-wave functions
7
TABLE III: Selected energy eigenvalues, in eV, at Γ, X and L for Si. Our results are compared with two other all-electron
implementations of the GW method and with experimental results. Data in parentheses are results when the denominator of
theGreen’s function is updated with QPenergies. Data in the last line correspond to theminimum energy gap Eg.
LDA GW approximation Expt.c
Present LAPWa Present PAW-PlPb LAPWa
Γ1v −11.97 −11.95 −11.85 (−11.89) −11.92 −12.21 -12.5±0.6
Γ′25v 0.00 0.00 0.00 (0.00) 0.00 0.00 0.00
Γ15c 2.54 2.55 3.09 (3.15) 3.16 3.30 3.40,3.05d
Γ′2c 3.23 3.17 4.05 (4.12) 4.09 4.19 4.23, 4.1d
X1v −7.82 −7.82 −7.74 (−7.78) −7.91 −8.11
X4v −2.85 −2.84 −2.90 (−2.92) −2.98 −3.03 -2.9e, -3.3±0.2f
X1c 0.61 0.65 1.01 (1.08) 1.10 1.14 1.25d
X4c 10.02 10.64 (10.72) 10.74
L′2v −9.63 −9.63 −9.57 (−9.60) −9.66 −9.92 -9.3±0.4
L1v −6.99 −6.98 −6.97 (−7.00) −7.15 −7.31 -6.7±0.2
L′3v −1.19 −1.19 −1.16 (−1.17) −1.24 −1.26 -1.2±0.2
L1c 1.44 1.43 2.05 (2.11) 2.08 2.15 2.1g,2.4±0.15
L3c 3.30 3.35 3.83 (3.90) 3.92 4.08 4.15±0.1h
Eg 0.55 0.52 0.92 (0.95) 0.97 1.01 1.17
aReference14
bReference15
cUnlessnoted,Ref35
dReference38
eReference39
fReference40
gReference41
hReference42
are constructed to reproduce the all-electron wave everonlyfewstudieshavebeenpublishedabouttheelec-
functions only in the interstitial region, i.e., outside the tronicstructures.48,49,52 Inthesematerials,hydrogenbe-
atomic spheres. This is however a good representation havesasanH−ion,leadingtopartlyionicmaterialswith
forstudyingpropertiesthatdependonlyonthebehavior a larger band gap than the studied ’sp’ semiconductors.
of the wave function in the bonding region. An error is It is, therefore, of interest to know whether the GWA is
then introduced in any PP-GWA calculation. This error capableofproducingsuchlargebandgaps. Inthisstudy
seems to fortuitously have a tendency to improve the we are only concerned with the determination of their
agreement with experiment, explaining the exceptional QP energies for the rocksalt crystallographic structure.
success of the PP-GWA. InFig. 3 wereportthe QPbandstructuresofNaHand
For InAs, the incorrect metallic behavior obtained KH along some high-symmetry directions,61 and present
within the LDA is corrected by our GWA calculation. a detailedoverviewoftheir LDA andQP energiesin Ta-
The GWA produces the true semiconducting state as ble V. In the following we detail and compare our re-
given by experiment. Since we did not account for the sults with existing experimental and calculated results:
spin-orbit coupling in our calculation, we have simply For NaH, the authors of Ref. 48 reported an LDA cal-
averaged out the spin-orbit-split experimental values to culation with an improved LMTO method in its atomic
make the comparison with our work possible. sphere approximation, however, from their band struc-
ture plotwe estimated thattheir band gapis only about
2.7 eV and is direct at the L point. This disagrees with
ourcalculation,since atthe LDA levelwe foundanindi-
rectbandgapof3.39eVfromW toL. Thiscouldbedue
B. Results for the alkali hydrides NaH and KH
to Ref. 48’s use of a different value of the lattice param-
eter of 8.90 a.u. A calculation by Kunz and Michish49
Alkali hydride materials exhibit a structural phase
based on the electronic polaron model produced a direct
transitionfromthe B1 (NaClstructure)to the B2 (CsCl
bandgapatX of1.52eV,a valuetoo lowfor theseionic
structure) type structure under hydrostatic pressure. A
materials,sothatitshouldbetakenonlyataqualitative
numberofstudieshavebeenperformedtounderstandthe level.62 However, our LDA results are in full agreement
equation of state of these materials,36,37,48,49 as well as
with results ofRef. 52, where an indirect bandgap from
the possibility ofaninsulator-metaltransition,50,51 how-
8
TABLE IV: Quasiparticle energies in eV at Γ, X and L for SiC, AlAs, and InAs. Data in parentheses are results when the
denominator of the Green’s function is updated with QP energies. Our results are compared with PP-GW method and with
experimental results (minimum band gaps are underlined).
LDA GW Expta
Present PPb Present PPb
SiC
Γ′25v →Γ15c 6.25 6.41 7.32(7.45) 7.35
Γ′25v →X1c 1.28 1.31 1.80(1.89) 2.34 2.39
Γ′25v →L1c 5.34 5.46 6.45(6.56) 6.53 4.2
AlAs
Γ′25v →Γ15c 1.94 1.77 2.72(2.79) 2.75 3.11c
Γ′25v →X1c 1.32 1.20 1.57(1.65) 2.08 2.24
Γ′25v →L1c 2.06 1.89 2.73(2.80) 2.79 2.49c;2.54d
InAs
Γ′25v →Γ15c -0.07 -0.39 0.46(0.49) 0.59 0.60
Γ′25v →X1c 1.48 1.57 (1.61) 2.10
Γ′25v →L1c 1.05 1.54 (1.58) 1.52 1.74
aUnless noted Ref 35. For InAs, data have been averaged to
account fortheneglectofspin-orbitcouplinginourcase.
bFor SiC, PP results are from Ref. 43, for AlAs from Ref. 44,
for InAs from Ref. 45 and averaged to account for the neglect of
spin-orbitsplitting.
cReference46.
dReference47.
9
TABLEV: QuasiparticleenergiesineVatΓ,X,L,W,andK
forNaHandKH.Thelastlineshowstheminimumbandgap
Eg. Data in parentheses are results when the denominator
of theGreen’s function is updatedwith QPenergies. Weare
notawareofanyexperimentalstudyconcerningtheelectronic
band structureof thesematerials.
NaH KH
LDA GW LDA GW
Γ1v -3.64 -3.39 (-3.59) -2.07 -2.02 (-2.11)
Γ1c 9.98 11.89 (12.38) 7.47 9.81 (10.43)
Γ25′c 14.98 16.71 (17.24) 7.57 10.16 (10.74)
Γ15c 15.72 17.87 (18.11) 9.33 11.98 (12.53)
X1v -0.13 -0.12 (-0.12) -0.18 -0.18 (-0.20)
X4′c 4.00 6.23 (6.66) 4.91 7.49 (8.01)
X3c 9.59 11.45 (11.96) 5.18 7.92 (8.47)
X5′c 9.95 12.15 (12.60) 8.76 11.68 (12.24)
X1c 10.84 12.95 (13.36) 9.27 12.31 (12.86)
L1v -1.03 -1.03 (-1.06) -0.15 -0.19 (-0.20)
L2′c 3.39 5.68 (6.11) 3.18 5.85 (6.35)
L1c 10.24 12.13 (12.59) 6.86 9.69 (10.23)
L3′c 12.89 15.14 (15.61) 8.46 10.49 (11.18)
L3c 15.42 17.34 (17.77) 10.98 13.57 (14.03)
W1v 0.00 0.00 (0.00) 0.00 0.00 (0.00)
W3c 5.43 7.79 (8.24) 4.96 7.60 (8.13)
W2′c 7.75 10.08 (10.55) 6.70 9.40 (9.94)
W1c 15.50 17.57 (17.98) 10.50 13.03 (13.55)
W3c 17.08 19.16 (19.57) 10.73 13.25 (13.78)
FIG. 3: Calculated LDA (doted lines) and GW (full lines)
electronic band structures of NaH and KH along some high- K1v -0.13 -0.12 (-0.12) -0.06 -0.06 (-0.07)
symmetrydirections. Thesecalculations areperformed using K3c 4.71 6.96 (7.41) 4.84 7.50 (8.02)
theparameters reported in Table II. K1c 5.84 8.11 (8.56) 4.94 7.62 (8.18)
K4c 10.57 12.85 (13.27) 8.13 10.72 (11.38)
K1c 11.02 13.09 (13.52) 8.36 10.94 (11.54)
W to L of about 3.3 eV was found. Table V shows that Eg 3.39 5.68 (6.11) 3.18 5.85 (6.35)
the energy level differences of NaH are substantially in-
creased by the use of the GWA compared with the LDA
results. In particular, the minimum band gap is 5.68 eV
within the GWA, whereas it’s only 3.39 eV within the
LDA. An other interesting point is that we found that
the most used scissors-operatorshift seems not to apply
for the computation of optical properties of NaH. This 2.71 eV at the L point.
is because the band gap is increased by 1.66 eV at the It is surprising to notice that, across the whole Bril-
Γ point, and as much as 2.34 eV for the direct transi- louin zone, the GW band gap shift of KH is larger than
tion at L point. A GWA calculation of the QP energies that of NaH despite that the LDA band gap of NaH is
in the full BZ is then required for the study the optical largerthanthatofKH.Thereasonforthispuzzlinglarge
properties of NaH. Table V shows also our LDA and QP shiftisthatthescreeningoftheCoulombpotentialinKH
results for potassium hydride. Our LDA results are in isfoundtobe lessefficientthaninNaH.Indeed,wehave
goodagreementwith the calculationin Ref. 52. We find found that the RPA static dielectric function of NaH is
thatinbothLDAandGWAcalculations,thebandgapis 3.43 much larger than that of KH which is about 2.62.
indirectfromW toL,andcouldeasilyswitchtoa direct ThehybridizationisalsolessstrongforKHthanforNaH,
gap with a small lattice parameter variation since the since the band width of hydrogen s-states is about 2.02
valence band (1s state of hydrogen) is very flat, i.e., the eV for KH (much smaller than the 3.64 eV band width
valences states at the X, L, W, and K have the same of NaH). As a consequence,the higher excited states are
energies within 0.2 eV accuracy. The previous remark lower in energy, across the Brillouin zone, by about 2 to
about the non-validity of the scissors-operator shift still 5 eV for KH than for NaH. We hope that our predictive
holds here for the same reasons. The LDA band gap in- resultswillstimulateexperimentaliststoperformphotoe-
creases by an amount of 2.29 eV at Γ, and as much as mission or optical studies of these interesting materials.
10
IV. CONCLUSION as the computation of QP lifetimes will be presented
in future work. The method is currently being applied
We have presented an implementation of the GWA to determine the excitation properties of LiH, and
using the all-electron PAW method where the screened the results will be reported elsewhere.53 Moreover,
Coulomb interaction is obtained using the RPA dy- the use of symmetry and an efficient implementation
namical dielectric function. Thus we avoided the use make us confident that we will soon be able to study
of the plasmon-pole approximation. We have applied systems with a large number of atoms per unit cell, like
it to study the QP energies of Si, SiC, AlAs, and surfaces or polymers. Furthermore, because we use a
InAs and found that a precise comparison with other mixed basis-set in our implementation we could study
availabletheoreticalandexperimentalresultsshowsthat systemswith localized’d’or’f’electronswith a reduced
sometimesthe GWAcanleadtonoticeablediscrepancies computational cost compared with methods based only
with experiments. Those discrepancies are generally not a plane wave basis-set.
so pronounced in the pseudopotential GWA calculations
using PlP models. We argued that the approximations
used in the pseudopotential method have a tendency to
fortuitously improve the agreement with experiment. Acknowledgments
We have presented detailed results for the first time
for NaH and KH alkali hydrides, and showed that Oneofus(S.L.)isparticularlygratefultoW.E.Pick-
the GWA enhanced substantially the LDA band gaps, ett since part of this work was done during a visit to
motivating further theoretical and experimental studies. the University of California Davis supported by DOE
Since our method can compute the imaginary part grant DE-FG03-01ER45876. Supercomputer time was
of the self-energy, we could then determine the QP provided by CINES (project gem1100)on the IBM SP3.
lifetimes, a task not possible using the PlP approxima- Thisresearchwassupportedinpartbythe NationalSci-
tion. Further inspection of spectral properties as well ence Foundation under Grant No. PHY99-07949.
1 P. Hohenberg and W. Kohn, Phys. Rev. 136 (1964); W. 121, Issues 9-10, 461 (2002).
Kohn and L.J Sham, Phys.Rev. 140, A1113 (1965). 18 W. Ku and A. G. Eguiluz, Phys. Rev. Lett 89, 12401
2 L. Hedin,Phys. Rev. 139, A796 (1965). (2002).
3 L. Hedin and S. Lundquist, in Solid State Physics, edited 19 H. N. Rojas, R. W. Godby, and R. J. Needs, Phys. Rev.
by H. Ehrenreich, F. Seitz, and D. Turnbull (Academic, Lett 74, 1827 (1995).
New York,1969), Vol.23, p. 1. 20 F.Bechstedt,M.Fiedler,C.Kress,R.DelSolePhys.Rev.
4 F.AryasetiawanandO.Gunnarsson,Rep.Prog.Phys.61, B 49, 7357 (1994).
237-312 (1998). 21 A. Fleszar and W. HankePhys.Rev.B 56, 10228 (1997).
5 W. G. Aulbur, L. J¨onsson, and J. W. Wilkins, ’Quasipar- 22 P.E Bl¨ochl, Phys. Rev.B 50, 17953 (1994).
ticle calculations in solids’, in Solid State Physics; edited 23 G.Kresse andD.Joubert, Phys.Rev.B. 59, 1758 (1999).
by H.Ehrenreich and F. Spaegen, vol 54. 24 N. A. W. Holzwarth, G. E. Matthews, R. B. Dunning, A.
6 G.Onida,L.Reining,andA.Rubio,Rev.Mod.Phys.74, R.Tackett, and Y. Zeng, Phys. Rev.B 55, 2005 (1997).
601 (2002). 25 N. A.W. Holzwarth, G. E. Matthews, A. R. Tackett,and
7 F. Aryasetiawan, ”Advances in Condensed Matter Sci- R.B. Dunning,Phys.Rev.B 57, 11827 (1998).
ence;”editedbyI.V.Anisimov(GordonandBreach2000). 26 S.Massidda,M.Posternak,andA.Baldereschi,Phys.Rev.
8 J-L. Li, G.-M. Rignanese, E. K. Chang, X. Blase, and S. B 48, 5058 (1993)
G. Louie, Phys.Rev. B 66, 035102 (2002). 27 S. D. Adler, Phys. Rev. 126, 413 (1962); N. Wiser, Phys.
9 P. van Gelderen, P. A. Bobbert, P. J. Kelly, G. Brocks, Rev.129,62(1963), seealsoD.L.Johnson,Phys.Rev.B
and R.Tolboom, Phys. Rev.B 66, 075104 (2002). 9, 4475 (1974).
10 A. Marini, G. Onida, and R. Del Sole, Phys. Rev. Lett. 28 M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390
88, 016403 (2002). (1986).
11 M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 29 W. G. Aulbur, PhD thesis, The Ohio State University,
(1986). 1996.
12 R. W. Godby,M. Schluter,and L. J. Sham, Phys. Rev.B 30 Y.-G.JinandK.J.Chang,Phys.Rev.B59,14841(1999).
37, 10159 (1988). 31 K.-H.LeeandK.J.Chang,Phys.Rev.B54,8285(1996).
13 M.Rohlfing,P.Kru¨gerandJ.Pollmann,Phys.Rev.B57, 32 H. J. Monkhorst and J.D. Pack, Phys. Rev. B 13, 5188
6485 (1998). (1976).
14 N. Hamada, M. Hwang and A.J. Freeman, Phys. Rev. B 33 F.GygiandA.Baldereschi,Phys.Rev.B34,4405(1986).
41, 3620 (1990). 34 B. Wenzien, G. Cappellini, and F. Bechstedt, Phys. Rev.
15 B. Arnaudand M.AlouaniPhys.Rev.B62,4464 (2000). B 51, 14701 (1995).
16 J. Furthmu¨ller, G. Cappellini, H.-Ch. Weissker, and F. 35 Numerical Data and Functional Relationships in Sci-
Bechstedt, Phys. Rev.B 66, 045110 (2002). ence and Technology, edited by K.H. Hellwege and O.
17 T. Kotani and M. Van Schilfgaarde, Solid State Comm. Madelung, Landolt-B¨ornstein, New Series, Group III,