5 0 0 2 Computational methods in Coupled Ele tron-Ion n a J Monte Carlo 4 ] a ∗ b h Carlo Pierleoni and David M. Ceperley p - p m 9th February 2008 o c . a s c Department of Physi s, University of L'Aquila, Via Vetoio, I-67010 L'Aquila, i s Italy y b h Dept. of Physi s and NCSA, University of Illinois at Urbana-Champaign, Ur- p [ bana, IL 61801, USA 1 v Abstra t 3 1 0 In the last few years we have been developing a Monte Carlo simula- 1 0 tion method to ope with systems of many ele trons and ions in the Born- 5 Oppenheimer (BO) approximation, the Coupled Ele tron-Ion Monte Carlo 0 / s Method (CEIMC). Ele troni properties in CEIMC are omputed by Quan- c i tum Monte Carlo (QMC) rather than by Density Fun tional Theory (DFT) s y based te hniques. CEIMC an, in prin iple, over ome some of the limita- h p tions of the present DFT based ab initio dynami al methods. Appli ation of : v thenewmethodtohighpressuremetalli hydrogenhasre entlyappeared. In i X thispaperwepresentanewsamplingalgorithmthatwehavedevelopedinthe r a framework of the Reptation QuantumMonte Carlo (RQMC) method hosen to sample the ele troni degrees of freedom, thereby improving its e(cid:30) ien y. Moreover,weshowherethat,atleastforthe aseofmetalli hydrogen,varia- tionalestimatesof theele troni energies leadto ana uratesampling ofthe protondegrees of freedom. ∗ Corresponding author: Carlo Pierleoni, Physi s Dept. University of L'Aquila, Via Vetoio, 67010Coppito,L'Aquila(Italy),fax:+39-0862433033;email: arlo.pierleoniaquila.infn.it 1 1 Introdu tion Modern theoreti al methods in ondensed matter physi s and hemistry rely heav- ily on numeri al simulations. The problem of solving the S hroedinger equation for many-body systems is too di(cid:30) ult to be addressed dire tly, even within the simpli(cid:28) ationprovidedbytheBorn-Oppenheimerapproximation. Inthemostpop- ular pra ti al approa hes (Hartree-Fo k (HF) and the Density Fun tional Theory (DFT)basedmethods[1℄)theoriginalproblemisrepla edbytheproblemofsolving the time independent S hroedinger equation for a single ele tron in the (cid:28)eld of the nu leiandthemean (cid:28)eldgeneratedbytheotherele trons. DFT is, in prin iple,an exa ttheory but the energyfun tional mustbe treated approximatelyforpra ti al purposes. In the simplest Lo al Density Approximation (LDA), this exa t theory be omes a self- onsistent mean (cid:28)eld theory. Extensions of LDA, su h as Gener- alized Gradient Approximation (GGA) provide more a urate results but remain essentially at the level of an e(cid:27)e tive mean (cid:28)eld treatment. Despite the mean (cid:28)eld hara ter, DFT s hemes have proved to provide quite a urate results for many di(cid:27)erent systems[1℄ In 1985, Car and Parrinello introdu ed an e(cid:30) ient method to ouple standard Mole ular Dynami s for lassi al nu lei with the ele troni stru ture al ulation at the level of LDA done (cid:16)on the (cid:29)y(cid:17) to extra t the nu lear for es[2℄. Be ause the method allowed study of the statisti al me hani s of lassi al nu lei with many body ele troni intera tions, it opened the way for the use of simulation methods for realisti systemswith an a ura ywell beyond the limits ofe(cid:27)e tive for e (cid:28)elds available. Inthelasttwentyyears,thenumberofappli ationsoftheCar-Parrinello ab-initio mole ular dynami s has ranged from simple ovalent bonded solids, to highpressurephysi s,materials ien eandbiologi alsystems. Therehavealsobeen extensionsoftheoriginalalgorithmtosimulatesystemsat onstanttemperatureand onstant pressure[3℄, (cid:28)nite temperature e(cid:27)e ts for the ele trons [4℄, and quantum nu lei [5℄. Despite re ent progress, DFT su(cid:27)ers from well-known limitations, for example, ex ited state properties su h as opti al gap and spe tra are less reliable. DFT showsseriousde(cid:28) ien iesin des ribingvan derWaalsintera tions, non-equilibrium geometries su h as rea tion barriers, systems with transition metals and/or luster isomerswith ompetingbonding patterns[1,6℄. Asa onsequen e, urrentab-initio predi tionsofmetallizationtransitionathighpressures,orevenpredi tionofphase transitionsareoftenonlyqualitative. Hydrogenisanextreme ase[7,8,9℄ but even β in sili on the diamond/ -tin transition pressure and the melting temperature are seriously underestimated[10℄. Another route to the ground state properties of a system of many ele trons in presen e of nu lei is the Quantum Monte Carlo method[11, 6℄. In its simplest form, an analyti many ele tron wave fun tion is hosen on the basis of the vari- ational prin iple (Variational Monte Carlo, VMC) and the quantum averages are obtained by a Metropolis Monte Carlo simulation of the ele troni oordinates. A morea uraterepresentationofthe groundstatewavefun tion anbeobtained by exp{−β H} H e proje ting the variational wave fun tion with the operator where β e is the many-body hamiltonian, and is the proje tion time. Provided that the variational wave fun tion is not orthogonal to the ground state wave fun tion, the proje ted fun tion tends exponentially fast to the ground state wave fun tion as β → ∞ e . Sin e matrix elements of the above proje tion operator at large values β P e of are unknown for non trivial systems, a Trotter breakup in many ( ) small τ =β /P) e e imaginary time intervals ( must be employed. In the on(cid:28)guration rep- 3N resentation, ea h proje tion orresponds to a -dimensional integral whi h an be performed by Metropolis Monte Carlo method provided that the propagator in imaginarytime anbe hosenrealand anbe interpreted asaprobabilitydistribu- tion. This is the essen e of the Di(cid:27)usion Monte Carlo method (DMC) whi h is an (cid:16)exa t(cid:17) methodforsystemsofbosonsorboltzmannons. Thismeansthatallsystem- ati errors in a simulation are under ontrol in the sense that they an be redu ed asmu hasdesired. Sin e ele tronarefermions,the aboves heme fails be ausethe imaginary time propagator must be ompletely antisymmetri under ex hange of two ele trons and therefore annot be hosen stri tly non-negative everywhere in on(cid:28)gurational spa e. This is the origin of the infamous (cid:16)fermion sign problem(cid:17). In order to avoid the sign problem the (cid:16)(cid:28)xed node approximation(cid:17) has been pro- posed and used routinely to perform fermion simulations[6℄. The energy al ulated withthis approximationisvariationalwithrespe ttothepositionofthenodalsur- fa es of the trial wave fun tion. Over the years, the level of a ura y of the (cid:28)xed 3He nodeapproximationforsimplehomogeneoussystems,su has andtheele tron gas,hasbeensystemati allyimprovedbyintrodu ingmoresophisti atednodalsur- fa es (ba k(cid:29)ow orbitals)[12, 13℄. In more omplex, inhomogeneous situations su h as atoms, mole ules and extended systems of ele trons and nu lei, progress have been somewhat slower. Nonetheless, in most ases, (cid:28)xed-node QMC methods have provedtobemorea uratethanmean(cid:28)eldmethods(HFandDFT)[6℄. Computing ioni for eswith QMCto repla ethe DFT for esin theab-initioMD, ismoredi(cid:30)- ult and a general and e(cid:30) ient algorithm is still missing. Moreover, the omputer time requiredforaQMCestimateofthe ele troni energyis, in general,morethan for a orresponding DFT-LDA al ulation. These problems have seriously limited the development of an ab-initio simulation method based on the QMC solution of the ele troni problem (cid:16)on the (cid:29)y(cid:17). In re ent years, we have developed a di(cid:27)erent strategy based entirely on the Monte Carlo method both for solving the ele troni problem and for sampling the ioni on(cid:28)gurationspa e[14,15℄. Thenewmethod, alledtheCoupledEle tron-Ion Monte Carlo method (CEIMC) has been applied so far to high pressure metalli hydrogen where it has found quite di(cid:27)erent e(cid:27)e ts of temperature than CPMD basedontheLDAfor es[16℄. Ourpresentinterpretationofthedisagreementisthat LDA providesaBorn-Oppenheimersurfa e quitesmootherthanthe morea urate T >0 QMC one and this strongly a(cid:27)e ts the stru ture of the protoni system at . Thepaperisorganizedasfollows. Thefollowingse tion2isdevotedtoanoutline of the CEIMC method. We will not go into all details sin e two long arti les have appearedon generalaspe ts and earlyimplementationsofthe method[14,15℄. One of the new aspe ts that we havere ently implemented in CEIMC, not des ribed in thosereferen es,istheReptationQuantumMonteCarloproje tionoftheele troni variationalwavefun tion[17℄. Sointhesubse tion2.1wereviewtheRQMCmethod andinthefollowingsubse tion2.2wefo usonthesamplingalgorithm,weintrodu e our new s heme to improve e(cid:30) ien y and reliability of RQMC and we provide an analyti al proof. In Se tion3 we report numeri al results on the onvergen eof the new s heme with the proje tion time and with the Trotter time step. Finally, in se tion 4, we on lude. 2 The Coupled Ele tron-Ion Monte Carlo method CEIMC method is based on the Born-Oppenheimer (BO) separation between the slow nu lei and the fast ele trons. This is in ontrast with other Quantum Monte Carlo methods, Di(cid:27)usion Monte Carlo (DMC)[11, 6℄ or (cid:28)nite temperature Path IntegralMonte Carlo(PIMC)[18,19℄ methods where ele tronsand ions aretreated on the same footing. As usual, the BO approximation allows to over ome the limitations of the other QMC methods, while introdu ing an often negligible error. InCEIMC,the on(cid:28)gurationalspa eoftheprotondegreesoffreedomatinverse β = (k T)−1 B temperature is sampled with a Metropolis algorithm in whi h the S S′ di(cid:27)eren e between the BO energy of a proton state and of a trial state is omputed by an ele troni ground state QMC al ulation. The QMC estimate of ∆=[E(S′)−E(S)] the energydi(cid:27)eren e hasstatisti alnoisewhi hwouldbias the standard Metropolis algorithm. Unbiased sampling of the proton on(cid:28)gurations is ∆ a hieved by the penalty method[20℄ whi h repla es the energy di(cid:27)eren e in the ∆+(βσ∆)2/2 σ∆2 a eptan eformulaby ,where isthevarian eoftheenergydi(cid:27)er- σ2 > 0 ∆ en e. Sin e , the noise always auses extra reje tions but this ompensates for (cid:17)uphill(cid:17) moves a epted be ause of a favorable energy (cid:29)u tuation. Severalmethodsfor omputingenergydi(cid:27)eren esinQMCareavailable.[14,15℄. Asimpleand e(cid:30) ientmethodisto sampletheele troni degreesoffreedom from a distributionfun tionwhi histhesumoftheele troni distributionfun tionsforthe S S′ and states(e. g. thesumofthethesquareofthetrialwavefun tionsinVMC). Averages of operators involving ele troni degrees of freedom and a single proton S on(cid:28)guration, say , (for instan e total energy, varian e, et .) are then omputed by orrelatedsampling[11,14,15℄. Forthetypi alsizeoftheprotonmoves(between 0.01 0.5 Å and Å for lassi al protons depending on density and temperature) and thetypi alsystemsize(upto54protons)wehaveinvestigated,thismethodismu h moree(cid:30) ientthanperformingtwoindependentele troni al ulationsforthestate S S′ and . In the ground state QMC methods, an ele troni trial wave fun tion must be hosena ordingtothephysi softhesystembeingstudied. Forthemetalli phase of hydrogen, we have re ently developed analyti fun tions whi h in lude ba k(cid:29)ow and three-body orrelations[21℄. These wave fun tions are parti ularly appropriate to CEIMC sin e they have a urate energies already at the variational level, they havenoadjustableparametersrequiringoptimization,andtheir omputational ost is mu h less than using orbitals expanded in a plane wave basis typi ally used in QMC al ulations[22℄. Inmetalli systems,(cid:28)nitesizee(cid:27)e tsarelargeandmustbesuitablytreated. The ommon pro edure is to repeat the al ulation for systems of in reasing size and extrapolatetothethermodynami limitbutthisisimpra ti alwithinCEIMC,sin e itwouldhavetobeperformedforanyproposedprotoni stepbeforeitsa eptan e. A mu h better strategy is to use Twist Averaged Boundary Condition (TABC) 1/N [23,15℄whi hredu esthe(cid:28)nitesizeerrorintheenergytothe lassi al behavior. It onsistsinaveragingtheenergyoverthephasethatthemanybodywavefun tion an pi k if a single ele tron wraps around the super- ell. This is equivalent to Brillouin zone sampling in the single ele tron approximation. Within CEIMC it does not ause a large in rease in required CPU time/step. Finally a re ent improvement of the method is the introdu tion of quantum e(cid:27)e ts for the protons, quite important in high pressure hydrogen. This is done by developing the thermal density matrix of protoni degrees of freedom on the BO surfa e in Feynman Path Integrals[18, 16℄. A similar te hnique in the ontext of Car-Parrinello method has appeared[5℄. We are not going to dis uss the last two aspe tsoftheCEIMC.WhileTABCimplementationinCEIMChasbeendes ribed in ref.[15℄, our implementation of PIMC for proton degrees of freedom in CEIMC will be the subje t of a future publi ation. 2.1 Reptation Quantum Monte Carlo Method To go beyond VMC ele troni energies, we implemented a Reptation Quantum Monte Carlo algorithm (RQMC)[17℄, rather than Di(cid:27)usion Monte Carlo algorithm (DMC). The implementation of the energy di(cid:27)eren e method is more straightfor- ward in RQMC, nor are averages of observables whi h do not ommute with the hamiltonian biased. In RQMC the ground state wave fun tion is obtained by onstru ting an imag- |Ψ0i inary time path integral for the ele troni degrees of freedom. If is the trial state,thetrialstateproje tedina(cid:16)time(cid:17) βe/2,|Ψβe/2i=e−βeH/2|Ψ0i,will onverge β e to the ground state for large . Let us de(cid:28)ne the (cid:16)partition(cid:17) fun tion Zβe =hΨ0|e−βeH|Ψ0i. (1) The energy is then de(cid:28)ned as d 1 E(βe)=−dβ lnZβe = Z hΨ0|e−βeHH|Ψ0i=hEL(R)i (2) e βe EL(R) = ℜ(Ψ−01(R)HΨ0(R)) where the averages of the lo al energy, , are with respe t to the path average. In pra ti e, the energy is omputed as the average ℜ of the lo al energy at the two ends of the path. Here indi ates the real part in E(β ) e the ase that the trial fun tion or Hamiltonian is omplex. The energy is an β e upper bound to the (cid:28)xed node energy for ea h value of , it onverges to this at β β e e large , and its derivative is stri tly negative. This latter quantity is, in fa t, minus the varian e of the total energy dE(β ) σ2(β )=− e =hE (0)E (β )i−hE(β )i2. e L L e e dβ (3) e β e The varian e tends to zero for large enough values of providing a useful signal for the onvergen e of the energy to the ground state. This is the zero varian e theorem in RQMC. On the other hand, the varian e of the lo al energy omputed σ2(β ) e at either end of the path is the mixed estimator of DMC for . In pra ti al β e implementations, it is desirable to keep as small as possible to maximize the e(cid:30) ien y of the energy di(cid:27)eren e method. To ompute the needed density matrix elements, we divide the proje tion time β P τ = β /P e e e into time sli es and make a semi- lassi al approximation for exp(−τeH) R={r1,...,rN} . Our notation for a single ele troni on(cid:28)guration is , s = {R0,R1,...,RP} while for the entire path is . The probability distribution for a path is P−1 Π(s)=exp −U(R0)−U(Rp)− Ls(Ri+1,Ri) ( ) (4) i=1 X U(R) = ℜ[lnΨ0(R)] Ls(R,R′) where and is the symmetrized link a tion for our approximationoftheshorttimepropagator. Wehaveusedtheimportan esampling Green's fun tion of the DMC propagator R|e−τeH|R′ = Ψ0(R) exp −τ E (R)− [R′−R−2λτeF(R)]2 . (cid:12)Ψ0(R′)(cid:12) (cid:20) e L 4λτe (cid:21) (5) (cid:10) (cid:11) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 2 F(R)=∇U(R) λ=h¯ /2m e wherethefor eis and ,whi hprovidesthesymmetrized link a tion τ L (R,R′) = e E (R)+E (R′)+λ F2(R)+F2(R′) s L L 2 (R(cid:2)−R′)2 (R−R′)·(F(cid:0)(R)−F(R′)) (cid:1)(cid:3) + + 4λτ 2 (6) e An alternative form for the link a tion ould be obtained through the pair a tion developed in (cid:28)nite temperature Path Integral MC[18℄. However, we have not im- plemented this form and do not have a omparison of its e(cid:30) ien y. Inordertoimposethe(cid:28)xedphase onstraintontheproje tedwavefun tion,we LFP(R,R′)=λτ 1dη|∇φ(X(η))|2 mustaddtothelinka tionatermoftheform s e 0 φ(X) R X where isthe phase of the trial wavefun tion at ele troni position and the X(η) X(0)= R, X(1)= integral is taken over all paths with boundary onditions R′ . We have taken an end-point approximation for this term ex ept for real wave fun tions in whi h ase (cid:28)xed-node boundary onditions were used. Note that in the expressions above, the dependen e on the nu lear degrees of freedomwasnotshowneventhoughallquantitiesdependonthem. Theprobability Π(s,S) S distribution of an ele troni path will be where indi ates the position of all nu lei. Be ause we have an expli it distribution of the ele troni paths, it is straightforwardtoapplytheimportan esamplings hemefortheenergydi(cid:27)eren es [Π(s,S)+Π(s,S′)] S S′ bysamplingtheprobabilitydistribution where and arethe Π(s,S) ∝ urrent and the trial protoni state, respe tively. Note also, for VMC 2 |Ψ0(s,S)| be omes the squareofmodulus ofthetrial wavefun tion(noproje tion R0 =RP and ). 2.2 The (cid:16)boun e(cid:17) algorithm In the original work on RQMC[17℄, the ele troni path spa e was sampled by a reptation algorithm, an algorithm introdu ed to sample the on(cid:28)gurational spa e of linear polymer hains. The slithering snake or reptation method seems to have s originated by Kron[24℄ and by Wall and Mandel[25℄. Given a path on(cid:28)guration , R0 RP amoveisdoneintwostages. Firstoneofthetwoends(either or )issampled R g with probability 1/2 to be the growth end . Then a new point near the growth R +2λτ F(R ) g e g end is sampled from a Gaussian distribution with enter at . In orderto keep the number oflinks on the path length onstant, the oldtail position isdis ardedinthetrialmove. Themoveisa eptedorreje tedwiththeMetropolis formula based on the probability of a reversemove. For use in the following, let us d d=+1 R = R d=−1 g P de(cid:28)ne the dire tion variable as for a head move ( ), and Rg =R0 d foratailmove( ). Instandardreptation,thedire tion is hosenrandomly at ea h attempted step. P(s→s′) In the standardreptation algorithm, the transition probability is the T (s→s′) a (s→ d d produ tofanattemptprobability andana eptan eprobability s′) . Note that the path distribution given in Eq.(4), be ause of the symmetrized d link a tion does not depend on the dire tion in whi h it was onstru ted. In the Metropolis algorithm, the a eptan e probability for the attempted move is Π(s′)T (s′ →s) a (s→s′)=min 1, −d d Π(s)T (s→s′) (7) (cid:20) d (cid:21) P (s→s′) d whi h ensures that the transition probability satis(cid:28)es detailed balan e Π(s)P (s→s′)=Π(s′)P (s′ →s) d −d (8) The auto orrelation time of this algorithm in Monte Carlo steps, that is the [(β /τ )2/A] e e numberofMCstepsbetweentwoun orrelated on(cid:28)gurations,s alesas , A β e where is the a eptan e rate, an unfavorable s aling for large . Moreover the o asionalappearan eofpersistent on(cid:28)gurationsboun ingba kandforthwithout reallysamplingthe on(cid:28)gurationspa ehasbeenpreviouslyobserved[26℄. Theseare two very unfavorable features, parti ularly in the present ontext, where we need to perform many di(cid:27)erent ele troni al ulations (at least one per protoni move). There is a premium for a reliable, e(cid:30) ient and robust algorithm. We have found that a minimal modi(cid:28) ation of the reptation algorithm solves both of these problems. The idea is to hose randomly the growth dire tion at the beginning of the Markov hain, and reverse the dire tion upon reje tion only, the (cid:16)boun e(cid:17) algorithm. Asfarasweareaware,(cid:17)boun e(cid:17) dynami shasnotbeen previ- ously investigated for RQMC, though Wall and Mandel[25℄ mentioned it without a detailedproofandsubsequentpolymersimulationsdidnotuseboun e,perhapsbe- ausethe a eptan e ratioin the polymer systems is mu h smallerthan in RQMC. There is a related algorithmfor dire ted loop algorithm on the latti e and for sim- ulations of trapped di(cid:27)usion[28, 27℄. What follows is the proof that the boun e algorithm samples the orre t prob- Π(s) d ability distribution . The variable is no longer randomly sampled, but, as before, the appropriate move is sampled from the same Gaussian distribution T (s→s′) d and a epteda ordingtothe Eq. (7). Tobeable tousethete hniques d of Markov hains, we need to enlarge the state spa e with the dire tion variable . {s,d} In the enlarged on(cid:28)guration spa e , let us de(cid:28)ne the transition probability P(s,d → s′,d′) of the Markov hain. The algorithm is a Markov pro ess in the extended path spa e, and it is ergodi as DMC method, hen e, it must onvergeto Υ(s,d) a unique stationary state, satisfying the eigenvalue equation: Υ(s,d)P(s,d→s′,d′)=Υ(s′,d′). (9) s,d X Π(s) We show that our desired probability is solution of this equation. Within the P(s,d→s′,d′)6=0 d=d′ imposed rule not all transitions are allowed, but for and s 6= s′ d′ = −d s = s′ (a epted move), or and (reje ted move) only. Without loss d′ = +1 ±1 of generality let us assume sin e we have symmetry between . Eq. (9) Υ(s,d) Π(s) with repla ed by is Π(s′)P(s′,−1→s′,1)+ Π(s)P(s,1→s′,1)=Π(s′). s6=s′ X Π(s)P(s,1→s′,1)=Π(s′)P(s′,−1→ Be auseof detailed balan e Eq.(8), we have s,−1) , whi h when substituted in this equation gives Π(s′) P(s′,−1→s′,1)+ P(s′,−1→s,−1) =Π(s′). " # s X s s = s′ Note that we have ompleted the sum over with the term be ause its probability vanishes. The term in the bra ket exhausts all possibilities for a move