Fundamental high pressure calibration from all-electron quantum Monte Carlo calculations

1.3 MB·English

by  K. P. Esler

Fundamental highpressure calibrationfrom all-electronquantum MonteCarlo calculations
K. P. Esler,1,2 R. E. Cohen,1 B. Militzer,3 Jeongnim Kim,2 R.J. Needs,4 and M.D. Towler4
1Carnegie Institution of Washington, Geophysical Laboratory
2University of Illinois at Urbana-Champaign, NCSA∗
3Universityof California, Berkeley, Dept. of Earthand PlanetaryScience and of Astronomy
4TheoryofCondensed MatterGroup, CavendishLaboratory, Cambridge CB30HE,UnitedKingdom
(Dated:January13,2010)
Wedevelopanall-electronquantumMonteCarlo(QMC)methodforsolidsthatdoesnotrelyonpseudopo-
0 tentials, and use it toconstruct a primaryultra-high pressure calibration based theequation of stateof cubic
1 boronnitride(c-BN).WecomputethestaticcontributiontothefreeenergywithQMC,andobtainthephonon
0 contributionfromdensityfunctionaltheory,yieldingahigh-accuracycalibrationupto900GPausabledirectly
2
inexperiment. Furthermore, wecompute theanharmonic Ramanfrequency shiftwithQMCasafunctionof

n pressureandtemperature,allowingopticalpressurecalibrationintable-topexperiments. Incontrasttopresent
a experimentalapproaches,smallsystematicerrorsinthetheoreticalEOSdonotincreasewithpressure,andno
J extrapolationisneeded. Thisall-electronmethodology isgenerallyapplicabletofirst-rowsolids, andcanbe

3 usedtoprovideanewreferenceforabinitiocalculationsofsolidsandtobenchmarkpseudopotentialaccuracy.
1
PACSnumbers:64.30.Jk,81.05.Je,02.70.Ss,63.20.dk

]
i
c
s Althoughthe numberofstudiesofmaterialsusingdensity ulus, in conjunctionwith X-ray diffraction measurementsof
- functional theory (DFT) continues to grow explosively, the volume, can be integrated to provide an EOS, but a correc-
l
r accuracyoftheirpredictionsisvariable,sometimesexcellent tionmustbemadetotransferfromanadiabatictoanisother-
t
m and sometimes poor, limiting the confidence one may place malpath[3]. Newapproachessuchasquasi-adiabaticZ-pinch
. in DFT results as a quantitative calibration for experimen- based experiments[4, 5] also hold future promise for a pri-
t
a tal studies. Quantum Monte Carlo (QMC), particularly dif- mary scale. There have been attempts to refine the ruby
m fusion Monte Carlo (DMC), is the highest-accuracy method scale[6], and new calibrations have also been suggested[3].
- forfindingthe groundstate ofa many-electronHamiltonian. Cubicboronnitride(c-BN)hasbeenidentifiedasapromising
d
ForsolidswithatomsheavierthanHe,however,theHamilto- materialfora newscale[7]. InthisLetter, we providea new
n
nianitselfisapproximatedusingpseudopotentials(PPs)based pressurescalebasedonDMCcalculationsoftheEOSandRa-
o
c onalower-accuracytheory,limitingitsreliability,andaswe manfrequencyofc-BN.Thistheoreticalapproachhasthead-
[ show,commonlyusedPPsgivedisparateresults. InthisLet- vantagethat, in contrastto presentexperimentalapproaches,

ter, we pushthe state ofthe artin accuracyby introducinga the methodworksequallywellunderhighcompression,and
1
v method for an all-electron DMC simulation of solids, elimi- uncertaintydoesnotgrowwithpressure.
9 natingthebiasfrompseudopotentials,andapplyittocreatea
PressureisthenegativevolumederivativeoftheHelmholtz
7 high-accuracypressurecalibrationscale.
0 free energy. In a wide-gap insulator such as c-BN, the free
2 The combination of ultra-high pressure mineralogy with energycanbewrittenasasumofthefrozenlatticeenthalpy,
. seismology has yielded a wealth of insight into the internal dependentonlyonvolume,andaphononthermalfreeenergy,
1
0 structure of our planet. Pressure is the key that links these which depends on both volume and temperature. Since the
0 disciplines,mappingphasetransitionsandmineralproperties static enthalpy is the dominantcontributionat ordinarytem-
1 to planetary depth. Establishing an absolute pressure cali- peratures, errors in a theoretical EOS can most often be at-
:
v bration at multi-megabarpressures, however, poses a funda- tributed to the static part. Previous calculations of the EOS
i mentalandcontinuingproblemforhigh-pressureexperiment. ofc-BNhavebeenbasedonDFT[7],whichuseapproximate
X
Current primary calibrations are based on data from shock- functionalstotreatelectronexchangeandcorrelation.Several
r
a waveexperiments,whichinferpressurefromconservationof functionalsareincommonuse,eachgivingrisetoadifferent
momentumandenergyastheshocktraversesthesample.Ex- EOS,andthereis noa prioriway to predictwhichwillgive
perimentaluncertainty,combinedwitherrorsinmodelextrap- themostreliableresult.
olation,yieldscaleswhichdifferfromeachotherbyasmuch
QuantumMonte Carlo simulationmethodsexplicitly treat
to7%atroomtemperature,withevengreaterdiscrepanciesat
electron exchangeand correlationinstead of resorting to ap-
hightemperature[1]. Suchdisparityremainsaseriousobsta-
proximatefunctionals. VariationalMonteCarlo(VMC)com-
cletoaquantitativeunderstandingofEarth’sinterior.
putespropertiesbyaMetropolissamplingofatrialwavefunc-
Apressurecalibrantisamaterialwithaknownequationof tion. DiffusionMonteCarlo samplesthe many-bodyground
state(EOS),asmallsampleofwhichmaybeincludedinhy- state of the Hamiltonian through a stochastic projection of
drostaticequilibriumwith the test subject(e.g. ruby[2]). As the trial function. In practice, a fixed node approximationis
an alternative to shock experiments, high-pressure Brillouin usedforfermions,suchaselectrons,whichtestshaveshown
scattering, which provides a measurement of the bulk mod- giverelativelysmallerrorwhenthe nodesare obtainedfrom
2
high-qualityDFT orbitalsfor electronicallysimple materials (a)Parametersfor300KEOS
suchasc-BN.DMCforsolidshasbeendemonstratedtogive Source PP/AEatoms V0(A˚3) B0(GPa) B0′ (nondim.)
significantly more accurate cohesive energies[8], equations
Trail-NeedsPP+AE 64/8 11.792(18) 381(6) 3.87(6)
ofstateandRamanfrequencies[9], andphasetransitions[10]
OPIUMGGAPP+AE 64/8 11.769(17) 385(6) 3.86(6)
than DFT. We have used both the CASINO QMC software BFDHFPP+AE 64/8 11.781(20) 382(7) 3.87(7)
suite[11]andQMCPACK[12]. BFDHFPP+AE 128/16 11.812(8) 378(3) 3.87(3)
QMCsimulationsofsolidsarecurrentlyperformedwithin Datchietal.[24] 11.8124 395(2) 3.62(5)
the pseudopotential (PP) approximation, in which the core Goncharovetal. 11.817(32) 387(4) 3.06(15)
electronsareeliminatedandtheireffectreplacedbyanonlo-
calpotentialoperator[8]. SincePPsarepresentlyconstructed (b)Parametersforthermalpressure
with a lower-accuracy theory, such as Hartree-Fock (HF) or
DFT,thisreplacementrepresentsanuncontrolledapproxima- θn0(KA˚−3n) αn(K−1A˚−3n) βn(KA˚−3n)
n=0 4.836656×103 2.598608×10−3 -4.869563×103
tion. To eliminate this error, we develop a method for all-
n=1 -6.929704×101 -1.200504×10−4 1.763443×102
electron(AE)QMCsimulationsofsolidsinQMCPACKusing
n=2 4.634278×10−1 2.789128×10−6 -2.186053×100
trialwavefunctionsderivedfromfull-potentiallinearizedaug-
n=3 -1.273468×10−3 -2.008994×10−8 9.303181×10−3
mented plane wave (FP-LAPW) calculations using the EX-
CITING code[13]. Space is divided into spherical muffin
(c)ParametersforRamancalibration
tin regions around the nuclei, and an interstitial region. Or-
bitalsarerepresentedinsidethemuffintinsasaproductofra- c0(cm−1) c1(cm−1) c2(K) b R0(GPa) R1(GPa) R2(GPa)
dialfunctionsandsphericalharmonics,andoutsideasplane- 1055.9 -144.24 1497.8 3.0155 349.87 1849.4 112.33
waves. Toensurethatawavefunctionsatisfiesthevariational
principle,itmustbebothcontinuousandsmooth.Weutilizea
TABLEI:Parametersforthec-BNEOSandRamancalibration.
super-LAPWformalismthatenforcescontinuityandsmooth-
ness atthe muffin tin boundary. For efficiency, we represent
theorbitalsas3DB-splinesintheintersticesandtheproduct Wecomputethefreeenergyforourc-BNsystemattwelve
ofradialsplinesandsphericalharmonicsinthemuffintins. unit-cellvolumes, spanningvolumecompressionratiosfrom
Since AE QMC simulations are computationally expen- 0.84 to 2.0, corresponding to pressures of about -50 GPa to
sive,weperformthesesimulationsin8-atomcubicsupercells. 900GPa. Whileitisnotfullycertainthatthecubicphaseof
Simulation cells this small would typically have significant BNisstableto900GPa,nootherstructurehasbeenobserved,
finite-sizeerrors.Weeliminatetheseerrorsbycombiningdata andtheoreticalstudieshavenotidentifiedanytransitionsbe-
fromAEsimulationswiththatfromPPsimulationsperformed low 1 TPa[22]. We use the Vinet form[23] for the isother-
inboth8-atomand64-atomsupercells. Thus,weareableto malEOS,whichwefindrepresentsourfree-energydatavery
simultaneouslyeliminatesystematicerrorsfrompseudpoten- well, yielding the bulk modulus, B0, its pressure derivative,
tials and from finite-size effects. The corrected static-lattice B0′,andtheequilibirumvolume,V0 (TableI(a)). Thestatisti-
energyisgivenateachvolumeas calerrorbarsforeachdatapointaredirectlydeterminedfrom
ourQMCdata.Wecomputestatisticalconfidenceranges,tak-
E =EPP+ EAE EPP +∆MPC+∆kinetic, (1) ing into account parameter cross-correlations with a simple
64 8 8 64 64
− MonteCarloprocedure.
where the term in(cid:2)brackets re(cid:3)moves the psuedopotential Figure 1 shows the EOS of c-BN at 300 K, with experi-
bias. ∆MPC and ∆kinetic are, respectively, potential[14] and mentaldatafromDatchietal.[24]andGoncharovetal.[7],as
64 64
kinetic[15] corrections for finite-size errors. We perform wellasthepresentworkfromsimulationswiththreedifferent
this procedure with three different PP sets commonly used PPs. The residuals in (c) are derived from DMC simulation
in QMC: HF PPs from Trail and Needs (TN)[16, 17]; HF with PPs alone, while those in (b) combine all-electron and
PPs fromBurkatzki, Filippi, andDolg(BFD)[18]; andDFT- PP data. There is significant discrepancy between the theo-
GGA[19] PPs generated with OPIUM (WC)[20]. Perform- reticalcurvesin(c),suggestingthatPPsimulationalonedoes
ing the same procedure with 128-atom PP simulations and notprovidesufficientaccuracy.OncethePPdataiscombined
16-atomAE simulationsyieldsstatistically indistinguishable with AE data, as in (b), all the theoretical curves come into
results, demonstratingthatfinite-size errorsare convergedat goodagreement. OurtheoreticalEOSagreesreasonablywith
thissize. Additionalmethodologicaldetailsandthefinite-size thatinRefs. [24]and[7]withintheexperimentallymeasured
datacanbefoundinEPAPSdocumentNo. #. pressurerange,buttheexperimentalextrapolationshowssig-
To compute the phonon free energy, we use density nificantdeviationathighpressure.
functional perturbation theory (DFPT) in the QUANTUM Wemaywritethethermalequationofstateintheform
ESPRESSO package[21] with the Wu-Cohen functionaland
P(V,T)=P300K(V)+Pth(V,T) Pth(V,T =300), (2)
the OPIUM PPs. The phonon density of states, from which −
wederivethermodynamicdata,isusuallyverywell-described whereP300K is the room-temperaturecontribution,fitto the
withDFT. Vinet form. The phonon contribution is written in an aug-
3
(a)CorrectedEOS (b)Correctedpressureresiduals (c)Uncorrectedpressureresiduals
FIG.1:Thec-BNEOSat300K,computedwithDMC,comparedtoexperiment.WeplotthefullEOSandthepressureresidualswithrespect
toDMCwiththeBFDPP.(a)givesthecorrectedEOSresultingfromcombiningtheEOSandPPdata,while(b)givesthepressureresiduals
for the same data. (c) gives the uncorrected pressure residuals from the PP simulations only. Shaded areas represent the one-σ statistical
confidenceregionfromQMCdata.
mentedDebyemodel, temperature. The matrix element for the transition from n
ton 1isproportionalto√n andisthermallyweightedby
∂F (θ,T) ∂θ −
P (V,T) = D (3) theBoltzmannoccupationofstaten,sotheintensity-averaged
th
− ∂θ ∂V frequency, ν ,canbegivenas
0 h i
θ(V,T) = θ (V)+β(V)exp( α(V)T) (4)
x(V) = x0+x1V+x2V2+x−3V3, x∈{θ0,α,β}(5) ν = ∞n=1I∞nEn−hEcn−1, where In =ne−kEBnT. (7)
h i I
P n=1 n
in which the Debye temperature, θ, is a function of both V
andT (TableI(b)). TheDebyefreeenergypertwo-atomcell, The excess thePrmal softening, i.e. beyondthat fromthermal
excludingthezero-pointterm,isgivenby expansion alone, is accounted for by the thermal average of
the anharmonic frequencies. Figure 2 shows the computed
θ θ 3 Tθ x3 Raman frequenciescomparedwith the experimentaldata re-
FD =6kBT "ln(cid:16)1−eT(cid:17)−(cid:18)T(cid:19) Z0 ex−1dx#. (6) poBrteodthinRRefesf.s.[[2255]–a2n7d].[26] give a ruby-like calibration for-
c-BN can be used to calibrate pressure optically by mea- mula,whichcanbeexpressedas
suring the frequencyshift of the TO Raman mode, allowing
bench-topexperiments.Wecomputethepressureandtemper- P =(R/b) (ν/ν¯)b 1 , (8)
−
aturedependenceofthisfrequency.Withinthequasiharmonic
(cid:2) (cid:3)
approximation,phononfrequencieshaveexplicitdependence where R, ν¯, and b have quadratic T-dependence. This de-
pendenceissufficientbelow2000K,butcannotrepresentour
onvolumeonly.Atconstantpressure,however,thesefrequen-
data at high temperature. We use a form which capturesthe
cies have an implicit temperaturedependenceresulting from
Boltzmannoccupationofphononexcitations,
thermal expansion. The dependence due to thermal expan-
sion accounts for only about half the total T-dependenceof
ν2(P)
theRaman mode,asobservedin [25] and[26]. Theremain- ν(P,T) = ν0(P)+ν1(P)exp (9)
− T
ingT-dependencecanbeattributedtosignificantanharmonic (cid:20) (cid:21)
effectsinc-BN,andisincludedinourcalculations. bP 1/b
ν (P) = c +1 , n=1,2,3 (10)
Sincetheopticalbranchhassmalldispersion,we treatthe n n R
anharmonicityasaone-dimensionalon-siteanharmonicoscil- (cid:20) n (cid:21)
lator,inasimilarapproachtotheQMCcomputationoftheTO withparametersinTableI(c)andplottedinFigure2(a). Note
Raman frequencyof diamond[9]. At each volume, we com- that this formula cannot be analytically inverted, but a very
pute the effective Born-Oppenheimer potential well for the simple iterative solutioncan be used for calibratingpressure
TO modewith DMC and the BFD PP at nine displacements fromν andT.
alongthe modeeigenvectorin the 64-atomsupercell. We fit The main axis of Figure 2(b) gives the room-temperature
thedatatoaquarticpolynomial,andnumericallysolvethe1D Raman frequencyversus pressure. There is good agreement
Schro¨dingerequationinthisanalyticpotential.Thisresultsin intherelativelylow-pressureregioninwhichtheRamanfre-
asetofsingle-phononenergylevels, E ,withnonuniform quency was measured. At very high pressure, the deviation
n
{ }
separation. From E , we compute an intensity-weighted withrespecttotheextrapolationin[25]increaseswithamax-
n
average Raman fre{quen}cy, ν¯, as a function of pressure and imum discrepancy of 38 cm−1 or, conversely, a deviation in
4
andcan be used to cross-calibratescales based on otherma-
terials. Sincetheaccuracyofourmethodsshouldnotdepend
oncompression,ourcalibrationcanbeusedupto900GPa.
ThisworkwassupportedundertheNationalScienceFoun-
dationGrantEAR-0530282. Thisresearchusedresourcesof
theNationalCenterforComputationalSciencesandtheCen-
ter for Nanophase Materials Sciences, which are sponsored
by the U.S. Department of Energy. This work was partially
supportedbytheNationalCenterforSupercomputingAppli-
cationsunderMCA07S016andutilizedtheAbeandLincoln
(a)FittopresenttheoryforRamanshift. clusters. We thank Alexander Goncharov, F. Mauri and M.
Lazzeriforhelpfulsuggestionsanddiscussions.
∗ Electronicaddress:[email protected]
[1] YingweiFeietal.,PNAS104,9182(2007).
[2] H.K. Mao, J. Xu, and P.M. Bell, J. Geophys. Res. 91, 4573
(1986).
[3] Chang-Sheng Zha, Ho-kwang Mao, and Russell J. Hemley,
PNAS97,13494(2000).
[4] J.R.Asayaetal.,Int.J.ofImpactEngineering23,27(1999).
[5] J.RemoandM.Furnish,Int.J.ofImpactEngineering35,1516
(2008).
[6] W.B.Holzapfel,HighPressureResearch25,87(2005).
[7] AlexanderF.Goncharovetal,Phys.Rev.B75,224114(2007).
[8] W.M.C.Foulkes,L.Mitas,R.J.Needs,andG.Rajagopal,Rev.
(b)ComparisonofRamanexperimentandtheory.
Mod.Phys.73,33(2001).
[9] RyoMaezono,A.Ma,M.D.Towler,andR.J.Needs,Phys.Rev.
FIG.2: TORamanfrequencyofc-BNasafunctionoftemperature
Lett.98,025701(2007).
and pressure. (a) gives compares our theoretical Raman frequency
[10] J.KolorencˇandL.Mitas,Phys.Rev.Lett.101,185502(2008).
withthefittedforminEqs.9and10. (b)comparesour calibration
[11] R.J.Needs,M.D.Towler,N.D.Drummond,andP.L.Rios,J.
at300Kwiththeexperimentalresultsfrom[25]and[26].Theblack
Phys.:CondensedMatter22,023201(2010).
dotswitherrorbarsaretheQMC frequencies, whilethebluesolid
[12] JeongnimKim,K.Esler,J.McMinis, B.Clark, J.Gergely,S.
line gives the fit. The lower inset gives an expanded view at low
Chiesa,K.Delaney,J.Vincent,andD.Ceperley.
pressure,inwhichtheshadedregiongivesthestatisticalconfidence
[13] J. K. Dewhurst et al., The EXCITING FP-LAPW code,
region.Theupperinsetgivesthevariationwithtemperature.
http://exciting.sourceforge.net.
[14] A.J.Williamsonetal.,Phys.Rev.B55,R4851(1997).
[15] S. Chiesa, D.M. Ceperley, R.M. Martin, and M. Holzmann,
the pressure calibration of 50 GPa at 900 GPa. The devia-
Phys.Rev.Letts.97,076404(2006).
tions with respect to [26] are 70 cm−1 and 120 GPa. The [16] J.R.TrailandR.J.Needs,J.Chem.Phys.122,174109(2005).
experimentalparameterscapturethecorrectqualitativehigh- [17] J.R.TrailandR.J.Needs,J.Chem.Phys.122,014112(2005).
pressurebehaviorupto900GPa,despitethefactthatdatawas [18] M. Burkatzki, C. Filippi, and M. Dolg, J. Chem. Phys. 126,
availableonlyto 20 and64GPa, respectively. Thissuggests 234105(2007).
[19] ZhigangWuandR.E.Cohen,Phys.Rev.B73,235116(2006).
theformforthefitwaswell-chosen.
[20] Opiumpseudopotentialgenerationproject,
We have presented a fully ab initio pressure calibration
http://opium.sourceforge.net.
based on quantumMonte Carlo simulations, and have intro-
[21] S.Baronietal.,J.Phys.:Condens.Matter21,395502(2009).
ducedamethodforall-electronsimulationsofsolidstoelimi- [22] J. Furthmu¨ller, J. Hafner, and G. Kresse, Phys. Rev. B 50,
natebiasfrompseudopotentials.Thismethodshouldbeappli- 15606(1994).
cabletoatleastfirst-rowsolids, allowingincreasedaccuracy [23] PVinetetal.,J.Phys.:CondensedMatter1,1941(1989).
inthestudyofothermaterialsandprovidinganewbenchmark [24] F. Datchi, A. Dewaele, Y. Le Godec, and P. Loubeyre, Phys.
Rev.B75,214104(2007).
forothermethods.Theonlyremainingsystematicerrorinthe
[25] F.DatchiandB.Canny,Phys.Rev.B69,144106(2004).
staticcontributiontotheEOSisfromthefixed-nodeapproxi-
[26] Alexander F. Goncharov, Jonathan C. Crowhurst, J. K. De-
mationusedinDMC.Forsimplematerialssuchasc-BNthis
whurst,andS.Sharma,Phys.Rev.B72,100104(R)(2005).
errorshouldbequitesmall, andtendstocancelbetweendif- [27] I. V. Alexandrov et al., High pressure study of diamond,
ferentvolumes. ThuswebelievetheEOSisrobustenoughto graphite and related materials (Terrepub AGU, Washington-
beuseddirectlyinexperimentasaprimarypressurecalibrant, Tokyo,1992),pp.409–416.