Table Of Content

Momentum-space finite-size corrections for Quantum-Monte-Carlo calculations R. Gaudoin1, I. G. Gurtubay2,3, and J. M. Pitarke2,4 1Faculty of Physics, University of Vienna, and Center for Computational Materials Science, Sensengasse 8/12, A-1090 Vienna, Austria 2Materia Kondentsatuaren Fisika Saila, Zientzia eta Teknologia Fakultatea, Euskal Herriko Unibertsitatea, UPV/EHU, 644 Posta kutxatila, E-48080 Bilbo, Basque Country, Spain 3Donostia International Physics Center (DIPC), E-20018 Donostia, Basque Country, Spain 4CIC nanoGUNE Consolider and Centro F´ısica Materiales (CSIC-UPV/EHU), Tolosa Hiribidea 76, E-20018 Donostia, Basque Country, Spain (Dated: January 31, 2012) 2 1 Extended solids are frequently simulated as finite systems with periodic boundary conditions, 0 which due to the long-range nature of the Coulomb interaction may lead to slowly decaying finite- 2 size errors. In the case of Quantum-Monte-Carlo simulations, which are based on real space, both n real-spaceandmomentum-spacesolutionstothisproblemexist. Here,wedescribeahybridmethod a which using real-space data models the spherically averaged structure factor in momentum space. J Weshowthat(i)byintegrationourhybridmethodexactlymapsontothereal-spacemodelperiodic 0 Coulomb-interaction(MPC)methodand(ii)thereforeourmethodcombinesthebestofbothworlds 3 (real-space and momentum-space). One can use known momentum-resolved behavior to improve convergence where MPC fails (e.g., at surface-like systems). In contrast to pure momentum-space ] methods,ourmethodonly dealswith a simplesingle-valued function and,hence,betterlendsitself r tointerpolationwithexactsmall-momentumdataasnodirectionalinformationisneeded. Byvirtue e h of integration, theresulting finite-sizecorrections can be written as an addition to MPC. t o PACSnumbers: 71.10.Ca,71.15.-m,73.20.-r . t a m I. INTRODUCTION length scale. Although QMC can, in principle, model - correlationswellatshortlengthscales,itisboundtofail d atlargerscalesduetothefactthatfinitesimulationcells Theissueoffinite-sizecorrectionsinQuantum-Monte- n must be used. o Carlo (QMC) calculations has recently been attracting c considerableattention.1–5 As QMC calculations of solids So far we have only touched on the issue of finite-size [ need to be carried out within a supercell, the Coulomb corrections for the interaction energy. The information 1 interactionistypicallyreplacedbytheso-calledEwaldin- to evaluate this is contained in the spherical average of v teractionthatiscompatiblewiththesupercellgeometry.6 the diagonal terms of the two-particle density matrix 4 This, in turn, is equivalent to dealing with an infinite (the two particle density) and hence Sk, which contains 1 systemwithaperiodicallyrepeatedexchange-correlation the same information, suffices. However, other quanti- 2 (xc) hole. In other words, using the Ewald interaction tiesmaybe computedusing QMCwhichalsosuffer from 6 includes the spurious effective interaction of an electron finite-size errors. E.g. Chiesa al.2 also deal with correc- . 1 with its periodically repeated xc hole. One solution to tions to the finite-size errors in the kinetic energy but 0 this drawback - the real-space solution - is to use QMC since the kinetic energy needs information that goes be- 2 dataatlengthscalessmallerthanthesupercellsizewhere yond the diagonal terms of the density matrix such an 1 QMCisexpectedtobeaccurateandsubstitutethemiss- analysisis beyondthe scope ofthis paper which is based : v ing terms implicitly or explicitly. E.g., the Model Peri- solelyonthe informationinSk. Similarly,finite-size cor- Xi odicCoulomb-interaction(MPC)method1dealswiththe rections to the momentum distribution3 cannot readily periodically repeatedxc hole by using the bare Coulomb be done using just Sk. r a term within the supercell and assuming that beyond the TheaimofthispaperistoshowthatmodellingS with k supercellanelectrononlyfeelstheHartreepotentialwith theMPCansatzyieldsamomentum-resolvedMPC.This no further correlations. This represents a good approxi- allowsforanintuitiveanalysisoffinite-sizeerrorsandfor mation in a bulk solid, where the xc hole decays rapidly. animprovementofMPC wheneverMPCfails, suchasin In contrast, Chiesa al.2 traced the Ewald error back to surface-likesystems. Firstof all, in SectionII we discuss an integration error and then added back the missing S in the context of QMC and apply the MPC ansatz, k contribution, yielding a momentum-space solution. thereby showing that via integration over k our method In a recent paper,4 the spherically-averaged structure exactly maps onto MPC. In Section III, we first look at factor S was introduced in order to analyze and reduce a non-interactinguniform electrongas in order to (i) see k Coulomb finite-size effects. This (S ) represents a natu- how MPC breaks down in this case and (ii) learn how k ralquantity to study finite size errors: Onthe one hand, to overcome the limitation of MPC by using our hybrid it is a simple one-dimensional function which via inte- analysis. We also report Diffusion Monte Carlo (DMC) gration yields the interaction energy and, on the other calculations of the structure factor of an interacting uni- hand, it naturally orders the QMC data according to form electron gas, which yields an intuitive understand- 2 ing of the MPC in interacting systems. This then leads L us to an easy way to understand the existing difficulties MPC that arise in the case of a slab geometry,7,8 and we pro- n nMrePdC pose a simple method to address this difficult problem. nTRUE TRUE Weroundoffthepaperwithasummaryandconclusions. red We use atomic units throughout (¯h=e2 =m =1). e nEWALD red II. MODELLING THE EXCHANGE CORRELATION HOLE IN MOMENTUM SPACE A. The interaction energy EWALD The interaction energy Uint of an arbitrary system of u=0 u=L many interacting electrons can be expressed as follows ′ FIG. 1. The reduced electron density n (r,r) seen by an red Uint =UH +Uxc, (1) electronatr,asafunctionofu=|r−r′|. Intheabsenceofex- changeand correlation, thiswould equaltheelectron density where UH represents the Hartree energy (which in the n(r′) (thin solid line) which in thecase of a uniform electron case of an infinite jellium model is exactly cancelled by gas would be constant. In the presence of exchange and cor- the presence of the positive background) relation, thetruereduced electron densityn (r,r′)is ofthe red form represented by the thick line marked “TRUE”, which e2 n(r′) for asupercell geometry (of lineardimension L) becomesthe U = drn(r) dr′ , (2) H 2 r r′ reduced Ewald density represented by the thick dotted line Z Z | − | marked “EWALD”. The MPC interaction energy Uint is ob- and U represents the so-called xc interaction energy tained by using a cutoff reduced electron density (thick solid xc correspondingto the attractive interactionbetweeneach linemarked“MPC”)thatcoincideswith(i)thereducedEwald electron and its own xc hole: electrondensitywithinthesimulationcelland(ii)theelectron ′ density n(r) beyond thesimulation cell. e2 n (r,r′) U = drn(r) dr′ xc . (3) xc 2 r r′ Z Z | − | thesumoftwoterms: TheHartreeenergyU calculated H Here, n(r) is the electron density and n (r,r′) repre- with the Ewald interaction, and the beyond-Hartree xc- xc sents the xc-hole density ofanelectronatr. For brevity, energy Uxc calculated with a cutoff Coulomb interaction we shallalsodefine a reduced electrondensity n (r,r′), usingtheminimumimageconvention,i.e.,translatingco- red whichrepresentstheelectrondensityatr′seenbyagiven ordinatesofelectronpairssuchthatr r′ lieswithinthe electron at r in the presence of (hence, reduced by) its simulationcell. Hence, the MPCinter−actionenergyUint xc hole: is obtainedby simply replacing the true reduced electron density n (r,r′) of Eq. (4) by the MPC reduced elec- nred(r,r′)=n(r′)+nxc(r,r′)= n2n((rr,)r′), (4) tbryonthdeetnhsiirctekydsnoMrleiddPCli(nre,rm′)arokfetdhe“MfoPrmC”d.isplayed in Fig. 1 wheren2(r,r′)istheso-calledtwo-particledensity.9Note that n (r,r′)<0. xc C. Spherically averaged structure factor S k Startingwiththexc-holedensityn (r,r′)atr′around B. Model periodic Coulomb interaction (MPC) xc an electron at r, one finds the following momentum- resolved form of the xc interaction energy of Eq. (3):10 The only periodic solution of Poisson’s equation for a periodic array of charges is the so-called Ewald interac- N U = (S 1)dk, (5) tion, which is of the Coulomb form 1/r only in the limit xc π k− Z of an infinitely large simulation cell. However, while the where S is the spherical average of the diagonal struc- Hartree energy is given correctly by this Ewald inter- k action, the part of the electron-electron energy coming ture factor Sk,k′:4 from the interaction of electrons with their own xc hole 1 sin(ku) yields a spurious contribution that is due to the inter- Sk =1+ drn(r) dr′ nxc(r,r′), (6) N ku action of an electron with its periodically repeated xc Z Z hole. This drawback was solved in Ref. 1 by replacing and u= r r′ . | − | the Ewald interaction by a model periodic Coulomb in- Equation(5)representsageneralexpressionforthexc teraction that yields an interaction energy consisting of interaction energy, which is formally exact not only for 3 homogeneous media but also for an arbitrary inhomoge- There aresomecaveats,however. The non-interacting neous many-electron system. uniform electron gas as well as semi-infinite systems contain a linear contribution to S at k 0 (due to k → the presence of a xc hole that is not negligible even D. Monte-Carlo (MC) sampling of the structure at large distances), which causes the failure of the factor MPC scheme. Nevertheless, in the framework of our momentum-resolvedapproachthereisroomforimprove- When performing MC sampling on a function f(r,r′) ment over MPC, since one has flexibility to go beyond such as the Coulomb potential 1/r r′ , what we are the MPC ansatz by replacing the low-k structure factor actually calculating is | − | SMPC by its known correct value. k f(r,r ) = drn(r) dr′n2(r,r′)f(r,r′), i j III. THE MPC STRUCTURE FACTOR IN * + Xi6=j MC Z Z PRACTICE (7) where the minimum image convention is used. In con- A. The non-interacting uniform electron gas junction with Eq. (4), this yields the MC sampling of the spherically averagedstructure factor of Eq. (6): The non-interacting (Hartree-Fock) uniform electron 1 sin(k r r ) gas was dealt with extensively in Ref. 4. Here we briefly S =1+ | i− j| SH, (8) discuss this system, with the aim of now introducing a k N *Xi6=j k|ri−rj| +MC − k correctionterm. The exact Hartree-Fock (HF) structure factor S of a uniform electron gas, which is easily de- k here SH being the Hartree contribution: rived analytically,11 is shown in Fig. 2 (dashed-dotted k red line) together with the result we obtain by using the 1 sin(ku) SH = drdun(r)n(r u) , (9) MPC ansatz (solid black line). This figure clearly shows k N − ku ZSC that the MPC result is in qualitative error at low k, due to the fact that the actual HF structure factor does not whereagaintheminimumimageconventionisbeingem- exhibit a quadratic behavior in the limit as k 0 but ployed,i.e.,the integralsarecarriedoutoverthe simula- → tion cell (SC). The densities n(r) here are standard MC alinearbehaviorinstead. Ourmomentum-resolvedtech- nique, however,has roomfor improvement,by going be- electron densities. yond the MPC ansatz. Beyond a system-size dependent momentum cutoff k c (i.e.,atk >k ),theQMCelectroncorrelationandhence E. MPC from Sk c structure factor are expected to be accurate. The momentum cutoff k should be of the order of the inverse If one applies the MPC ansatz in the calculation of c of the characteristic length L of the supercell. A sim- the spherically averaged structure factor S , by simply k ple way to model a correction proceeds as follows: After assuming (as in Fig. 1) that beyond the supercell corre- choosinganappropriatecutoffk (verticallineofFig.2), lations are only due to variations in the density, then c we model the corrected HF structure factor as a straight Eq. (5) yields exactly the MPC interaction energy of line between S = 0 and SMPC and the MPC struc- Ref. 1. This can be seen by introducing Eq. (6) into k=0 kc turefactorasaquadraticcurvebetweenthesamepoints. Eq. (5) and performing the k integration. Hence, the ThecorrectiontoMPCisthengivenbytheareabetween momentum-space based method of Ref. 4 reduces, un- these two curves. The dotted blue line and the dashed der the MPC ansatz, to what we may call a momentum- green line of Fig. 2 illustrate this modelling. More real- resolved MPC. istic models are of course feasible. E.g., one could take The MPC ansatz implies a finite extentofthe xc hole, into account derivatives of S at k =0 or k =k , etc. whichinturnresultsinaquadraticbehaviorofS ask k c k → Hence, as long as we are interested in the xc interac- 0. In the case of finite systems (e.g., atoms, molecules, tionenergyU ,thereisnoneedtoactuallyevaluatethe and clusters) and bulk solids, this quadratic behavior of xc spherically averagedstructure factor for all k. In fact, if the structure factor is qualitatively correct because of we assume that k is already in the asymptotic regime theshortrangeofthexcholeinthosesystems. Hence,in c one can base the entire correctionon the value of k and thosesystems,aslongasthesimulationcellissufficiently c the asymptotic form of the structure factor. Let us as- large for the spherically averagedstructure factor at the sume the MPC structure factor below k is of the form cutoff momentum k 1/L to already be in or close to c c the asymptotic low-k∼behavior, the MPC ansatz yields SMPC =βk2+γk3 , (10) k accurateresults: SinceinthosesystemsthetrueS isan k while the true functional dependence (in the asymptotic essentially quadratic function as k 0, constrained by S =0,andthe same constraints→holdforSMPC, only region) ought to be k=0 k higher-order terms contribute to any residual error. S =αk , (11) k 4 1 N=18 N=54 N=102 1 N=178 MPC asymptotic exact 0.5 ) k 0 0.2 0.4 0.6 0.8 1 1.2 S(0.5 S(k)0.5 (k) 0.4 S 0.2 0 0 1 2 0 0 0 2 4 6 8 10 0 0 4 k (au) 2 k FIG. 3. (Color online) DMC spherically averaged structure factor S of an interacting uniform electron gas of r =1 (in k s FIG.2. (Coloronline)Thesphericallyaveragedstructurefac- an fcc simulation cell with N=18, 54, 102, and 178, N being tor S of a non-interacting (Hartree-Fock) uniform electron k thenumberofelectrons),attheMPClevel. Thethickdotted gas with an electron-density parameter r = 1 and 54 elec- s blue line marked “asymptotic” represents the true structure tronsinaface-centered-cubic(fcc)simulationcell. Theexact factor of a uniform electron gas at k → 0: S = k2/2ω .13 k p HFstructurefactorS isrepresentedbyadashed-dottedred k Theverticallinesrepresentthecutoffk =1/LforN=18,54, c line and the MPC result by a solid black line. These two 102, and 178. curves are also plotted in the inset, but now together with Eq. (11) (nearly atop the exact HF structure factor), represented by a blue dotted line, and the modelling of Eq. (10) (nearlyatoptheMPCHFstructurefactor),representedbya the value of the MPC HF structure factor at the cutoff green dashed line. The difference between the corresponding where it ought to coincide with the exact HF structure integralsrepresentsthecorrectiontotheMPCHFdata. The factor. vertical line indicates thecutoff kc. Aspointedoutabove,sofarthisrepresentsacrudeway of devising the finite-size correction to the MPC result. More complex functional forms can be chosen, by using α coming from the linear term in the HF S .11 At the k thestructurefactororevenitsderivativesatk =0ork = cutoff, both should of course coincide: 6 0 to estimate the corresponding parameters. Analytic αk =βk2+γk3 , (12) integration would then yield a more accurate correction c c c term. as should their derivatives: α=2βk +3γk2, (13) c c B. The interacting uniform electron gas which yields β =2α/k and γ = α/k2. In this model, the ccorrected st−ructucre factor Sk only InRef.4,VariationalMonteCarlo(VMC)calculations differs from the MPC structure factor at momenta be- of the structure factor S and the xc interaction energy k tween zero and kc, so the correctionin Uxc can easily be Uxc ofaninteractinguniformelectrongaswerereported. evaluated using Eq. (5); it turns out to be: Here we report DMC calculations of these quantities, as obtained by following the Hellmann-Feynman sampling 1 ∆U =Ucorrected UMPC = αk2. (14) introduced in Ref. 12. xc xc − xc 12 c Figure 3 exhibits the DMC structure factor S of an k With k 1/L, we immediately get that the MPC error interacting uniform electron gas of r = 1, as obtained c s scales as∼ 1/L2 1/V2/3. byfollowingtheMPCansatzforvariousvaluesofN (the ∼ ∼ We now alsosee how limiting the constraintof a finite number of electrons in the simulation cell). As pointed cutoff when using an MPC like analysis actually is: The out in Ref. 4, the true interacting structure factor S is k quadratic model HF structure factor of Eq.(10) (see the quadratic at k 0.13 Fig. 3 shows that our MPC DMC → dashed green line of Fig. 2) and the actual MPC HF structure factor nicely reproduces this low-k limit (thick structure factor (solid black lines of Fig. 2) are nearly dottedblue line marked“asymptotic”). The MPC DMC identicalandquite differentfromthe exactHFstructure structure factor of Fig. 3 could be improved by using at factor (dashed-dotted red line and dotted blue line of low wavevectorsthe structure factorthat one can obtain Fig. 2). This is despite the model MPC HF structure numerically in the random-phase approximation (RPA), factor only using the asymptotic quadratic shape and which is known to be accurate in the low-k regime. 5 C. Semi-infinite electron gas givesScor.,thefinite-sizecorrectedstructurefactor. This k will differ from SMPC only between k = 0 and k = k ; k c Itiswellknownthatinthe caseofasemi-infinite elec- the fraction in Eq. (18) is therefore well-defined. But trongasthecorrectbehaviorofthestructurefactorS at the integral of the kernel, (2/π) kcsin(ku)/ku=1/u, is k 0 low k consistsof a linear term accountingfor the surface precisely the implied Coulomb potential VC. Changing R contribution that is augmented by the usual quadratic the kernel thus changes the implicit Coulomb potential: bulk term. For a finite simulation cell of width L and z surface area A, one finds (in the long-wavelength limit ∆Vc = 2 kcdk Skcor.−SkMPC sin(ku). (19) k 0):14 π SMPC 1 ku → Z0 (cid:18) k − (cid:19) S =αk+βk2, (15) This is a short-range correction in the sense that when k performing the MC run u is effectively restricted to where now within the simulation cell. On the other hand, the k π in the integral is smaller than k , which is inversely pro- α= [1/ω 1/2ω ] (16) c p s L − portional to the simulation-cell size, and so ku entering z the sin function in Eq. (19) is generally below 1. and In the case of the model of Eqs. (10) and (11) we get, β =1/2ω . (17) for example, p Hfroerme,Eβqsr.e(p2r.e1s3e)natsndth(e3.u3s4u)aolfbRuelkf.1te4r,mω,13=aωnd/α√2foblleoiwngs ∆Vc = 2 kcdk α[1− k2kf + kkf22] sin(ku). (20) the surface-plasmon energy. As the sysstempgets larger πuZ0 −1+α2kkf2 −αkkf23  (Lz ) the relative surface contribution shrinks, as   expe→cted∞. Nonetheless, surface energies, which have no SkH mustthenalsobeevaluatedusing∆VC,asitisused contribution from the bulk part of the structure factor, in Sk via nxc. In contrast, UH continues to be the stan- areentirelydominatedbythe linearsurface-contribution dard Ewald Hartree energy. αk. The MPC ansatz will never yield the linear term of IV. SUMMARY AND CONCLUSIONS Eq. (15). In order to solve this shortcoming, one can perform a simple analysis similar to the one leading to Eq. (14), now employing the linear surface term of First of all, we have demonstrated that modelling the Eq. (16). This would yield an MPC error that scales as spherically averaged structure factor Sk with the MPC 1/V. However,for this to makesense one needs a well ansatz (momentum-resolved MPC) yields exactly, after ∼definedcutoffk belowwhichthequadraticMPCbehav- integration, the MPC interaction energy of Ref. 1. This c ior should be replaced by the correct linear behavior of allows us to see explicitly that in the case of solids and Eq. (15) and above which the MPC structure factor is finite systems MPC improves convergence considerably essentially exact. Nevertheless, in the case of realistic (overthemoretraditionalEwaldscheme),duetothefact calculations8 of surface energies where the surface area that for such systems MPC yields the correct quadratic needs to be varied for a fixed slab width and vice versa, behavior of Sk as k 0. We also see explicitly how → finding a well-defined cutoff k might not be possible, so the MPC ansatz breaks down in the case of the non- c that one would need to explore more complex functional interacting uniform electron gas and, in general, in the formsofthestructurefactorinvolvingtheMPCstructure case of all systems exhibiting a similar pathology, i.e., factor itself and possibly its derivatives at one or more anextendedxc hole(atsurfaces,forexample), wherethe values of k, such as k = kc. In conjunction with the leading term of Sk at k =0 is proportional to k. known surface contribution α of Eq. (16), such schemes Astheexplicitk dependenceofSk atlowwavevectors could correct significantly the existing MPC surface cal- can usually be known, one can look at QMC systems at culations. different length scales separately enabling us to analyze thexcinteractionenergyatdifferentlengthscalesandto deriveacorrectionterm. Onintegration,thistermyields D. Implicit correction to the Coulomb kernel a correction to QMC calculations that are based on the model periodic-Coulomb interaction MPC. Aninterestingobservationisthatanycorrectiontothe structurefactorimpliesacorrection∆VctotheCoulomb potential VC =1/u: For 0<k <k , SMPC is given by ACKNOWLEDGMENTS c k evaluating the kernel sin(ku)/ku entering Eq. (6). In- stead, sampling the quantity Thisworkhasbeensupportedbythe BasqueUnibert- sitate eta IkerketaSaila(GrantNo. GIC07IT36607)and Skcor.−1 sin(ku) (18) the Spanish Ministerio de Ciencia e Innovacioń (Grants SMPC 1 ku No. FIS2009-09631and No. CSD2006-53). (cid:18) k − (cid:19) 6 1 L.M.Fraser,W.M.C.Foulkes,G.Rajagopal,R.J.Needs, 9 D.PinesandP.Nozieres,TheTheoryof Quantum Liquids S.D.Kenny,andA.J.Williamson,Phys.Rev.B53,1814 (Addison-Wesley,Reading, MA, 1989). (1996); A. J. Williamson, G. Rajagopal, R. J. Needs, L. 10 Notethatonlythesphericallyaveragedxcholeentershere, M. Fraser, W. M. C Foulkes, Y. Wang, and M.-Y.Chou, which is due to the fact that the Coulomb interaction is Phys. Rev.B 55, R4851 (1997). spherically symmetrical and, hence, the angular informa- 2 S. Chiesa, D. M. Ceperley, R. M. Martin, and M. Holz- tion of the xc hole drops out even for inhomogeneous sys- mann, Phys.Rev. Lett.97, 076404 (2006). tems. 3 M. Holzmann, B. Bernu, C. Pierleoni, J. McMinis, D. M. 11 Using the HF xc-hole (i.e., the exact exchange-hole) den- Ceperley, V. Olevano, and L. Delle Site, Phys. Rev. Lett. sity entering Eq. (6) in the limit of an infinitely large su- 107, 110402 (2011); M. Holzmann, B. Bernu, V.Olevano, percell, and performing thereal andreciprocal space inte- R. M. Martin, and D. M. Ceperley, Phys. Rev. B 79, grals,yields1atk>2k andSHF =3k/4k −k3/16k3 at F k F F 041308(R) (2009). k<2k , wherek representsthemagnitude of theFermi F F 4 R. Gaudoin and J. M. Pitarke, Phys. Rev. B 75, 155105 wavevector. (2007). 12 R.GaudoinandJ.M.Pitarke,Phys.Rev.Lett.99,126406 5 N. D. Drummond, R. J. Needs, A. Sorouri, W. M. C. (2007). Foulkes, Phys. Rev.B 78, 125106 (2008). 13 Thetruestructurefactorofauniformelectrongasatk→0 6 W.M.C.Foulkes,L.Mitas,R.J.Needs,andG.Rajagopal, is S = k2/2ω ; here, ω = (4πn)1/2 represents the bulk- k p p Rev.Mod. Phys.73, 33 (2001). plasmon energy, n being theelectron density. 7 B. Wood, W. M. C. Foulkes, M. D. Towler, N. D. Drum- 14 D. Langreth and J. P. Perdew, Phys. Rev. B 15 2884 mond, J. Phys.: Condens. Matter 16, 891 (2004). (1977). 8 B.Wood,N.D.M.Hine,W.M.C.Foulkes,andP.Garcia- Gonzalez, Phys.Rev.B 76, 035403 (2007).