Levy-Lieb principle: The bridge between the electron density of Density Functional Theory and the wavefunction of Quantum Monte Carlo
aa r X i v : . [ phy s i c s . c h e m - ph ] N ov Levy-Lieb principle: The bridge between the electron density ofDensity Functional Theory and the wavefunction of QuantumMonte Carlo ∗ Luigi Delle Site † Institute for Mathematics, Freie Universit¨at BerlinArnimallee 6, D-14195 Berlin, Germany.
Abstract
The constrained-search principle introduced by Levy and Lieb, is proposed as a practical, thoughconceptually rigorous, link between Density Functional Theory (DFT) and Quantum Monte Carlo(QMC). The resulting numerical protocol realizes in practice the implicit key statement of DFT:“
Given the three dimensional electron density of the ground state of a system of N electrons withexternal potential v ( r ) it is possible to find the corresponding N -dimensional wavefunction ofground state. ”From a numerical point of view, the proposed protocol can be employed to speed up the QMCprocedure by employing DFT densities as a pre-selection criterion for the sampling of wavefunc-tions. ∗ Chemical Physics Letters, in press (2015) † Electronic address: [email protected] . INTRODUCTION The Hohenberg-Kohn (HK) formulation of DFT and the subsequent development of theKohn-Sham (KS) scheme made DFT the most popular tool for electronic structure cal-culations (see, e.g. [1]). The profound meaning of the HK theorems for many-electronsystems was further strengthened by the Levy-Lieb formulation of the problem [2, 3]. Inthe next section the relevant aspects of the Levy-Lieb formulation will be reported. Such aformulation removes limitations of the HK theorem, such as the degeneracy of the groundstate, and sets the exact correspondence between the 3-dimensional electron density and the3N-dimensional wavefunction of the ground state of a system of N electrons. In practicalcalculations, this formulation is never explicitly used and it is usually considered only aconceptual proof of validity of DFT. In this paper I present the Levy-Lieb formulation ina different perspective; I propose such a formulation as the central algorithm which con-sistently merge two complementary but distinguished, approaches, namely KS DFT andQMC. The merging occurs in such a way that the numerical efficency of DFT contributes toenhance the numerical efficency of QMC and that the poor accuracy of DFT is enhanced bythe accuracy of QMC. A similar idea, restricted to Orbital-Free DFT and to Ground StatePath Integral QMC was presented in a previous paper [4]; here the idea is extended to allQMC methods based on the explicit calculation of the ground state wavefunction (above allthe most recent characterized high accuracy of the wavefunction e.g. [5, 6]), and to the KSDFT scheme. Another (suggestive, but yet not fully proved) implication of the procedureproposed is that the Levy-Lieb principle can be interpreted as a particular case of a decodingkey of a set of 3-dimensional data into a set of 3 N -dimensional data within the more generalframework of Information Theory [7]; this interpretation may have important implicationsfor problems of “inverse chemistry” (see e.g. Ref.[8]). However in this paper I wish only tomention the possible “decoding” interpretation and a detailed discussion can be found inRef.[7]. 2 I. LEVY-LIEB FORMULATION AND THE EXCHANGE-CORRELATIONFUNCTIONAL
The Levy-Lieb constrained-search principle is formulated as: E gs = M in ρ (cid:20) M in ψ → ρ D ψ (cid:12)(cid:12)(cid:12) b T + b V ee (cid:12)(cid:12)(cid:12)E + Z v ( r ) ρ ( r ) d r (cid:21) . (1)with E gs the ground state energy, ψ ( r ...... r N ) an antisymmetric wavefunction of an N -electron system, b T the kinetic operator, b V ee the electron-electron Coulomb potential operatorand v ( r ) the external potential (usually electron-nucleus Coulomb potential). The innerminimization is restricted to all antisymmetric wavefunctions ψ leading to ρ ( r ), i.e. ρ ( r ) = N R Ω N − ψ ∗ ( r , r , ..... r N ) ψ ( r , r , ..... r N ) d r ....d r N (where Ω N − is the N − ψ and searches over all the ρ ’s whichintegrate to N . Within this framework, the universal functional of the Hohenberg-Kohnformulation of DFT is written as: F [ ρ ] = M in ψ → ρ D ψ (cid:12)(cid:12)(cid:12) b T + b V ee (cid:12)(cid:12)(cid:12) ψ E . (2)The universal functional F [ ρ ] is unknown, however the accuracy of its approximations is thekey to the accuracy of DFT calculations. More specifically, within the KS formulation ofDFT, that is the most widely used approach of DFT for electronic structure calculations, itis not the entire F [ ρ ] that needs to be approximated, but only the so called exchange andcorrelation part of it, E xc [ ρ ]. Consistently with the definition of F [ ρ ], E xc [ ρ ] can be definedas follows: E xc [ ρ ] = F [ ρ ] − T nint [ ρ ] − E H [ ρ ] (3)where T nint [ ρ ] is the kinetic functional for non-interacting particles (whose explicit expressionis also unknown, see e.g. [9])and E H [ ρ ] = R ρ ( r ) ρ ( r ′ ) | r − r ′ | d r d r ′ is the Hartree energy. III. SPEED UP THE INITIAL STAGE OF QMC BY SAMPLING ψ AT GIVEN ρ Let us suppose we employ a given expression of E xc [ ρ ] (let us refer to it as: E xc [ ρ ] = R ρ ( r ) ǫ xc ( r ) d r ) in a DFT-KS calculation with external potential v ( r ). The KS calculationdone with E xc [ ρ ] will deliver a density of ground state ρ ( r ). In parallel, let us considerthe problem: M in ψ D ψ (cid:12)(cid:12)(cid:12) b T + b V ee (cid:12)(cid:12)(cid:12) ψ E within a QMC approach. The search for ψ min , in this3ase, would correspond to find the ground state wavefunction of a gas of N electrons, that iswe sample the space of antisymmetric ψ ’s, and a wavefunction is preferred to another if itsenergy is lower in absence of any other external constraints. However, if, within the QMCprocedure, we add the constraint that ψ must integrate to a given ρ ( r ), (e.g. to ρ ( r )), thenwe numerically realize the Levy-Lieb principle, that is we practically realize the operationof M in ψ → ρ of Eq.1.The relevant, and non trivial, consequence of the argument above is that, if E xc [ ρ ] was exact,then ρ ( r ) would be exact and the ψ delivered by QMC as a solution of M in ψ D ψ (cid:12)(cid:12)(cid:12) b T + b V ee (cid:12)(cid:12)(cid:12) ψ E with the constraints that ψ integrates to ρ ( r ), would correspond to the exact (within theaccuracy of QMC) wavefunction of the ground state for a system of N electrons with externalpotential v ( r ). Based on the considerations above, the proposed procedure to speed up QMC,consists of the following steps: • Employ E xc [ ρ ] for a KS calculation and obtain ρ ( r ). • Perform a QMC calculation where ψ is sampled in such a way that it integrates to ρ ( r ) and minimizes D ψ (cid:12)(cid:12)(cid:12) b T + b V ee (cid:12)(cid:12)(cid:12) ψ E . Let us call the corresponding wavefunction ψ ρ . • If the DFT functional employed is accurate enough, then ψ ρ would represent a verygood initial guess for solving the full problem in QMC.The essence of the scheme is that one, computationally trivial, DFT calculation can always be used to simplify a complex QMC calculation, at least in its initial stage. The practicalquestion concerns the construction of an efficient numerical algorithm for sampling ψ atgiven ρ ( r ) (suggestions for a specific case are given in Ref.[10]). Here for efficiency ismeant that the computational cost for the convergence of the sampling at given ρ ( r ) mustbe negligible compared to the computational cost of convergence of the standard QMCprocedure for searching ψ . The suggested procedure is well defined in each of its step andits success/failure depends only on the numerical efficiency in sampling ψ at given ρ ( r ).However, one may aim to extend the procedure beyond its utility for the “initial stage”of a QMC calculation and ask whether it may be possible to design a combined KS DFT-QMC scheme to iteratively solve the ground state problem; in the next section I propose aniterative scheme. 4 V. ITERATIVE DFT-QMC PROCEDURE: SELF-CONSISTENT DERIVATIONOF THE GROUND STATE
The key issue for the efficiency of the scheme suggested in the previous section is theaccuracy of E xc [ ρ ], the more accurate is E xc [ ρ ] the closer ψ ρ from QMC is to the realground state wavefunction. In general E xc [ ρ ] may not be enough accurate and in this sectionI suggest a procedure that updates E xc [ ρ ] and ψ in iterative manner. Interestingly, uponconvergence, the procedure goes beyond the “initial stage” of a QMC calculation and leads,in principle, to the ground state wavefunction. For this purpose I define: F [ ρ ] = Z ρ ( r ) f ( r ) d r . (4)Next I need to introduce a “universal” kinetic and Hartree energy reference for non inter-acting electrons. When this reference energy is subtracted to F [ ρ ] one obtains the exchangeand correlation energy as in Eq.3. E H [ ρ ] = Z ρ ( r ) e H ( r ) d r ; e H ( r ) = Z ρ ( r ′ ) | r − r ′ | d r ′ (5) T nint [ ρ ] = Z ρ ( r ) t nint ( r ) d r (6)where t nint ( r ) = t HF ( r ) (7) HF stays for Hartree-Fock, t HF ( r ) = K HF ( r ) ρ ( r ) , and K HF ( r ) is such that the total kineticenergy writes, T HF = R K HF ( r ) d r . This means that the kinetic energy density for noninteracting electrons is determined by the Hartree-Fock calculation of the N -electron systemconsidered. For consistency with the definition of t nint ( r ) we define: E H [ ρ ] = Z ρ ( r ) e HFH ( r ) d r (8)that is the Hartree energy density is calculated by the Hartree-Fock method for the specific N -electron system considered. In this way we consider t nint ( r ) and e HFH ( r ) as referencequantities in the limit in which electron correlations, due to the explicit electron-electroninteraction, are absent. The reason for calculating e H ( r ) and t nint ( r ) with HF is due to theneed of having a general reference energy for non interacting particles. Such a quantity shallbe universal and independent of a specific ρ ( r ) generated by a specific choice of a E xc [ ρ ],5therwise the iterative procedure shown below would be biased by the specific choice ofthe initial E xc [ ρ ]. It must be noticed that the exchange term of the HF calculation is notconsidered and used since this latter is automatically calculated in f ( r ) within the QMCapproach. The HF calculation is not needed for the first step of the procedure (“initial stage”of the QMC calculation) outlined in the previous section, but, as it is reported next, it isrequired to update/improve E xc [ ρ ] at the successive steps. In any case the HF calculation isneeded only once and is numerically not expensive. The iterative procedure consists of thefollowing steps (pictorially illustrated in Fig.1): • Employ a chosen E xc [ ρ ] for a KS calculation and obtain ρ ( r ). • Perform a QMC calculation where ψ is sampled via minimization of D ψ (cid:12)(cid:12)(cid:12) b T + b V ee (cid:12)(cid:12)(cid:12) ψ E with the additional request that ψ integrates to ρ ( r ), let us indicate the resultingwavefunction as ψ ρ . • From ψ ρ we calculate f ( r ): f ( r ) = 1 ρ ( r ) Z D ψ ρ (cid:12)(cid:12)(cid:12) b T + b V ee (cid:12)(cid:12)(cid:12) ψ ρ E d r ....d r N (9) • At this point we can numerically update the exchange and correlation energy density: ǫ xc ( r ) = f ( r ) − e H ( r ) − t nint ( r ) (10)However rather than ǫ xc ( r ) what is needed is the exchange and correlation potential: v cx = δE xc [ ρ ( r )] δρ ( r ) . In order to have v xc , one possibility is that of fitting numerically ǫ xc ( r )in terms of ρ ( r ). The numerical fit delivers an analytical expression for ǫ xc ( ρ ( r )).The analytic expression is taken as an approximation of ǫ xc ( ρ ( r )) and can be employedin the KS procedure to derive (analytically) v xc . The fitting procedure could be per-formed, for example, with approaches similar to the “kernel ridge regression” (KRR) amethod for non-linear regression that prevents over-fitting [11] and that is successfullyemployed in machine learning approaches (see e.g.[12–15]). • We can now employ ǫ xc ( ρ ( r )) (or better the corresponding v xc ) for a KS calculationwhich will deliver an updated electron density ρ ( r ).Next, starting from ρ ( r ) the scheme above can be repeated and continued until certaincriterion of convergence on ρ ( r ), chosen a priori, is satisfied.6 [ ρ ] i.e. ε ( r ) xcxc0 0QMCDFT−KS CalculationsCalculationsMin ψ < ψ | T+V ee | ψ > ψ integrates to ρ ρ ψ Calculate E xc [ ρ ] E xc [ ρ ] WITHDFT−KS Calculations ρ { reiterate ρ | _ ρ | δ | ρ _ ρ | <_ δ end Ground State ρ & ψ ρ (1)(2) (3) (4)> { (5) (6a)(6b) WITH ρ ρ ε ( xc ) r FIG. 1: The schematic representation of the iterative procedure. • Once the procedure converges, the QMC wavefunction obtained at the final i − th iteration, ψ ρ i , corresponds to the “true” ground state wavefunction of the system.Finally it must be reported that here the convergence of the iterative procedure, has notbeen proved mathematically. However the idea remains valid as long as a check that the totalenergy decreases at each iterative step is performed. As long as the total energy decreases,the DFT-QMC procedure is useful and at any successive step brings us closer to the exactground state. If the energy increases at the i -th iteration step, then, the wavefunction, ψ i − ,of the i − V. CONCLUSION
We have proposed a scheme to link in a consistent way KS-DFT and QMC via theLevy-Lieb principle. The electron density of KS-DFT, that is a quantity solution of a 3-dimensional problem can be used to significantly restrict the 3 N -dimensional space of wave-functions sampled by the QMC procedure. The procedure can be divided in two parts: first,given any exchange and correlation functional, the corresponding electron density, solutionof the KS equation, can be used as a pre-selection criterion for sampling a wavefunction thatcan be employed as initial guess in a QMC calculation. One could stop here and by em-ploying such initial guess proceed with the standard full QMC procedure. The second part7f the scheme proposed is a continuation of the first part, that is a procedure to iterativelysolve the ground state problem by combining at each step KS calculations and subsequentQMC sampling. This second part requires an additional HF calculation, in order to define ageneral energy of reference for non interacting electrons, and a numerical fitting procedure.It must be underlined that the whole procedure is based on the assumption that there existsan efficient numerical procedure of sampling ψ at a given electron density ρ ; despite someproposals were put forward in previous work, up to now there is no numerical implementa-tion of such a scheme. Moreover the application of the KRR method for the fitting proceduremay require a considerable amount of computational resources, thus its application may alsoneed to be optimized so that the total cost of the proposed scheme is convenient comparedto that of a full QMC calculation. Thus this work must be considered as a protocol thatmay be worth to explore further and eventually implement in a computational code. VI. ACKNOWLEDGMENTS
I would like to thank Luca Ghiringhelli for a critical reading of the manuscript, for theclarifying lecture about the KRR method and for several valuable suggestions. I would alsolike to thank Anatole von Lilienfeld for providing references for the KRR method and itsapplications. This work was supported by the Deutsche Forschungsgemeinschaft (DFG), viathe Heisenberg Grant, No. DE 1140/5-2. [1] W.Yang and R.G.Parr,
Density Functional Theory of Atoms and Molecules , Oxford UniversityPress, New York, 1989[2] M.Levy, Proc.Natl.Acad.Sci.U.S.A. , 6062 (1979); see also: M.Levy, Phys.Rev.A , 1200(1982)[3] E.Lieb, Int. Jour. Quant. Chem. , 243-277 (1983). An expanded version appears in DensityFunctional Methods in Physics, R. Dreizler and J. da Providencia eds., Plenum Nato ASISeries , 31-80 (1985)[4] L.Delle Site, L.M.Ghiringhelli and D.M.Ceperley, Int.J.Quant.Chem. , 155 (2013)[5] M.Holzmann, D.M.Ceperley, C.Pierleoni and K.Esler, Phys.Rev.E , 046707 (2003)
6] G.H. Booth, A.J.W. Thom and A.Alavi, J. Chem. Phys. , 054106, (2009)[7] L.Delle Site, Int.J.Quant.Chem. DOI:10.1002/qua.24823 (2015).[8] T.Weymuth and M.Reiher, Int.J.Quant.Chem. , 823-837 (2014).[9] S.B. Trickey, V.V. Karasiev, and A. Vela, Phys. Rev. B, , 075146 (2011)[10] L.Delle Site, Levy-Lieb Principle meets Quantum Monte Carlo , pg 361-375 in
Many-ElectronApproaches in Physics, Chemistry and Mathematics: A Multidisciplinary View ,V.Bach and L.Delle Site Eds. Studies in Mathematical Physics, Springer International Pub-lishing Switzerland, 2014[11] T.hastie, R.Tibshirani and J.Friedman,
Elements of Statistical Learning, Data Mining, Infer-ence and Prediction , 2nd ed. Springer, New York, 2009.[12] J.C. Snyder, M.Rupp, K.Hansen, K-R. M¨uller, and K.Burke, Phys. Rev. Lett. , 253002(2012).[13] J.C. Snyder, M.Rupp, K.Hansen, L.Blooston, K-R. M¨uller, and K.Burke,J. Chem. Phys. ,224104 (2013)[14] M.Rupp, A.Tkatchenko, K-R. M¨uller, and O. Anatole von Lilienfeld, Phys. Rev. Lett. ,058301 (2012).[15] A.Lopez-Bezanilla and O. Anatole von Lilienfeld, Phys. Rev. B , 235411 (2014)., 235411 (2014).