A Full Configuration Interaction Perspective on the Homogeneous Electron Gas
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p A Full Configuration Interaction Perspective on the Homogeneous Electron Gas
James J. Shepherd, ∗ George Booth, Andreas Gr¨uneis, and Ali Alavi † University of Cambridge, Chemistry Department, Lensfield Road, Cambridge CB2 1EW, U. K. (Dated: August 6, 2018)Highly accurate results for the homogeneous electron gas (HEG) have only been achieved to datewithin a diffusion Monte Carlo (DMC) framework. Here, we introduce a newly developed stochastictechnique, Full Configuration Interaction Quantum Monte Carlo (FCIQMC), which samples theexact wavefunction expanded in plane wave Slater determinants. Despite the introduction of a basisset incompleteness error, we obtain a finite-basis energy which is significantly, and variationally lowerthan any previously published work for the 54-electron HEG at r s = 0.5 a.u., in a Hilbert space of10 Slater determinants. At this value of r s , as well as of 1.0 a.u., we remove the remaining basisset incompleteness error by extrapolation, yielding results comparable or better than state-of-the-artDMC backflow energies. In doing so, we demonstrate that it is possible to yield highly accurateresults with the FCIQMC method in sizable periodic systems. PACS numbers: 71.10.Ca,31.15.V-, 71.10.-w, 71.15.-m
The homogeneous electron gas (HEG), described bya simple model Hamiltonian, encapsulates many of thedifficulties with modern electronic structure theory. Todate the only truly successful ab initio methods to yieldaccurate ground state energies at a range of densitieshave been quantum Monte Carlo techniques, in partic-ular diffusion Monte Carlo (DMC) [1–6]. The most fa-mous of these was the results of Ceperley and Alder fromwhich the LDA functionals of Density Functional Theorywere parameterised [1, 7]. Diffusion Monte Carlo wouldbe an exact technique but for the fixed-node approxima-tion, which requires the nodes in the wavefunction dueto fermionic exchange to be specified in advance by sometrial wavefunction. In general, the fixed-node approxima-tion lacks a method of being systematically improved tofind the exact result. Attempts to go beyond the fixed-node approximation have been met with some success,however complete elimination of this error has not beenachieved[1, 2, 4, 6]. In particular, the release node (RN)method is practical only in systems for which the Bosonicground state is close in energy to the Fermionic one. Inthe HEG this is only true at low density. At high den-sities, the RN-DMC is unstable, and fixed-node DMCwith backflow corrections is the most viable option. Thisleaves open the question of the magnitude of the remain-ing fixed-node error.Full configuration interaction (FCI) aims to find thewavefunction expressed as a linear combination of Slaterdeterminants, formed from rearrangements of N elec-trons in an underlying one-electron basis of M spatialorbitals[8, 9]. This is equivalent to an exact diagonaliza-tion of this space. Since such a basis set of Slater deter-minants scales as (cid:0) MN/ (cid:1) , benchmarks from FCI are ex-tremely difficult to produce. There has been surprisingly ∗ Electronic address: [email protected] † Electronic address: [email protected] little work undertaken with polynomially-scaling high-accuracy quantum chemical techniques, even though ithas recently been shown that finite systems ranging fromas few as 54 electrons can begin to capture the physics ofthe 3D HEG accurately[10, 11]. In part this might be dueto the required size of the one-electron basis and that, onapproaching the thermodynamic limit for metals, manyapproximate methods find divergent energies[12]. In con-trast, truncated configuration interaction will tend to-wards zero correlation energy.We present the application of a new method, FCI quan-tum Monte Carlo[13–16], which stochastically samplesthe exact wavefunction providing the accuracy of ex-act diagonalization at a greatly reduced computationalcost, to the high-density 54-electron HEG at r s =0.5 and1.0 a.u. This is the regime in which backflow correctionsto FN-DMC are the largest[4, 6]. The Model.-
We seek to find the ground-state wave-function and total energy of the N -electron HEGsimulation-cell Hamiltonian:ˆ H = X α − ∇ α + X α = β
12 ˆ v αβ + 12 N v M (1)where the two-electron operator ˆ v αβ is the Ewald inter-action,ˆ v αβ = 1Ω X q v q e i q · ( r α − r β ) ; v q = (cid:26) π q , q = , q = (2) v M is the Madelung term, which represents contribu-tions to the one-particle energy from interactions betweena point charge and its own images and a neutralisingbackground[10, 17, 18], and Ω is the real-space unit cellvolume.We use an expansion of the wavefunction in a Slaterdeterminant basis, Ψ = X i C i | D i i , (3)where each determinant is a normalised, antisymmetrizedproduct of plane waves, D i = A [ ψ i ( x ) ψ j ( x ) ...ψ k ( x N )] (4) ψ j ( x ) ≡ ψ j ( r , σ ) = r e i k j · r δ σ j ,σ . (5)The i index, which uniquely labels each determinant, isits normal-ordered string[19]. The wavevectors k are cho-sen to correspond to the reciprocal lattice vectors of areal-space cubic cell of length L , k = 2 πL ( n, m, l ) , (6)where n , m and l are integers. The Hartree-Fock determi-nant is the determinant occupying N plane waves withthe lowest kinetic energy. The full basis set for our calcu-lation is constructed of all Slater determinants that canbe made from M plane waves (2 M spin orbitals) forminga closed-shell of orbitals in k-space with a kinetic energylower than an energy cutoff k c . Plane waves are conve-nient because taking a single cutoff parameter to infinitymakes the one-electron basis set complete. Moreover,plane waves are natural orbitals for the electron gas, im-plying that a FCI expansion is rapidly convergent in thisbasis[20].The determinant expansion given in Eq. 3 can beinserted into the imaginary-time Schr¨odinger equation,yielding a set of coupled equations for the determinantcoefficients − dC i dτ = ( H ii − S ) C i + X j = i H ij C j . (7)Setting dC i /dτ = 0 and solving for S by exact diagonal-ization yields the total energy for the problem in a givenbasis.In a novel and recently developed quantum MonteCarlo, termed Full Configuration Interaction QMC(FCIQMC)[13], Eq. 7 is regarded as a set of master equa-tions governing the dynamics of the evolution of the de-terminant coefficients in imaginary time, with elementsof H being non-unitary transition rates. These dynamicsare simulated by introducing a population of N w ‘walk-ers’ distributed over the determinants, which are signedto represent the sign of the coefficients within the simu-lation, C i ∝ h N i ( τ ) i . The walker population is then al-lowed to evolve through discretized imaginary time-stepsby spawning, death/cloning and annihilation events ac-cording to Eq. 7 until a steady-state is reached. Theexact rules for this can be found in [13].The parameter S , termed the shift, is a population con-trol parameter which can be updated self-consistently atequilibrium to oscillate around the total energy. How-ever, throughout this work, the projected energy is used (cid:0) ) / a.u.3.213.223.233.243.253.263.27 T o t a l e n e r g y p e r e l e c t r o n / a . u . E i (cid:1) FCIQMC ( (cid:2) ) (a) A typical i -FCIQMC run. At τ ≃ . S wasallowed to vary to keep the walker number at an average of 20million. From this point, an average was taken of the total energy.(b) i -FCIQMC calculations with n add =3 are run at increasing N w values, with the aim that the limit N w → ∞ is found by thesimulation at maximum walker number. FIG. 1: Plots showing calculation of the i -FCIQMC energy for N = 54, 2 M = 682, r s = 0.5 a.u.The result is reached in the limit of long time (iterationnumber) and large walker number N w .as a stochastic correlation energy estimator, E FCIQMC = h E ( τ ) i = X j h D j | H | D i h N j ( τ ) ih N ( τ ) i , (8)where D is taken as the Hartree-Fock determinant andthe sum j need only be taken over the O (cid:2) N M (cid:3) doublyexcited determinants of D .Typically the system is initially grown by setting S tosome positive value and allowing evolution from a singledeterminant to allow an unbiased evolution of the popula-tion. Only populations above a critical system-dependentsize are able to converge to the FCI distribution, and thissize scales linearly with the size of the Hilbert space[13].In order to alleviate this problem, an adaptation of thismethod has been developed, called initiator-FCIQMC( i -FCIQMC)[14–16]. The determinant space is instanta-neously divided into those determinants exceeding a pop-ulation of n add walkers, termed initiator determinants,and those that do not. When considering a determinantwhose current population is zero, the sum in the second E rr o r ( E i (cid:3) F C I Q M C ( N w ) - E i (cid:4) F C I Q M C ( (cid:5) ) ) / a . u .
162 spinorbitals246 spinorbitals514 spinorbitals2378 spinorbitals
FIG. 2: Convergence with walker number up toapproximately 100 million walkers (taken to be theinfinite walker limit) is shown for a variety of basis setsfor N = 54, r s = 1.0 a.u. Each line is labelled with thespinorbital number, 2 M and was calculated with n add =3. As the basis set size grows, so the size of spaceand the number of walkers required to sample the spaceaccurately grows.term of Eq. 7, the term describing net flux of walkersonto that determinant, is taken to be only over initiatordeterminants. This effectively introduces a survival ofthe fittest criterion for survival of newly spawned walk-ers. If a walker has been spawned from a determinantwith an instantaneous population exceeding a parame-ter n add , the child is allowed to survive. However, ifthe parent walker is a determinant with a populationsmaller or equal to n add then the child only survives if ithas been spawned to a currently occupied determinant.This i -FCIQMC has been shown to dramatically acceler-ate the convergence of FCIQMC with respect to walkernumber. Note that in the large walker number limit, the i -FCIQMC tends to the FCIQMC algorithm, which it-self converges rigourously to the FCI energy. Figure 1illustrates an i -FCIQMC energy calculation in this way.Previous work has shown that the rapid convergence tothis limit can be examined by finding correlation energiesat increasing walker numbers (Fig. 1b)[16].As the basis set size grows, so the number of walkersrequired to recover the total energy to a given level ofaccuracy increases (Fig. 2). Further results are taken tobe converged with respect to this initiator error. Basis set extrapolation.-
Although i -FCIQMC is ableto produce exact results in a finite basis set, these are onlyupper-bound estimates to the true ground-state energy.This error in the energy, termed the basis set incomplete-ness error, is absent in DMC results which do not have a r s FCI 2 M = 2838 FCI M = ∞ / a.u. / a.u. per electron / a.u. per electron0.5 3.22086(2) 3.2202(2)1.0 0.53073(4) 0.5300(3) TABLE I: i -FCIQMC total energies for 2 M spinorbitals. The error estimate for the finite-basiscorresponds to stochastic error. The M = ∞ result isbased on extrapolations shown in Fig. 3, from whichthe error estimate derives.substantial dependence on a basis set[21].Extrapolation of the correlation energy to the com-plete basis set limit is performed regularly in molecu-lar systems for which scaling laws have been investi-gated extensively[22]. In plane wave systems, a 1 /M extrapolation is used for the basis set incompleteness er-ror in methods employing the random phase approxima-tion and second-order Møller-Plesset theory[23, 24]. Wenote that analytic expressions can be derived with thesemethods for the HEG that also show a 1 /M relation-ship. Figure 3 illustrates that by using this fit at highbasis set sizes, complete basis set exact diagonalizationenergies can be obtained that compare well with mostrecent high-accuracy DMC results[6]. The unquantifiedinitiator error is sub-mE h , and is the largest source oferror in these results but these are nonetheless thoughtto be upper-bound estimates of the exact energy, due tothe observed variationality of the initiator error in largebasis sets (see Fig. 2). Results and Conclusions.-
Results of i -FCIQMC calculations performed on the 54-electrongas, for basis sets containing between 162 and 2838spinorbtials, are shown in Fig. 3 and Table I. In thesecalculations, Hilbert spaces ranging from 10 to 10 Slater determinants (Fig. 3) were sampled to producehigh-accuracy energies for the HEG using approximately300,000 core hours in total.For r s =1.0 a.u., we obtain a variational finite-basisresult that lies just below the fixed-node DMC result,with the extrapolated energy agreeing well with back-flow DMC energies. At this r s , the remaining fixed nodeerror can be estimated by variance extrapolation of thebackflow DMC results, and is thought to be sub-mE h [4].The results presented here, containing a similar order ofmagnitude of error, seem to substantiate this claim. Toresolve the remaining fixed-node error in the backflow re-sults, we would need to further reduce the leading sourceof error in the i -FCIQMC results, which is due to basisset incompleteness.For r s =0.5 a.u., we obtain a variational finite-basisresult that lies below the backflow DMC result. How-ever, our extrapolated energy falls significantly belowthe lowest DMC result found to date, suggesting resid-ual fixed-node error in the backflow DMC energies is of (cid:6) (cid:7) (cid:8) (cid:9) (cid:10) (2M) (cid:11) T o t a l e n e r g y p e r e l e c t r o n / a . u . DMC SJ Fixed NodeDMC SJ3 BackflowFCIQMC extrapolationFCIQMC (cid:12) (cid:13) (cid:14)(cid:15) Size of space (a) r s =0.5 (cid:16) (cid:17) (cid:18) (cid:19) (cid:20) (2M) (cid:21) T o t a l e n e r g y p e r e l e c t r o n / a . u . DMC SJ Fixed NodeDMC SJ3 BackflowFCIQMC extrapolationFCIQMC (cid:22) (cid:23) (cid:24)(cid:25) Size of space (b) r s =1.0 FIG. 3: i -FCIQMC total energies for a basis of 2 M spin orbitals. Each basis set corresponds to a kinetic energycutoff, with 2 M = 2838 corresponding to 208 Ryd at r s =0.5 a.u. and 52.1 Ryd at r s =1.0 a.u. Each calculation used40 million walkers for r s = 0 . r s = 1 . M → ∞ based on the expected form 1 /M using the data set with the largest number of walkers,shown with error bars in the inset. The DMC results, taken from R´ıos et al. [6], do not suffer from basis set error andare shown as two horizontal lines representing the mean plus and minus one standard deviation. Almost identicalbackflow results can be found for r s = 1 . et al. [4].the same order as the backflow corrections themselves.It has been commented that the wavefunctions for the3D HEG change much more at higher densities than atlower densities under backflow transformations[4]. Ourresults suggest that in spite of this greater change, thebackflow transformation still comparatively struggles at r s = 0 . Acknowledgements.-
The authors gratefully acknowl-edge Neil Drummond for many useful conversations andto Mike Towler for providing us the CASINO code forcomparison benchmark values at the early stages of thisstudy[25]. This work was supported by a grant from theDistributed European Infrastructure for SupercomputingApplications under their Extreme Computing Initiative. [1] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. , 566(1980).[2] M. Holzmann, D. M. Ceperley, C. Pierleoni, and K. Esler,Phys. Rev. E , 046707 (2003).[3] G. Ortiz and P. Ballone, Phys. Rev. B , 1391 (1994).[4] Y. Kwon, D. M. Ceperley, and R. M. Martin, Phys. Rev.B , 6800 (1998).[5] I. G. Gurtubay, R. Gaudoin, and J. M. Pitarke, J. Phys.:Condens. Matter , 065501 (2010).[6] P. L´opez R´ıos, A. Ma, N. D. Drummond, M. D. Towler,and R. J. Needs, Physical Review E , 066701 (2006).[7] J. P. Perdew and A. Zunger, Phys. Rev. B , 5048(1981).[8] J. Olsen, B. Roos, P. Jørgensen, and H. Jensen, J. Chem.Phys. , 2185 (1988).[9] P. J. Knowles and N. C. Handy, Chem. Phys. Lett. ,315 (1984).[10] N. D. Drummond, R. J. Needs, A. Sorouri, and W. M. C.Foulkes, Physical Review B , 125106 (2008).[11] M. Holzmann, B. Bernu, C. Pierleoni, J. McMinis, D. M.Ceperley, V. Olevano, and L. Delle Site, Phys. Rev. Lett. , 110402 (2011).[12] L. Onsager, L. Mittag, and M. J. Stephen, Annalen derPhysik , 71 (1966).[13] G. H. Booth, A. J. W. Thom, and A. Alavi, J. Chem.Phys. (2009).[14] D. Cleland, G. H. Booth, and A. Alavi, J. Chem. Phys. , 041103 (2010).[15] D. Cleland, G. Booth, and A. Alavi, J. Chem. Phys. , , 084104 (2011).[17] P. Ewald, Ann. Phys. (1921).[18] L. M. Fraser, W. M. C. Foulkes, G. Rajagopal, R. J.Needs, S. D. Kenny, and A. J. Williamson, Phys. Rev. B , 1814 (1996).[19] W. Kutzelnigg and D. Mukherjee, J. Chem. Phys. ,432 (1997).[20] E. Davidson, Rev. Mod. Phys , 451 (1972).[21] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra- jagopal, Rev. Mod. Phys. , 33 (2001).[22] W. Kutzelnigg and J. D. Morgan, J. Chem. Phys ,4484 (1992).[23] J. Harl and G. Kresse, Physical Review B , 045136(2008).[24] A. Gr¨uneis, M. Marsman, and G. Kresse, J. Chem. Phys. , 074107 (2010).[25] R. J. Needs, M. D. Towler, N. D. Drummond, and P. L.R´ıos, J. Phys.: Condens. Matter22