A variational polaron self-interaction corrected total-energy functional for charge excitations in wide-band gap insulators

0.57 MB·

by  Babak Sadigh

Avariationalpolaronself-interactioncorrected
total-energyfunctionalforchargeexcitationsininsulators
Babak Sadigh,1,∗ Paul Erhart,2 and Daniel A˚berg1
1LawrenceLivermoreNationalLaboratory,PhysicalandLifeSciencesDirectorate,Livermore,California,94550,USA
2ChalmersUniversityofTechnology, DepartmentofAppliedPhysics, S-41296Gothenburg, Sweden
(Dated:July29,2015)
We conduct a detailed investigation of the polaron self-interaction (pSI) error in standard approximations
totheexchange-correlation(XC)functionalwithindensity-functionaltheory(DFT).ThepSIleadstodelocal-
izationerrorinthepolaronwavefunctionandenergy,ascalculatedfromtheKohn-Sham(KS)potentialinthe
nativechargestateofthepolaron.ThisconstitutestheoriginofthesystematicfailureofDFTtodescribepolaron
5 formationinbandinsulators.Itisshownthatthedelocalizationerrorinthesesystemsis,however,largelyabsent
1 intheKSpotentialoftheclosed-shellneutralchargestate.ThisleadstoamodificationoftheDFTtotal-energy
0 functionalthatcorrectsthepSIintheXCfunctional. TheresultingpSIC-DFTmethodconstitutesanaccurate
2 parameter-freeabinitiomethodologyforcalculatingpolaronpropertiesininsulatorsatacomputationalcostthat
isordersofmagnitudesmallerthanhybridXCfunctionals. Unlikeapproachesthatrelyonparametrizedlocal-
l
u izedpotentialssuchasDFT+U,thepSIC-DFTmethodproperlycapturesbothsiteandbond-centeredpolaron
J configurations.Thisisdemonstratedbystudyingformationandmigrationofself-trappedholesinalkalihalides
(bond-centered) as well as self-trapped electrons in an elpasolite compound (site-centered). The pSIC-DFT
7
2 approachconsistentlyreproducestheresultsobtainedbyhybridXCfunctionalsparametrizedbyDFT+G0W0
calculations. Finally,wegeneralizethepSICapproachtohybridfunctionals,andshowthatinstarkcontrastto

] conventionalhybridcalculationsofpolaronenergies,thepSIC-hybridmethodisinsensitivetotheparametriza-
i tionofthehybridXCfunctional.Onthisbasis,wefurtherrationalizethesuccessofthepSIC-DFTapproach.
c
s
-
l
r I. INTRODUCTION The present work builds on two ideas: (i) the LDA/GGA
t
m failuretopredictpolaronformationismainlyduetotheself-
interaction (SI) of the localizing carrier6,12,13 and (ii) while
. Polaronsforminmanywide-bandgapmaterialsasaresult
t LDA and GGA severely underestimate the band gaps of in-
a of charge carrier localization due to strong electron-phonon
sulators, they accurately predict band edge variations, in-
m coupling. Theirformationandtransportarecrucialtounder-
duced by configurational changes.14,15 Based on these con-
standing important quantum mechanical processes occurring
-
d in electrochemical devices such as batteries1 and scintilla- siderations, we propose a simple modification to the DFT
n tiondetectors.2Conventionalapproximationstotheexchange total-energy functional, the polaron self-interaction correc-
o tion(pSIC),thataccountsforthespuriouspolaronSIwithout
and correlation (XC) potential within density-functional the-
c explicitadditionofEX,andwithoutanimplicitlimitationto
ory (DFT), such as the local density approximation (LDA)
[
and generalized gradient approximations (GGA),3 fail qual- site-centeredlocalpotentials. TheresultingpSIC-DFTfunc-

3 tional is parameter-free, can be used for ab initio molecular-
itatively to describe polaron formation. Hybrid functionals,
v which combine exact exchange (EX)4,5 with LDA/GGA to dynamics simulations of hole as well as electron polarons in
7 arbitrarybandinsulatorswithanaccuracythatisonparwith
correct for the band gaps of insulators, are capable of pre-
3
thebesthybrid-DFTparametrizations,atacomputationalcost
1 dicting polaron formation but are computationally orders of
comparablewithstandardDFT.Itsstrengthandversatilityis
7 magnitudemoreexpensivethanDFT.
demonstrated by studying formation and migration of hole-
. Forsmallpolaronslocalizedatspecificatomicsites,which
1
polarons in the alkali halide family of compounds, and elec-
0 is the case in in many oxides,6–9 an orbital-dependent local
tron polarons in Cs LiYCl (CLYC), which is an elpasolite
4 potentialinthespiritoftheDFT+U approach10 canbeadded 2 6
compound that has received much attention in recent years
1 tolocalizetheexcesschargeandstabilizethepolaron. While
: thisenablesonetosignificantlyreducethecomputationalcost duetoitspotentialasascintillatormaterial.
v
i for calculating polaron energies, the approach is not viable Itisimportanttopointoutthatthepolaronself-interaction
X if the polaronic state involves multiple atomic sites and/or correction can be applied to any XC functional. Naturally,
r largeatomicdisplacements. Thisisunfortunatelyaverycom- whenappliedtohybridfunctionals,thepSICmethoddoesnot
a
mon occurrence, e.g. in halides, where polaron formation decreasethecomputationalcost.However,alsointhesecases,
proceedsviadimerization(bond-centeredpolarons,so-called pSICleadstoasystematicimprovementandsignificantlyre-
V -centers11) or in the case of polaron migration. Further- duces the sensitivity of the results to hybrid parametrization.
K
more, both hybrid-DFT and local potential methods depend ThislendsfurthercredibilitytothepSICmethodasageneral
onadjustableparameters. ThefractionofEXorthestrength formalismapplicabletopolaronsincrystallineanddisordered
of the local potential are usually determined a priori for a insulatorsaswellasmolecularsystems,whenevertheremoval
reference configuration and then kept fixed throughout con- or addition of an electron from/to a frontier orbital (highest
figurational changes. In summary, numerous complications occupied or lowest unoccupied quasi-particle state) drives a
havelimitedthesystematicinvestigationofpolaronsandtheir structuralchange.
transportfromfirstprinciples. This paper is organized as follows. The pSIC method is
2
Vk−center ideal iondimerI−(aV -center). Thedimerizationisaccompanied
V) 1.2 2 K
ce (e 1.0 Na I PopBtEim0ized hybrid bcoymdmisopdlaacteemtheenltasttoicfethsterasiunr.rFouignudrieng1sahtoowmsstthheaetnheerlgpytaoloancg-
ct latti 0.8 DFT abapsaerdticounlathrepaPtBhwEaXyCtofsuenlcf-titoranpapl,i2n0g(ciia)lcthuelaPteBdEf0rohmyb(ir)idD-EFXT
erfe 0.6 pSIC-DFT functional with 0.25 fraction of EX,4 and (iii) an optimized
p
o 0.4 hybrid-EX functional in which the fraction of EX (0.31) has
e t beenchosentoreproducethebandgapofNaIcalculatedwith
ativ 0.2 PBE+G W . Further computational details are listed in the
el 0 0
y r 0.0 Appendix. Figure 1 clearly shows that when a hole charge
g
er −0.2 is present in the system, the PBE XC functional predicts the
n
E undistortedlatticetobethelowestenergyconfiguration. Ac-
2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6
cordingtoPBE,theholeinthevalenceshellisthusahomo-
Iodine separation (Å)
geneouslydistributeddelocalizedcarrier. Ontheotherhand,
FIG.1. Energyasafunctionofthetwoiodineionsthatformthecore Fig. 1 also demonstrates that hybrid functionals, which in-
oftheV -centerinNaI.The“optimizedhybrid”hasanEXfraction clude some fraction of EX, allow the system to lower its en-
K
of0.31. InagreementwiththehybridfunctionalspSIC-DFTyields ergyviaalatticedistortionthatcouplestotheholechargeand
astablepolaronconfiguration. UncorrectedDFTqualitativelyfails triggers its localization. Hybrid functionals are thus capable
toreproduceself-trapping. ofstabilizingholepolaronsinNaI.FromFig.1,itisalsoev-
ident that the polaron binding energy sensitively depends on
thechoiceoftheEXfractioninthehybridfunctional.
derived in Sect. II. Section III presents a comprehensive in-
Itisimportanttomentionherethatintheabsenceofcharge
vestigationoftheaccuracyofthepSICmethodincomparison
excitations, DFT predicts accurate ground state (GS) ener-
with the hybrid technique. In particular, formation and mi-
gies as well as phonon spectra for NaI. Furthermore, Fig. 2
gration of hole polarons in alkali halides, and formation of
(left panel) shows that in the neutral charge state, PBE and
electron polarons in CLYC are studied. Section IV summa-
G W quasiparticle (QP) energies for a supercell containing
rizesanddiscussesthefindingsofthepresentwork. Finally, 0 0
atomic displacements of a V -center, are in good agreement
computationaldetailsaregivenintheAppendix. K
witheachotherapartfromarigidrelativeshiftofthevalence
andtheconductionbands,theso-calledscissorshiftthatcom-
pensatesfortheabsenceofderivativediscontinuityinthePBE
II. THEPSICMETHOD
functional.AmoredetailedcomparisonbetweentheQPspec-
trafromPBE,optimizedhybrid,andG W isgiveninFig.3.
0 0
In this section, we motivate and develop the pSIC method The close agreement of the QP energies implies that for the
byrevisitingthefailureofstandardDFTfunctionalsforhole neutralchargestate, thePBE-KSpotentialisagoodapprox-
self-trappinginthealkalihalidecompoundNaI.3Itisanionic imation to the optimized effective potential (OEP) obtained
compound in the rocksalt crystal structure with an equilib- from the exact XC functional.21 As a result, for this charge
rium lattice constant of 6.48A˚.16 The crystal consists of two statePBEaccuratelypredictsboththeGSenergyoftheV -
K
interpenetrating fcc sublattices, each consisting of either Na center, as well as the QP energy of the hole state relative to
or I ions only. The ionic bonding in NaI leads to an insu- thevalencebandmaximum(VBM).Incontrast,Fig.2(right
lating ground-state electronic structure with a band gap of panel) shows that in the charged state, the unoccupied hole
5.9eV.17–19 Thevalencebandisdominatedbyhybridizedio- level virtually overlaps with the VBM, which implies a de-
dinep-orbitals,whiletheconductionbandconsistsmainlyof localized hole character. This explains why PBE predicts no
sodium s-orbitals. When a hole is introduced in the valence bindingfortheV -centerinFig.1.
K
band,itmediatescovalentbondingbetweennearestneighbor
Beforeembarkingonaformaldiscussion,itisbeneficialto
iodineionsleadingtodimerization;thisimpliesholelocaliza-
inspect the failure of DFT for the charged V -center in NaI
tion and thus the formation of a polaron.11 In the following K
fromtheperspectiveoftheSIerrorofthelocalizedhole. This
section, the pathway to polaron formation in NaI is studied
is illustrated in the center panel of Figure 2, where the de-
from various angles using DFT as well as hybrid XC func-
viationoftheDFTenergyfromlinearityasafunctionofthe
tionals.
fractionalholechargeinthesystemisshownfor(i)theperfect
crystal and (ii) the dimerized V -center configuration. It is
K
apparentthatthedelocalizedholecarrierintheperfectcrystal
A. Basicconsiderations (withnegligibleSI)exhibitsverylittledeviationfromlinear-
itywhilesignificantnon-linearityisobservedintheenergetics
WhenaholeisintroducedinthevalencebandofNaI,itin- ofthelocalizedpolaronintheV -centerconfiguration. Inac-
K
ducesasubstantiallatticedistortion,causingapairofnearest- cordwithrecentliterature,6–9,12,13,22–24weascribethisnonlin-
neighbor I− ions to move along (cid:104)110(cid:105) toward each other. earitytotheDFTSIerrorforthelocalizedpolaron. Notethat
Their separation is reduced from 4.5A˚ in the perfect crystal thesolidbluelineinFig.2(centerpanel)curvesinsuchaway
to about 3.3A˚, effectively leading to the formation of an an- as to raise the energy of the charged polaronic configuration
3
Neutral (Δ=0) Charged (Δ=+1)
0
V) ΔEDFT[0]
e
nergy ( 4 y (eV) -1 coVnKf-igceunratetiron
e g
cle 2 ner ΔEDFT[+1]
arti 0 al e
si-p Tot -2 ideal lattice
a
Qu μD±FT
ΔΠ [+1]
DFT
0 0.5 1
DFT G0W0 Hole charge ΔEf (VK-center) DFT G0W0
FIG.2. Thecenterpanelshowsthevariationofthetotalenergyoftheperfectlattice(dottedredline)andtheV -center(solidblueline)
K
configurationswithfractionalholecharge.TheSIcanbecorrectedbyremovingthecurvatureofthesolidbluelineasillustratedbythedashed
greenline.Theoffsetbetweentheuncorrected(solidblue)andSIcorrectedinthelimit∆=+1thencorrespondstothetermΠ [∆=+1]
DFT
in Eq. (1). The left and right panels show side-by-side views of quasi-particle energies calculated within DFT and G W for neutral and
0 0
charged,respectively,64-atomsupercellsofNaIcontainingaV -center.
K
8 charge state by Ne. Let EDFT(Ne ±∆;R) denote the DFT
(a) DFT (b) opt. hybrid (c) G0W0 energy of this system as a function of the fraction ∆ of ex-
V) 6 cess electrons/holes. The ion positions are specified by the
e
y ( 3Np-dimensional vector R, where Np is the number of ions
g 4
er inthesystem. Webreakdownthecontributionstotheenergy
n
e atarbitrary∆asfollows
e 2
cl
arti 0 EDFT[Ne±∆;R]=EDFT[Ne;R]±∆·µ±DFT[R]
p
si− +ΠDFT[±∆;R], (1)
ua −2
Q whereµ± istheelectronchemicalpotential,definedas
DFT
−4
∂E
1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 µ± [R]= lim DFT[N ±η;R], (2)
k−point index DFT η→0 ∂N e
andΠ quantifiesthedeviationofE fromlinearityand
FIG.3. EnergylevelsfortheV -centerconfigurationcalculatedus- DFT DFT
K
ingDFT,optimizedhybrid,andG W (basedonPBE0calculations) thuscanbeusedasameasureofthepolaronSI.Figure2(cen-
0 0
forNaI.Redandbluelinesrepresentoccupiedandunoccupiedlev- ter panel) shows that the VK-center configuration becomes
els, respectively. Thepolaroniclevelisshownbyorangelinesand stabilizedfor∆=1ifthetermΠDFT(thepolaronSI)iselim-
opencircles. inatedfromtheaboveequation.
Atfirstsight,thisseemstocounterthewell-knownfactthat
theelectronchemicalpotentialsµ± ,calculatedwithinDFT,
DFT
relativetothatoftheideallattice,andthuspreventsbinding. areinaccuratepredictorsofthehighestoccupied(HOMO)and
The effect of the SI error on the stability of polarons within lowest unoccupied molecular (LUMO) levels in molecules
DFThasalsobeendiscussedinRef.12. andofbandgapsinsolids.ThecurvaturesofLDA/GGAfunc-
tionals resulting from their SI are, however, such that they
correct the values of the differential quantities µ± for the
DFT
B. ThepSIC-DFTtotalenergyfunctional HOMO/LUMO levels in small molecules. This constitutes
thebasisfortheso-calleddelta-SCFmethod,whichhasbeen
Wenowformalizetheabovediscussionwiththegoaltode- quite successful in molecular applications. It is however, a
rive a robust method for structural relaxations as well as ab precarious practice to correct for the lack of derivative dis-
initiomoleculardynamicssimulationsofpolaronsingeneral continuityinLDA/GGAXCfunctionals25,26bytheirSIerror.
insulating systems with a computational efficiency compara- In particular, Π [±∆;R] can have a spurious dependence
DFT
bletoLDA/GGA,whichdoesnotrelyonaprioriknowledge ontheionicconfigurationRandthusleadtowrongstructural
orguesswork. Consideranextendedmany-electronsystemin predictions. This can most dramatically be demonstrated in
a periodic supercell containing N ions, at a reference elec- insulatingsolidssuchasNaI,wheretheperfectcrystalstruc-
p
tronicconfiguration,say,theneutralgroundstate. Inthefol- ture has negligible SI while Π is quite substantial for the
DFT
lowing, we denote the number of electrons in this reference self-trappedpolaron. Thestrongconfigurationdependenceof
4
ΠDFTleadstothefailuretopredictpolaronformationinthese 2
systems. A systematic correction should rather correct µ± V) NaF pSIC
by adding the contribution from the derivative discontinuDiFtyT y (e NaI pSIC
of the exact exchange-correlation functional, which is miss- erg 1 NaF conventional
n NaI conventional
e
ing in LDA/GGA. Hence we propose the following polaron on 0
SIcorrected(pSIC)totalenergyfunctionalforasystemwith ati
m
anaddedelectron/hole: n for −1
o
Ep±SIC[R]=EDFT[Ne;R]±µ±DFT[R]±∆±XC[R], (3) Polar −2 αPBE0 αoNpatI αoNpatF
0.0 0.1 0.2 0.3 0.4 0.5
where ∆± is the missing contribution of the derivative dis-
XC Fraction of exact exchange
continuityoftheexchange-correlationfunctionaltothechem-
icalpotentials.27Itcanbeformallyexpressedas FIG.4. Variationof thecalculatedpolaron formationenergiesof
NaIandNaFwiththemixingparameterofthehybridXCfunctional.
∆± [R]= lim ∂Eexact[N ±η;R]−εexact-OEP[N ;R], Notethatfortheconventionalhybrid-XCfunctionalsthepolaronfor-
XC η→0 ∂N e ± e mationenergyisstronglydependentonthemixingparameter. This
(4) situationisalleviatedalmostentirelybyusingthepSIC-DFT/pSIC-
hybrid approach. The deviation between conventional hybrid and
where εexact-OEP and εexact-OEP denote the lowest unoccupied pSIC-hybridissmallestfortheoptimizedhybrids,i.e. αm,assug-
+ − gestedbyEq.(9).
andthehighestoccupiedeigenvaluesoftheexactOEPpoten-
tial in the reference charge state, respectively. Note that for
the E± functional in Eq. (3) to be accurate, it is required the total-energy functional Eq. (3). We thus arrive at the ex-
pSIC
that pression
EDFT[Ne;R]≈Eexact[Ne;R]+constant (5) ∆Ep±SIC-DFT[R]=∆EDFT[Ne;R]±∆µ±DFT[R]. (7)
µ± [R]≈εexact-OEP[N ;R]+constant. (6) Thenotation∆E and∆µisusedabovetohighlightthefact
DFT ± e
thatweareonlyinterestedinrelativeenergies.
Theaboveisguaranteed,ifforthereferencechargestate,the
KSpotentialobtainedfromapproximateXCfunctionalsused
in practice closely resembles the exact OEP. We have shown C. GeneralizationtopSIC-hybridfunctionals
in Sect. IIA that this is indeed the case for the closed-shell
neutral charge state of NaI. In fact, it is reasonable to argue WhileRef.14justifiestheassumptionofR-independence
thatfordefect-freebandinsulatorsintheclosed-shellneutral of ∆± , it also suggests that the explicit R-dependence of
XC
chargestate, commonapproximateXCfunctionalsarelikely ∆± istolargeextentdeterminedbythestaticCOHSEXap-
XC
toproducereasonablyaccurateKSpotentials. proximation, and thus little is gained by employing full dy-
The derivative discontinuity ∆± , can be calculated accu- namic screening within the GW approximation. A reason-
XC
rately via many-body perturbation theory, e.g., in the GW ablewaytoincorporatecontributionsfromstaticscreeningof
approximation.28 In this formalism, ∆± correspond to the exchangeinthetotal-energyfunctionalisthroughthehybrid
XC
so-called scissor shifts, which constitute the difference be- technique. Of course the fraction of EX in the hybrid func-
tweentheelectronaddition/removalenergiescalculatedusing tional, or in other words the magnitude of static screening,
DFT+G0W0andtheirrespectiveDFTvalues. Furthermore,it needstobedeterminedapriori.
hasbeenshownbyBaldereschi,GygiandFiorentini14,29 that Equation (7) describes a general formalism for correcting
∆± isdominatedbythelong-rangecomponentsofthestatic polaron self-interaction in any XC functional. Up to now,
XC
Coulomb-hole plus screened exchange (COHSEX) contribu- our discussion has been primarily concerned with calcula-
tions. Therefore∆± variesweaklywithlocalizedlatticedis- tionsbasedonsemilocalXCfunctionalssuchasPBE,butev-
XC
tortions such as those typically associated with polarons, so erything said still holds if one were to correct a hybrid-XC
longastheDFTproducesaccurateKSpotentialsfortheionic functional for polaron self-interaction. Bear in mind that for
configurationsofinterest. Hence,wecansimplifyEq.(3)sig- orbital-dependentfunctionals,suchashybrid-XCfunctionals,
nificantlybyneglectingthedependenceof∆± onR. Next, thechemicalpotentialsµ± ,correspondtotheeigenvaluesof
XC hyb
weobservethatwithregardtodefectsingeneralandpolarons ageneralizedKSHamiltonianwithanon-localself-consistent
in particular, one is interested in energy differences. For ex- potential. Hence the XC derivative discontinuity ∆± , is al-
XC
ample, for polarons in insulators, all energies are calculated readyincorporatedinthechemicalpotentials,22 providedthe
relativetotheundistortedcrystal,inwhichtheexcesscharge hybridfunctionaliscorrectlyparametrized.
is delocalized. The detailed expressions for calculating po- For a class of hybrid functionals, where PBE is combined
laron binding energies can be found in Sect. 2 of the Ap- with a fraction 0 ≤ α ≤ 1, of EX, Equation (7) takes the
pendix. Since we are only interested in energy differences, forms
it is possible to move the zero of energy in such a way as to
∆E± [R;α]=∆E [N ;R;α]±∆µ± [R;α], (8)
remove the contribution from the derivative discontinuity in pSIC-hyb hyb e hyb
5
The energy of a polaron in configuration R relative to the bydeducingµ± viareplacingEq.(2)byafinitedifference
DFT
undistorted crystal R , can be estimated in two ways: ei- formulasuchas
0
ther through the conventional approach using the functional
∆funEchtyibo[nNael±∆1E;±R;α],[oRr;thαr]o.ugThhteheteprmSIiCnoalpopgryoaccohnvuesnintigonthael µ±DFT =21δ(cid:0)4EDFT[Ne±2δ]−EDFT[Ne]−EDFT[Ne±δ](cid:1)
pSIC-hyb (10)
versuspSICapproachwillbeusedthroughouttheremainder
ofthispaper.ThedetailedexpressionscanbefoundinSect. 2
whereδ (cid:28) 1andthedependenceonRhasbeenomittedfor
oftheAppendix.
brevity. This relation, which is correct to the third order in
Thedifferencebetweentheconventionalandthehybridap-
δ, in conjunction with Eq. (7) allows for the atomic forces
proach is illustrated in Fig. 4. It shows polaron energies as
to be calculated from a linear combination of three separate
a function of the hybrid mixing parameter α, for two alkali
Hellmann-Feynmanforcesasfollows:
halidecompoundsNaIandNaF,calculatedusingtheconven-
tional and the pSIC methods. Note that the polaron ener- ∂E± (cid:18) 1 (cid:19)∂E
giesvarysignificantlywithα,whencalculatedusingthecon- pSIC = 1∓ DFT[N ]±
∂R 2δ ∂R e
ventionalenergyfunctional∆E [±1;R;α ],whereasthey
hyb m (cid:18) (cid:19)
depend very weakly on α, when calculated using the pSIC 1 4∂EDFT[N ±2δ]− ∂EDFT[N ±δ] .
functional ∆E± [R;α]. The insensitivity of the pSIC- 2δ ∂R e ∂R e
pSIC-hyb
hybridfunctionaltoparametrizationisthemainadvantageof (11)
the pSIC method. It allows for accurate predictions of po-
laronenergieswithαsettozero, i.e. pSIC-DFT.Inthisway Theaccuracyofthisexpressioncanbeestablishedbycompar-
thecomputationalcostofthecalculationscanbereducedsig- ing the chemical potential calculated via Eq. (10) to the one
nificantly compared to hybrid functionals, and brought to be obtainedfromtheeigenvaluespectrum. Thenumericalstabil-
onparwithsimpleDFTcalculationswithaminorimpacton ityoftheatomicforcesfromEq.(11)mainlydependsonthe
accuracy. sizeofthefinite-differenceincrementδ. Wehavefoundthat
Figure4showsthatforacertainmixingparameterα the avalueofδ = 0.025issufficienttoobtainaccurateenergies
m
conventional and the pSIC methods are equivalent. This can andconvergeforcestobelow1meV/A˚.
beformallyexpressedasfollows: We have implemented the pSIC method in the Vienna ab
initiosimulationpackage. Ateachionicstep,thewavefunc-
∆E± [R;α ]=∆E [N ±1;R;α ]. (9) tions of three electronic configurations are converged inde-
pSIC-hyb m hyb e m
pendently,representingtheelectronicstatesN ,N +δ,and
e e
Thisimpliesthatatα ,thetotalenergyasafunctionoffrac- N +2δforelectronpolarons,andN ,N −δ,andN −2δfor
m e e e e
tional hole-polaron charge has negligible curvature and the holepolarons. Theenergiesandforcesarecollectedfromall
hybridXCisthuspolaronself-interactionfree;9 inthissense, replicas and combined according to Eqs.( 10) and (11). It is
α canbeconsideredtheoptimalmixingparameter. Wewill apparentthatinthisparallelimplementation,apSIC-DFTcal-
m
seeinSect.IIIthatforthematerialsstudiedinthiswork,the culationrequiresthreetimesasmuchwalltimeasastandard
choiceofαthatreproducestheDFT+G W bandgaps,nearly DFTcalculation.
0 0
coincideswithα .
m
III. RESULTS
D. Numericalimplementationofenergyandforces
A. V -centerformationinalkalihalides
K
WenowdiscussthenumericalimplementationofthepSIC
energyfunctional. Equations(7)and(8)comprisetwoterms, LetusstartthissectionbyrecapitulatingthepSIC-DFTre-
namely (i) the total energy in the reference charge state and sultsfortheV -centerinNaI.ItisapparentfromFig.1,that
K
(ii)thechemicalpotentialinthereferencechargestate. Byin- in contrast to conventional PBE calculations, the new pSIC-
vokingKoopman’stheorem,thelatterquantitiescanbecalcu- PBEfunctionalpredictsastableV -centerconfiguration,the
K
latedfromtheeigenvaluespectrumoftheKSHamiltonianin geometryandenergeticsofwhichareinexcellentagreement
thecaseofDFTcalculations,andthegeneralizedKSHamil- withtheoptimizedhybridfunctional. Thebondlengthofthe
tonianwithanon-localself-consistentpotentialinthecaseof I− dimer that forms the core of the V -center is calculated
2 K
hybrid calculations. In particular, for electron polarons, µ+ to be 3.38A˚ (3.23A˚) when using the pSIC-PBE (optimized
correpondstotheconductionbandminimum(CBM)andfor hybrid) scheme while the energy gain due to self-trapping is
holepolarons,µ−correspondstotheVBM. foundtobe0.267eV(0.345eV).
Structuralrelaxationsandmoleculardynamicssimulations In the following, we further evaluate the accuracy of the
require the calculation of atomic forces. Equations (7) and pSIC-PBEmethodbyanextensivestudyofitspredictionsfor
(8)are,however,notdirectlysuitableforevaluationofatomic theV -centersinalkalihalides. Thesecompoundsconstitute
K
forces via the Hellman-Feynman theorem, since µ± is not adiverseclassofwide-gapinsulatorswithbandgapsranging
DFT
avariationalquantity. Explicitcalculationsofwavefunction from 5 to 14eV and lattice constants from 4.0 to 7.3A˚, see
derivatives with respect to ionic coordinates can be avoided Table II. We compare the results of pSIC-PBE calculations
6
Å) 4.0 (a) energiesonthehybridparametrizationisfurtherillustratedfor
e ( Lithium halides Sodium halides NaI and NaF in Fig. 4, highlighting a major disadvantage of
c theconventionalapproach.
an 3.2
st
di
on 2.4 In contrast, the pSIC functional is quite insensitive to the
e i fraction of EX chosen for the hybrid calculations. To illus-
d
ali tratethis,wehaveperformedpSIC-hybridcalculationsofthe
H 1.6
polaronbindingenergiesusingEq.(8). Foreachcompound,
V) 0.0 (b) calculations were carried out in the relaxed ionic configura-
e tion obtained from a conventional optimized hybrid calcula-
y (
erg −0.6 tion,usingpositivelychargedsupercells. Theresultsarecom-
n piled in Fig. 5(b), where for every compound, pSIC-hybrid
n e pSIC−opt. hybrid calculationsareshownsidebysidewithconventionalhybrid
matio −1.2 oppStI.C h−yPbBridE results. Figure5(b)showsnearlyperfectagreementbetween
or PBE0 the pSIC and the conventional approaches for the optimized
F −1.8 hybridparametrization. Thissuggeststhatforthecompounds
F Cl Br I F Cl Br I
inthisstudy,thehybridparametrizationschosentoreproduce
PBE-G W bandgaps, yieldnearlypolaron-self-interaction-
FIG.5. Comparisonof(a)halogen–halogenseparationatthecore 0 0
free XC functionals, for which Eq. (9) suggests equality of
oftheV -centerand(b)V -centerbindingenergiesbetweenpSIC-
K K
the pSIC and the conventional approaches. However, a far
PBEandhybrid-DFTcalculations.
more important conclusion can be reached when recollect-
ing that the pSIC-DFT calculations of the polaron energies
shown in Fig. 5(b), are already within 10-20% of the opti-
with conventional hybrid-XC calculations using two distinct
mizedpSIC-hybridestimates: Whilethepolaronbindingen-
parametrizations: (i)thePBE0functional,4whichusesafrac-
ergiesarehighlysensitivetohybridparametrizationwithinthe
tionof0.25ofEX,and(ii)asetofoptimizedhybridfunction-
conventionalapproach,thisvariationisanorderofmagnitude
als,whichhavebeenparametrizedtoreproducePBE+G W
0 0 smaller in the pSIC approach. This is clearly illustrated for
calculations. For these compounds, the fraction of EX used
NaI and NaF in Fig. 4, where the polaron binding energies
intheoptimizedhybridcalculationsvaryfrom0.25inLiIto
ascalculatedwithinthepSICandtheconventionalhybridap-
0.41inNaF.FurtherdetailsregardingtheG W calculations
0 0 proaches are shown in the same graph. Note that again the
areprovidedinSect. 3oftheAppendix.
two curves cross at α , see Eq. (9), which is the optimized
m
TheresultsofpSIC-PBEandhybridcalculationsarecom-
fractionofEXinthehybridfunctional.
piledinFig.5,whichshows(a)thehalogen–halogensepara-
tionatthecoreoftheV -centerand(b)theV -centerbinding
K K
energies. Withregardstothehalogendimerseparations,itis Finally, it is in order to discuss the sources of error in the
evidentthatthepSIC-PBEmethodyieldsexcellentagreement pSIC-hybridfunctionalEq.(8),whichleadtothevariationof
with hybrid calculations with a typical deviation of less than ∆E± [R,α]withtheparameterα. Forthispurpose,we
pSIC-hyb
0.15A˚ ((cid:46) 5%). focusonNaF,forwhichpSIC-PBEhasthelargestdiscrepancy
with optimized hybrid, of any of the compounds considered
Figure5(b)showsthatthepolaronbindingenergiescalcu-
inthispaper. ForNaF,pSIC-PBEpredictsapolaronbinding
latedwithinpSIC-PBEagreewellwiththeresultsfromopti-
energyof0.92eVversus1.19eVforoptimizedhybridinthe
mizedhybridcalculations. Overall,pSIC-PBEyieldspolaron
conventional approach, and 1.24 eV for optimized hybrid in
binding energies that are 10 to 20% smaller than those from
thepSICapproach.
optimizedhybrids. Ontheotherhand,formostmaterialscon-
sideredhere,especiallythosewithlargerbandgaps,thePBE0
functional, in the conventional mode of calculation, predicts Therelativelylargeerror(inexcessof20%)ofthepolaron
VK-centerbindingenergiesthataremuchsmallerthaneither bindingenergyinNaF,whencalculatedwithinpSIC-PBEcan
pSIC-PBE or optimized hybrid calculations. This is consis- berelatedtotheindividualerrorofeachterminEq.(8). For
tent with a systematic underestimation of the band gaps of thecaseofholepolaronsinNaF,Eq.(8)consistsoftwocon-
these wide-gap materials by PBE0, see Table II. As already tributions: (i)neutralgroundstateenergy,and(ii)theenergy
discussed,furtherreducingthefractionofEXfrom0.25to0, of the VBM. Both quantities are calculated as the difference
i.e. PBE,leadstoanevenlargerunderestimationofbandgaps betweenthepolaronconfigurationandtheperfectlattice. The
and,moredramatically,thefailuretopredictholelocalization first term is calculated within PBE (optimized hybrid) to be
andpolaronformationinallalkalihalides. −2.99 (−3.06), while the second term is 3.85 (4.30). The
In summary, the potential energy landscapes of polarons, relativeground-stateenergiesdifferbylessthan3%between
calculatedwithinconventionalhybridschemesarehighlyde- PBE and optimized hybrid, while the error in the relative
pendent on the particular parameterization of the hybrid XC VBM levels from PBE is no more than 10%. The relatively
functional employed, i.e. the fraction of EX included in the largeerrorinexcessof20%, inthefinalpolaronbindingen-
calculations. Thepronounceddependenceofpolaronbinding ergyiscausedbytheoppositesignsofthetwocontributions.
7
TABLEI. ActivationbarriersforV -centermigrationinNaIinunits (a) CLYC: pSIC+DFT (b) CLYC: PBE0
K
ofeVfromcalculationandexperiment.30 8
V) 6
Rotationangle Optimizedhybrid pSIC-PBE Experiment e
60◦ 0.23 0.20 0.2 gy ( 4 6
90◦ 0.28 0.24 >60◦ ner 4
e 2
120◦ 0.23 0.20 0.2 e
180◦ 0.22 0.19 articl 0 2
−p 0
si −2
a
u −2
Q −4
ideal polaron polaron ideal
  FIG.7. Energylevelsfortheidealandpolaronicconfigurationcal-
culatedusing(a)DFTand(b)PBE0forCLYC.Redandbluelines
represent occupied and unoccupied levels, respectively. The pola-
roniclevelisshownbyorangelinesandopencircles.


C. ElectronpolaronsinCLYC


Lastly,wedemonstratetheaccuracyofthepSICmethodfor
electron polarons. For this purpose, we consider Cs LiYCl
2 6
(CLYC), which is one of the most thoroughly investigated
FIG. 6. Illustration of the electron polaron in CLYC as obtained
compounds in the elpasolite family due to its potential as a
frompSICcalculations. Thearrowsindicatetherelaxationpattern.
neutron detector. Electron polarons have been found in this
Theyellowisosurfaceiscenteredatanyttriumatom.
systembothexperimentally34 andtheoretically35 usingPBE0
calculations. Itscrystalstructure(double-perovskite)hascu-
bicsymmetrywithtenatomsintheprimitivecell. Here,elec-
tron polaron calculations were performed using 80-atom su-
B. V -centermigrationinNaI percells with a 2×2×2 k-point sampling. The calculated
K
bandgapwithinPBE0is7.1eV,incloseagreementwithex-
periment(7.5eV).36 HencePBE0andoptimizedhybridyield
Thanks to its computational efficiency, the pSIC-DFT nearly identical results for this compound. By comparison,
methodprovidesanunprecedentedpotentialforstudyingpo- DFT-PBEyieldsamuchsmallerbandgapof5.0eV.Wefind
laron dynamics. To demonstrate this aspect, we also bench- inaccordwithRef.35thattheelectronpolaronlocalizesona
markcalculatedpolaronmigrationbarrierswithmeasureddif- Y-site and the polaron level has strong d-character as shown
fusivitiesinNaI.Bysymmetry,themigrationofanI− dimer inFig.6. BothpSIC-PBEandPBE0predictastretchingofY-
to a nearest-neighbor site in the rocksalt lattice can2involve Cl bonds by about 0.1A˚. The binding energy of the electron
rotations through 60◦, 90◦, 120◦ or 180◦, the last one corre- polaronis0.25eV(0.32eV)accordingtopSIC-PBE(PBE0).
sponding to a pure translation. The polaron diffusivity in al- ThePBE0valueforthepolaronbindingenergyincludesasig-
kali halides can be accurately measured in laser pump-probe nificant contribution from image charge corrections of about
experiments.30Tothisend,V -centersarecreatedandaligned 0.1eVduetothesmallsizeofthesupercells(80atoms)that
K
alongaparticulardirection(“polarized”)usingapumplaser. was computationally affordable at this time. (Note that the
Subsequently,thegraduallossofpolarizationovertimeisfol- hithertomostrecentcalculationfortheelectronpolaroninthis
lowedusingaprobelaser,fromwhichrotationalbarrierscan systemwasperformedina40-atomsupercell35). Infact,per-
beinferred. forming pSIC-PBE0 calculation of the polaron energy in the
ionic confguration obtained from the conventional PBE0 ap-
Experimentscannotdistinguishbetween60◦and120◦rota- proach, yields a polaron binding energy of about 0.28 eV, in
tionsandareinsensitivetomigrationviatranslation. InNaIat verygoodagreementwiththepSIC-PBEresult. Thisdemon-
50K, one observes almost exclusively 60◦/120◦ rotations,30 stratesthatpSICisalsocapableofpredictingelectronpolaron
indicative of the fact that their activation barriers are notice- geometriesandbindingenergiesincloseagreementwithhy-
ablysmallerthanfor90◦rotations.Itwassimilarlyshownthat bridcalculations.
60◦/120◦ rotations dominate in RbI.31,32 Table I shows acti- In order to illustrate that this agreement is not fortuitous,
vationbarrierscalculatedwithinthepSIC-PBEandoptimized we show in Fig. 7 the quasi-particle spectra of CLYC cal-
hybrid approaches using the climbing image nudged elastic culated within DFT, as well as PBE0 for the ideal and po-
bandmethod.33Bothmethodsagreeverywellwitheachother laron configurations in their charge-neutral states. Note that
andwithexperimentaldata. the scissor shift is independent of configuration and that the
8
calculated position of the localized electron level relative to andlittlecompromiseonaccuracywillresultfromdrastically
the conduction-band minimum within DFT agrees well with simplifyingthecalculationsbysettingα = 0. Onthisbasis,
PBE0. theparameter-freepSIC-PBEmethodemerges.
The pSIC formalism can be applied to more general sys-
tems than has been presented in this paper. For the sake
IV. CONCLUSIONSANDOUTLOOK of clarity, we have intentionally limited the present scope to
intrinsic polaron self-trapping in crystalline insulators. We
have shown that in these systems, the neutral charge state
In this paper we have introduced a self-interaction correc-
isdescribedaccuratelywithinPBE(ormoregenerallyDFT).
tion technique, the so-called pSIC method, that corrects the
Hence this reference state has been used to calculate the en-
DFTfailuretopredictstablepolaronsininsulators.ThepSIC-
ergyofchargeexcitationsinthepresenceoflatticedistortions
DFT method is parameter-free and easy to implement in ex-
usingEq.(9). Moregenerally,thepSICformalismcanbeap-
isting electronic structure codes. We have validated the new
pliedtobothdefectsininsulatingsolids,aswellasmolecular
methodology by benchmarking pSIC-PBE results to hybrid-
systemswheneverstructuralrearrangementsaredrivenbyad-
DFT functionals for both electron and hole polaron forma-
dition/removal of electron or holes in frontier orbitals. The
tion and transport in a number of materials. Thereby, we
pSICmethodallowsanysupercellcalculationcontainingN
have shown that (i) polaron trapping energies depend sensi- e
electrons, whetherperiodicornot, tobeperformedinanyof
tively on the parametrization of the hybrid functionals, (ii)
thethreechargestatesN ,N +1,orN −1.Wehavedemon-
the parameter-free pSIC-PBE method exhibits similar accu- e e e
stratedinthispaperthatinthecaseofperfectcrystallineinsu-
racytohybridfunctionalsthatarelaboriouslyparametrizedto
lators,thereislargevariabilitywithchargestatewithrespect
fitGW bandgapsonacasebycasebasis,and(iii)isanorder
to the quality of the KS potentials produced from PBE. By
ofmagnitudefaster.Thispresentsabreakthroughinmodeling
choosing to perform the calculations in the charge state with
thestructureanddynamicsofpolaronsininsulatorsfromfirst
thehighestqualityKSpotential,thepSICmethodenablesdra-
principles.
maticenhancementsoftheaccuracyofDFTXCfunctionals,
Fromafundamentalphysicsperspective,thepSICmethod
for prediction of potential energy landscapes of polarons in
emerges as a result of the decoupling of the two main prob-
thesesystems.
lematic features of DFT for band insulators: (i) severe un-
The subject of defects in insulators and semiconductors is
derestimationofbandgaps, and(ii)absenceoflocalizedpo-
vastandcomplex, andofhugeinteresttomanyapplications.
larons. The counterintuitive aspect of such a decoupling is
Recently,theeffectofGW correctionsontheformationener-
thatitiswell-knownthatpolaronstabilityisgreatlyenhanced
giesofchargeddefectsinsemiconductorswasexamined.37,38
forlarge-gapsystems. Hencesignificantcorrelationbetween
In these works, structural relaxations were performed using
the band gap of an insulator and its polaron binding energy
standard DFT, only to keep the computational costs down to
should be expected. Based on this, it is then logical to ar-
a manageable level. It will be the subject of future to estab-
gue that the lack of stable polarons in insulators within the
lish, if the pSIC method can also be used to obtain accurate
DFT must stem from its band gap error. In fact, the two er-
defect configurations and energies at a reasonable computa-
rors are unrelated as shown by the following arguments: (i)
tionalcost.
fortheexactXCfunctional,theelectronaddition/removalen-
ergiescanbeextractedfromdifferentialquantities,i.e. chem-
icalpotentials,(ii)formostbandinsulators,thereisacharge
ACKNOWLEDGMENTS
state, usually the closed-shell neutral charge state, for which
DFTcanpredictaccurateground-stateenergiesaswellasva-
We thank Michael Surh at LLNL for very helpful discus-
lenceandconductionbandstructures,and(iii)whileLDAand
sions. Lawrence Livermore National Laboratory is operated
GGAseverelyunderestimatethebandgapsofinsulators,they
byLawrenceLivermoreNationalSecurity,LLC,fortheU.S.
accurately predict band edge variations with configurational
DOE-NNSA under Contract DE-AC52-07NA27344. Fund-
changes. Notingthatpolaronbindingdoessensitivelydepend
ing for this work was received from the NA-22 agency. P.E.
on the variation of the ground-state energy and band edges
acknowledges funding from the Knut and Alice Wallenberg
with configuration rather than the absolute value of the scis-
Foundation and the Area of Advance – Materials Science at
sorshift,itfollowsthatpolaronbindingisindependentofthe
Chalmers. Computer time allocations by the Swedish Na-
absolutevalueofthebandgaperror.
tional Infrastructure for Computing at NSC (Linko¨ping) and
InordertorationalizethesuccessofthepSICapproach,we
C3SE(Gothenburg)aregratefullyacknowledged.
haveappliedittocorrectforthepSIinhybridXCfunctionals.
Inthisprocess,wehaveshownthatinbandinsulators,hybrid
functionalstunedtoreproducetheDFT+G W bandgapsare
0 0
Appendix:Computationaldetails
practicallypolaronself-interactionerrorfree,andthusconsti-
tutetheoptimalparametrization. Ourkeyfindinghowever,is
that calculations of polaron energies within the conventional 1. Total-energycalculations
approach, i.e. using charged supercells, are very sensitive to
thefractionαofEXincludedinthehybridfunctional. Onthe Calculations were carried out using the projector aug-
otherhand,thepSIC-hybridmethodisalmostinsensitivetoα, mentedwavemethod39 asimplementedintheViennaabini-
9
tio simulation package,40 using the supplied standard plane localfieldeffectsinε forthehybridfunctionalsbydirectly
∞
waveenergycutoffs. ThreeclassesofXCfunctionalsarejux- applyingthedifferencebetweenDFTdielectricconstantswith
taposedinthiswork: (i)PBE,(ii)hybrid-DFT,and(iii)pSIC- andwithoutlocalfieldeffects. Computationalinputparame-
DFT.ThetwolatterfunctionalsalsobuilduponthePBEXC ters and calculated properties such as lattice constants, band
functional.20Forthealkalihalidesweemployed216-atomsu- gaps, and dielectric constants are listed in Table II. We fi-
percellsandtheBrillouinzonewassampledusingtheΓ-point nallynotethatEq.(A.2)isequivalenttotherelaxationenergy
only. FortheCLYCsystemweemployed80-atomsupercells ∆E [N ±1;R;α]intheinfinitedilutionlimit.9
hyb e
witha2×2×2k-pointsampling.
2. Polaronbindingenergies
3. G W calculations
0 0
Polaron formation (or binding) energies are computed in
thepSICformalismsimplyas
Quasiparticle G W energies28 were calculated for the al-
0 0
∆E [R ]=E± [R ]−E± [R ], (A.1) kalihalidesintherocksaltcrystalstructureontopofPBEus-
f pol pSIC pol pSIC ideal
ingΓ-centered6×6×6Monkhorst-Packgridsasimplemented
whereE± isdefinedinEq.(3). intheViennaabinitiosimulationpackage.44Weincluded192
pSIC
Conventional hybrid calculations determine polaron bind-
TABLEII. Latticeconstants,bandgaps,anddielectricconstantsfor
ingenergiesthroughenergydifferencesbetweenchargedand
thealkalihalidesandCLYC.Thedielectricconstantsforthelatterare
neutralsupercells(seeRefs.9and41)accordingto
slightly anisotropic and has been averaged over the three principal
axes.
∆E =E [±1,R ]−E [0,R ]±εKS, (A.2)
f PBE0 pol PBE0 ideal ±
Bandgap(eV)
which requires image charge corrections to be applied. This a(A˚) PBE PBE0 G0W0 αEX ε
causes further uncertainties as there is no exact correction LiF 4.065 8.84 11.85 13.24 0.37 10.55
procedure.42 Note that in contrast the pSIC method relies on LiCl 5.090 6.46 8.57 9.08 0.31 9.91
LiBr 5.452 5.07 6.97 7.30 0.29 14.45
the(electron)chemicalpotentialsintheneutralreferencestate
LiI 5.984 4.32 5.95 5.98 0.25 19.43
forobtainingpolaronenergies, andthusisnotsubjecttoim-
NaF 4.537 6.64 9.60 11.60 0.41 4.25
agechargecorrections. Inthecaseofthepresenthybridcal-
NaCl 5.578 5.22 7.33 8.42 0.38 5.66
culations, thespuriousimagechargeinteractionistakeninto
NaBr 5.916 4.24 6.15 6.94 0.35 5.89
accountviathemodifiedMakov-PaynecorrectionofLanyand
NaI 6.407 3.72 5.39 5.80 0.31 7.07
Zunger41,43
CLYC 10.485 4.89 7.12 − 0.25 9.85
2 M
∆E = , (A.3)
corr 32εL
where M is the Madelung constant and L is the lattice
bandstoconvergethedielectricfunctionusedinthescreened
constant of the supercell. The dielectric constants used in
interactionW togetherwithafrequencygridwith48points.
0
Eq.(A.3)isgivenby
Quasiparticle G W calculations were also performed for
0 0
ε=ε +ε , (A.4) the self-trapped hole polaron in NaI. They employed a 64-
∞ ion
atomsupercellandcalculatedG W correctionsontopofthe
0 0
where the first and second terms denote the electronic and PBE0 hybrid functional using 1080 bands for the dielectric
ionic contribution, respectively. We approximately include functiontogetherwith24pointsinthefrequencygrid.
∗ [email protected] 6 S.LanyandA.Zunger,Phys.Rev.B80,085202(2009).
1 T. Maxisch, F. Zhou, and G. Ceder, Phys. Rev. B 73, 104301 7 B.J.MorganandG.W.Watson,Phys.Rev.B80,233102(2009).
(2006). 8 P.R.L.Keating,D.O.Scanlon,B.J.Morgan,N.M.Galea, and
2 R. T. Williams and K. S. Song, J. Phys. Chem. Solids 51, 679 G.W.Watson,J.Phys.Chem.C116,2443(2012).
(1990). 9 P.Erhart, A.Klein, D.A˚berg, andB.Sadigh,Phys.Rev.B90,
3 J.L.Gavartin,P.V.Sushko, andA.L.Shluger,Phys.Rev.B67, 035204(2014).
035108(2003). 10 V.I.Anisimov,J.Zaanen, andO.K.Andersen,Phys.Rev.B44,
4 J.P.Perdew, M.Ernzerhof, andK.Burke,J.Chem.Phys.105, 943(1991).
9982(1996). 11 W.Ka¨nzig,Phys.Rev.99,1890(1955).
5 J.Heyd,G.E.Scuseria, andM.Ernzerhof,J.Chem.Phys.118, 12 P.Zawadzki,K.W.Jacobsen, andJ.Rossmeisl,Chem.Phys.Lett.
8207(2003),erratum:ibid.124,219906(2006). 506,42(2011).
10
13 P.Zawadzki,J.Rossmeisl, andK.WedelJacobsen,Phys.Rev.B 28 L. Hedin, Phys. Rev. 139, A796 (1965); L. Hedin and
84,121203(2011). S. Lundqvist, Solid State Phys. 23, 1 (1970); W. G. Aulbur,
14 F.GygiandA.Baldereschi,Phys.Rev.Lett.62,2160(1989). L.Jo¨nsson, andJ.W.Wilkins,SolidStatePhys.54,1(2000).
15 B. Sadigh, P. Erhart, D. A˚berg, A. Trave, E. Schwegler, and 29 V.FiorentiniandA.Baldereschi,Phys.Rev.B51,17196(1995).
J.Bude,Phys.Rev.Lett.106,027401(2011). 30 R. D. Popp and R. B. Murray, J. Phys. Chem. Solids 33, 601
16 P.Cortona,Phys.Rev.B46,2008(1992). (1972).
17 K.TeegardenandG.Baldini,Phys.Rev.155,896(1967). 31 C.J.Delbecq,J.P.Stokes, andP.H.Yuster,Phys.Stat.Sol.B80,
18 F.C.Brown, C.Ga¨hwiller, H.Fujita, A.B.Kunz, W.Scheifley, 471(1977).
andN.Carrera,Phys.Rev.B2,2126(1970). 32 A. L. Shluger, L. N. Kantorovich, E. N. Heifets, E. K.
19 P.Erhart,A.Schleife,B.Sadigh, andD.A˚berg,Phys.Rev.B89, Shidlovskaya, andR.W.Grimes,JournalofPhysics:Condensed
075132(2014). Matter4,7417(1992).
20 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 33 G.Henkelman,B.P.Uberuaga, andH.Jo´nsson,J.Chem.Phys.
3865(1996),erratum,ibid.78,1396(E)(1997). 113,9901(2000).
21 NotethattheeigenvaluesoftheexactOEParenotequivalentto 34 T.PawlikandJ.-M.Spaeth,J.Phys.Cond.Matter9,8737(1997).
thederivativeoftheenergyfunctionalwithrespecttotheelectron 35 K.BiswasandM.-H.Du,Phys.Rev.B86,014102(2012).
occupations.Thedifferencemakesuptheso-calledderivativedis- 36 A. Bessie`re, P. Dorenbos, C. W. E. van Eijk, L. Pidol, K. W.
continuity,seeRefs.27. Kra¨mer, andH.U.Gu¨del,J.Phys.Cond.Matter16,1887(2004).
22 P.Mori-Sanchez,A.J.Cohen, andW.T.Yang,Phys.Rev.Lett. 37 P.Rinke,A.Janotti,M.Scheffler, andC.G.vandeWalle,Phys.
100,146401(2008). Rev.Lett.102,026402(2009).
23 I.Dabo,A.Ferretti,N.Poilvert,Y.Li,N.Marzari, andM.Co- 38 S.LanyandA.Zunger,Phys.Rev.B81,113201(2010).
coccioni,Phys.Rev.B82,115121(2010). 39 P. E. Blo¨chl, Phys. Rev. B 50, 17953 (1994); G. Kresse and
24 I. Dabo, A. Ferretti, and N. Marzari, in First Principles Ap- D.Joubert,Phys.Rev.B59,1758 (1999).
proachestoSpectroscopicPropertiesofComplexMaterials,Top- 40 G.KresseandJ.Hafner,Phys.Rev.B47,558(1993); G.Kresse
ics in Current Chemistry No. 347, edited by C. D. Valentin, andJ.Furthmu¨ller,Comput.Mater.Sci.6,15(1996).
S.Botti, andM.Cococcioni(Springer,2014)pp.193–233. 41 S.LanyandA.Zunger,Phys.Rev.B78,235104(2008).
25 M.S.HybertsenandS.G.Louie,Phys.Rev.Lett.55,1418(1985). 42 H.KomsaandA.Pasquarello,Phys.Rev.B84,075207(2011).
26 R.W.Godby,M.Schlu¨ter, andL.J.Sham,Phys.Rev.Lett.56, 43 G.MakovandM.C.Payne,Phys.Rev.B51,4014(1995).
2415(1986). 44 M. Shishkin and G. Kresse, Phys. Rev. B 74, 035101 (2006);
27 L.J.ShamandM.Schlu¨ter,Phys.Rev.Lett.51,1888(1983);J.P. Phys.Rev.B75,235102(2007);F.Fuchs,J.Furthmu¨ller,F.Bech-
PerdewandM.Levy,Phys.Rev.Lett.51,1884(1983). stedt,M.Shishkin, andG.Kresse,76,115109(2007).