An auxiliary-Field quantum Monte Carlo perspective on the ground state of the dense uniform electron gas: An investigation with Hartree-Fock trial wavefunctions

Describing electron correlation in a scalable way that can handle hundreds of electrons is a grand challenge in quantum chemistry and condensed matter physics. State-of-the-art methods include coupled-cluster (CC) methods,^1–7 density matrix renormalization group (DMRG) methods,^8,9 and quantum Monte Carlo (QMC) approaches.^10–12 Each method exhibits different weaknesses and strengths, and therefore they have been applied to solve a different class of problems in chemistry and condensed matter physics.

In this work, we will focus on a projector QMC method, namely, the auxiliary-field QMC (AFQMC) approach.^11,13 Projector Monte Carlo methods, while formally exact, typically impose a constraint in the imaginary time propagation in order to overcome the fermion sign problem and achieve a polynomial scaling algorithm. Both, diffusion Monte Carlo (DMC)¹⁰ and AFQMC enforce this constraint using a trial wavefunction, which can, in principle, be systematically improved toward the exact result. These constraints lead to the phaseless-AFQMC (ph-AFQMC) and fixed-node (FN-DMC) algorithms both of which scale like $O(N3)−O(N4)$ with the number of electrons N.

Although the formalism of DMC and AFQMC is very similar, there are some key differences between the two. First, AFQMC works in the second-quantized framework common to most quantum chemical methods and introduces a finite basis set. Therefore, AFQMC energies need to be extrapolated to the complete basis set (CBS) limit in order to compare directly with experiments. This is in contrast with DMC, which works in real space and directly in the CBS limit. Second, incorporating widely used Jastrow factors (JFs) into AFQMC is quite challenging. JFs are economical ways to incorporate residual electron correlation by enforcing cusp conditions either between electrons and nuclei or among electrons. Finally, unlike FN-DMC, ph-AFQMC is not variational.¹⁴ 

Despite these issues, AFQMC offers a number of promising advantages precisely because it works directly in an orbital-based basis. In particular, all-electron, frozen core and nonlocal pseudopotential calculations can be performed with no additional approximations. Furthermore, as most quantum chemistry methods are performed with a finite basis set, many tricks used in quantum chemistry can be used to improve AFQMC as well. For instance, tensor hyper contraction approaches^15–20 have recently been employed to reduce the memory requirement of AFQMC.^21,22 Employing explicitly correlated basis functions (similar in spirit to JFs for DMC) should also be possible to reduce the basis set incompleteness error of AFQMC.²³ In addition to this, computing properties other than the total energy, which has historically been a challenge for projector QMC methods, can be more straightforwardly achieved in AFQMC.²⁴ Recent examples include one- and two-particle reduced density matrices,²⁴ imaginary time correlation functions,^25–27 as well as forces.²⁸ 

AFQMC has been successfully applied in recent years to a number of challenging problems in both quantum chemistry^29–33 and solid state physics.^34–36 However, the broad applicability of the method is not as well understood as more traditional quantum chemistry approaches, which have seen decades’ worth of sustained development and benchmarking. The primary limiting factor of AFQMC is the choice of trial wavefunction. Single-determinant trial wavefunctions from Hartree-Fock (HF) or density functional theory calculations have shown remarkable accuracy for a broad range of applications including the two-dimensional Hubbard model,³⁷ dipole-bound anions,³⁸ and solid state applications,^39–42 with total energies often approaching the accuracy of coupled cluster singles and doubles with perturbative triples CCSD(T).²² However, for more strongly correlated systems such as transition metal containing complexes, single determinant trial wavefunctions are not sufficiently accurate and multideterminant trial wavefunctions become necessary.^{32,33,43–46}

Often short determinantal expansions from complete active space self-consistent field (CASSCF) wavefunctions^31,33 or nonorthogonal multi-Slater deteraminant trial wavefunction³² can help to restore the accuracy of the method. However, since multideterminant wavefunctions scale exponentially with system size, this approach to improving the trial wavefunction is ultimately limited, particularly for large scale applications. The search for more economical and accurate trial wavefunctions for AFQMC (and also FN-DMC) is an active area of research, and no single approach can achieve both polynomial scaling and broadly consistent accuracy at the same time.

Our goal is to assess the quality of HF trial wavefunctions for the Uniform Electron Gas (UEG) model.⁴⁷ HF wavefunctions are the simplest possible reference states to perform subsequent correlation calculations in quantum chemistry. Often, artificial symmetry breaking caused by HF wavefunctions causes confusion in understanding electron correlation.^48,49 It is artificial because the symmetry breaking occurs due to the lack of treatment for weak correlation, not strong correlation. Artificial symmetry breaking affects the performance of subsequent correlation calculations greatly (especially for weakly correlated systems).^50–57 In such cases, it is preferred to use spin-restricted Hartree-Fock (RHF) orbitals as opposed to variationally preferred broken symmetry HF orbitals. Therefore, one has to be cautious when choosing a proper HF state for correlation calculations. Nonetheless, HF is not only the simplest but also scalable to hundreds of electrons. With HF trials (which we name as HF+ph-AFQMC), AFQMC scales strictly as $O(N3)−O(N4)$ for a fixed statistical accuracy. Therefore, it is crucial to assess the accuracy of HF+ph-AFQMC and understand the scope of it in simulating large-scale chemical and solid-state systems. After all, the gold-standard quantum chemistry method, CCSD(T),⁵⁸ is performed on top of HF states.

In the context of this paper, the UEG model provides us with a simplified version of the full ab initio Hamiltonian for a solid, essentially omitting the electron-ion interaction term and all of the material complications it entails. One can tune the magnitude of dynamic and static correlation at the Hamiltonian level using a single parameter (for a fixed number of electrons) through the dimensionless Wigner-Seitz radius r_s. For low r_s values, the electron density is high, the distance between electrons is short, and electrons are all paired. In this regime, diagrammatic resummation techniques such as random phase approximation (RPA) and low-order perturbation theory based on RHF references are highly accurate.⁵⁹ We call this regime “weakly correlated.” On the other hand, for high r_s values, electrons are unpaired and spatially well-separated and closed-shell RHF references provide qualitatively wrong picture (i.e., it does not describe the open-shell nature of the system). In this case, RPA or other low-order perturbation theory based on RHF references fails. Furthermore, the use of broken symmetry HF references does not provide accurate description either due to spin contamination. We call this regime “strongly correlated.” This tunability allows for an unambiguous comparison between the strengths and weaknesses of various methods. Moreover, there exist a number of benchmark results both for intermediate system sizes within the reach of traditional quantum chemistry approaches, as well as results for much larger system sizes, and also results extrapolated to the thermodynamic limit.

Recently, the formally exact full configuration interaction (CI) quantum Monte Carlo method^12,60 (FCIQMC) provided benchmark results for a range of densities for a 14- and 54-electron system.^61,62 These FCIQMC studies also motivated recent coupled-cluster Monte Carlo^63–65 studies on the 14-electron UEG model by Neufeld and Thom where they provided CCSD, CCSD and triples (CCSDT), CCSDT, and quadruples (CCSDTQ) with RHF reference results for a wide range of r_s and basis sets.⁶⁶ As the scope of truncated spin-restricted CC (RCC) approaches is relatively well understood, comparing HF+ph-AFQMC against these results will lead us to a better understanding of the scope of HF+ph-AFQMC. Although the ph-AFQMC has been applied to 3D UEG before to construct the KZK functional^67,68 and to small 2D UEG models,^25,26 we believe this is the first published extensive benchmarking study of the 3D UEG using ph-AFQMC.

This paper is organized as follows: (1) we briefly review the formalism of ph-AFQMC and the UEG model, (2) we analyze the basis set convergence of AFQMC in the 14-electron UEG model and compare its result against FCIQMC and CCMC, and (3) we study larger systems (54-electron and 114-electron) and discuss the AFQMC perspectives on simulating the ground state of solids.

In this section, we briefly summarize the basics of the AFQMC method and the phaseless approximation (ph) which leads to the phaseless AFQMC algorithm (ph-AFQMC). We use n_occ to denote the number of occupied molecular orbitals (MOs) and n_vir to denote the number of unoccupied MOs.

The zero-temperature AFQMC algorithm is a stochastic realization of power methods that target the lowest root of the Hamiltonian Ĥ. The algorithm is based on the following identity:

where |Ψ₀⟩ is the exact ground state and |Φ₀⟩ is an initial wavefunction satisfying ⟨Φ₀|Ψ₀⟩ ≠ 0. Although the initial wavefunction |Ψ₀⟩ can differ from the trial wavefunction |Ψ_T⟩, for the purpose of this work, we will assume |Ψ₀⟩ = |Φ_T⟩ unless mentioned otherwise. Equation (1) is implemented stochastically by repeatedly applying a propagator, exp(−ΔτĤ), to a set of random walkers until the ground state is reached. Each walker is comprised of a Slater determinant, |ψ_n(τ)⟩, and a weight $wn$(τ) such that the statistical representation of the wavefunction is given by |Ψ(τ)⟩ = ∑_n$wn$(τ)|ψ_n(τ)⟩.

In order to practically realize the projection, we first split the Hamiltonian into one-body and two-body operators (i.e., Ĥ = Ĥ₁ + Ĥ₂). For the two-body terms, we write them in the sum of squared operators

Then, we apply the Hubbard-Stratonovich⁶⁹ (HS) transformation to rewrite the imaginary-time propagator in terms of only one-body operators. With the symmetric Trotter decomposition, the propagator reads

where p(x) is the standard normal distribution, x is a vector of N_α auxiliary fields, and $B^$ is defined as

At each time step, each walker draws Gaussian random numbers to sample one instance of x and provides a sample to the HS transformation in Eq. (3). The application of a one-body operator such as Eq. (4) to a Slater determinant yields yet another single Slater determinant.^70,71

For a generic ab initio Hamiltonian, the propagators appearing in Eq. (4) will, in general, be complex and the weights of the walkers will acquire a phase that will be distributed uniformly in the complex plane in the long imaginary time limit.³⁰ This “phase problem” is analogous to the notorious fermion sign problem encountered in DMC and has no known solution, in general. The phase problem can be somewhat mitigated through mean-field subtraction⁷² (i.e., redefining $v^α′=v^α−⟨v^α⟩0$⁠) in Eq. (4), but the statistics will be eventually swamped by the phase problem. Note that mean field subtraction is essentially identical to normal-ordering $v^α$ which ensures $v^α′|Φ0=0$ for all α.

It is possible to eliminate this phase problem entirely at the sake of introducing biases into the results using the so-called phaseless approximation.¹¹ This is achieved by first performing an importance sampling transformation to the propagator such that walkers now undergo the modified propagation

where the importance function (in hybrid form) is defined as

where S_n is the overlap ratio of the nth walker

and $x¯n$ is an “optimal” force bias which is a shift to the Gaussian distribution, given as

The phaseless approximation (ph) is then defined as a modification to this importance function

where the phase θ_n(τ) is given by

The walker weights and Slater determinants are then updated as

Evidently, the phaseless approximation ensures that the walker weights remain real and nonnegative throughout the simulation and therefore removes the phase problem completely.

The mixed estimate for the local energy estimator can be computed with the generalized Green’s function (or one-particle reduced density matrix) P,

where $Cψn$ is the occupied MO coefficient of |ψ_n(τ)⟩ and $CΨT$ is the occupied MO coefficient of |Ψ_T⟩. Once the simulation has equilibrated, we will have a statistical representation of the ground state wavefunction given by

from which we can compute the mixed estimator for the energy as

where ϵ_n(τ) is the local energy of a walker,

We will see how the local energy evaluation is done specifically for the UEG model later.

Size-consistency is a property of a wavefunction for isolated systems A and B that asserts the product separability of a supersystem wavefunction (|Ψ_AB⟩ = |Ψ_A⟩|Ψ_B⟩) and also the additive separability of energy (E_AB = E_A + E_B). Configuration interaction (CI) based quantum chemistry methods are, in general, not size-consistent.⁷³ In particular, the only size-consistent CI methods are CI with singles (CIS) and FCI. On the other hand, single-reference CC methods are size-consistent as long as the form of wavefunction is parametrized by an exponential of the cluster operator. To reliably obtain the thermodynamic limit of large systems, size-consistency is necessary. At the thermodynamic limit, size-inconsistent methods approach just mean-field total energy and estimate no correlation energy. Therefore, CC methods, due to size-consistency, have stood out as a unique tool for simulating bulk systems.⁷⁴ 

We will show that ph-AFQMC is also size-consistent as long as the trial wavefunction is product separable. For isolated systems A and B, the supersystem Hamiltonian separates into Ĥ_A and Ĥ_B. Furthermore, these two operators commute since these systems are isolated. Therefore, the propagator is also product separable,

The HS transformation can be performed on exp(−ΔτĤ_A) and exp(−ΔτĤ_B) separately so that we have $B^AB=B^AB^B.$ This proves the size-consistency of free-projection AFQMC.

It can be also shown that the phaseless constraint is product separable. The overlap function in Eq. (7) can be written as

where the only assumptions we are making are (1) the product separability of the trial wavefunction: $|ΨTAB=|ΨTA|ΨTB$ and (2) the product separability of the slater determinant of nth walker: $|ψnAB=|ψnA|ψnB$⁠. The assumption (2) can be satisfied as long as we start from a product separable wavefunction since the propagator is product separable. With this overlap function, one can show that the importance function also obeys the product separability, and therefore, we conclude that ph-AFQMC is size-consistent.

The Hamiltonian for the uniform electron gas (UEG) is given simply as the sum of the kinetic energy and electron-electron interaction operator (up to a constant)

We will work with a basis of planewave spin orbitals $⟨rσ|Giσi⟩=1L3/2eiGi⋅rδσ,σi$⁠, where L is the length of the simulation cell, $Gi=2πLni$ for n_i a vector of integers, and σ_i is a spin index (either α or β). We impose a kinetic energy cutoff E_cut and work with a finite basis of 2M spin orbitals. In this basis, the kinetic energy is written as

and the electron-electron interaction operator is given by

where Ω = L³ is the simulation cell volume, Q is a momentum transfer vector that lives in an enlarged basis of size 4E_cut, and we have dropped the subscript index on G for simplicity. Finally, the Madelung energy E_M is included to account for the self-interaction of the Ewald sum under periodic boundary conditions.⁷⁵ For simplicity, we use the formula proposed by Schoof and co-workers⁷⁶ 

where N is the number of electrons in the simulation cell and $rs=3L34πN1/3$ is the dimensionless Wigner-Seitz radius.

The local energy ϵ_n(τ) for the UEG then reads

where the Coulomb two-body density matrix Γ_Q is

and the exchange two-body density matrix Λ_Q is

The formation of Γ_Q costs $O(M2)$⁠, whereas the formation of Λ_Q takes $O(M3)$ amount of work. Therefore, the evaluation of the exchange contribution is the bottleneck in the local energy evaluation. As noted in Ref. 39, the evaluation of the energy (and propagation) can be accelerated using fast Fourier transforms; however, we did not use this optimization here.

The two-body Hamiltonian $V^ee$ needs to be rewritten as a sum of squares to employ the AFQMC algorithm. It was shown in Ref. 39 that

where

and

with the momentum transfer operator $ρ^$ defined as

where Θ is the Heaviside step function. The Hubbard-Stratonovich operators $v^$ are now Â(Q) and $B^(Q)$⁠, and the rest of the AFQMC algorithm follows straightforwardly.

In ph-AFQMC, the main source of error is the bias introduced by the phaseless constraint. The magnitude of this bias is heavily dependent on the quality of trial wavefunctions. Although there are advanced options available for these such as multideterminantal trials and self-consistently determined single-determinantal trials,⁷⁷ we will employ a simple single determinant RHF trial wavefunction in most cases. In the UEG model, this is an M × N matrix (where N is the number of electrons and M is the number of planewaves) with 1’s on the diagonal entries.

Typically, for strongly correlated systems, it is useful to exploit essential symmetry breaking with HF wavefunctions. It is essential (as opposed to artificial) in the sense that the property of a single determinant wavefunction is qualitatively wrong without it. An attempt to exploit essential symmetry breaking typically leads to either spin-unrestricted HF (UHF) or spin-generalized HF (GHF) trial wavefunctions which have a lower energy than RHF. Indeed, such essential symmetry breaking was shown to be powerful when applying ph-AFQMC to bond dissociation of F₂.⁷⁸ 

An example of artificial symmetry breaking is buckminsterfullerene (C₆₀), a stable, electron paramagnetic resonance silent (EPR-silent) molecule.⁷⁹ There is a complex, GHF (cGHF) solution⁴⁸ for C₆₀ which was characterized to be an artifact due to the lack of treatment for weak correlation at the HF level.⁴⁹ In other words, orbital optimization in the presence of weak correlation such as the second-order Møller-Plesset theory would restore artificial symmetry breaking. Since both experiments^80,81 and computations⁴⁹ suggest that C₆₀ is a stable closed-shell molecule, the RHF state is more qualitatively correct than other broken-symmetry HF states. A detailed study of artificial vs essential symmetry breaking in ph-AFQMC will be published in our forthcoming paper.

The instability of RHF solutions is expected at all r_s values of the UEG model at the thermodynamic limit as proven by Overhauser.^82–84 As mentioned in Ref. 85, however, the R to U spin-symmetry breaking may not occur in the UEG model with a finite number of electrons. Instead, there is a critical Wigner-Seitz radius (⁠$rsc$⁠) below which no UHF solution exists. This is not surprising in the context of quantum chemistry since this is the same concept as “Coulson-Fischer points”⁸⁶ in molecules. Namely, when dissociating molecules, there exists a critical bond length where the R to U instability occurs. At bond distances closer than this, there is no genuine UHF solution.

In order to perform HF calculations, one must compute the effective one-body operator called the “Fock” operator defined as

where $D=CoccCocc†$ with C_occ being the occupied MO coefficient matrix. After some straightforward algebra, we find

where the kinetic energy matrix, T, reads

the Coulomb matrix, J, is

with

and the exchange matrix is given by

with

A similar derivation can be found in Ref. 85.

With the above Fock matrix, one can perform an HF calculation by optimizing the HF energy expression with respect to the orbital rotation parameter Θ_vo (a matrix of n_vir-by-n_occ) which relates two different MO coefficients via a unitary transformation,

where the antihermition matrix Δ, which is parametrized by $Θvo$⁠, is

The subscript of each matrix block denotes the dimension of the corresponding block, o = n_nocc and $v$ = n_nvir. An HF solution is defined as a stationary point that satisfies

where

The local stability of a given stationary point can then be tested by diagonalizing the orbital Hessian, M,

We will see whether there is essential symmetry breaking in the low r_s regime in the UEG model and try to utilize this essential symmetry breaking when appropriate.

Unless otherwise noted, the AFQMC calculations in this work were performed by a development version of QMCPACK.⁸⁷ PAUXY⁸⁸ was also used in the initial testing stages. Unless noted otherwise, AFQMC results below are obtained using QMCPACK. HANDE^89,90 was used to crosscheck our numbers for small systems that are not presented in this work. We used 0.005 a.u. for the time step Δτ throughout the paper. This was found to be enough for systems we considered here. A total of 2880 walkers were used, and the population bias from this was found to be negligible in the results reported here. The comb⁹¹ and pair branching⁹² population control algorithms are used in PAUXY and QMCPACK, respectively. The demonstration of the convergence of these parameters is available in the supplementary material.

All broken symmetry HF calculations were performed with a development version of Q-Chem,⁹³ and the details for the implementation can be found in Refs. 49, 94, and 95. The optimizer used for those HF calculations is geometry direct minimization (GDM) developed by Van Voorhis and Head-Gordon.⁹⁶ The internal stability analysis was performed for all HF solutions to ensure the local stability of each solution⁹⁷ where we used a finite-difference orbital Hessian using analytic orbital gradient.⁹⁸ All calculations were performed with periodic boundary conditions; no twist averaging was performed.

The UEG model has been explored by multiple methods at T = 0, and an extensive amount of benchmark data is already available. We compare our ph-AFQMC results against other methods and discuss whether the use of RHF trial wavefunction is reliable for r_s ≤ 5.0. It is expected that the quality of an RHF wavefunction degrades as r_s increases (approaching the atomic limit) since electrons tend to localize. r_s = 5.0 is a commonly investigated intermediate Wigner-Seitz radius, so it will be interesting to see how ph-AFQMC performs without employing more sophisticated trial wavefunctions.

For simplicity, we will refer to ph-AFQMC with an RHF trial wavefunction (RHF+ph-AFQMC) to as ph-AFQMC unless mentioned otherwise.

We summarize some interesting aspects in the HF solutions of the UEG model due to its simple form of Hamiltonian:

1. The MO coefficient matrix of an RHF solution is an identity matrix. This makes obtaining an UHF solution from solving an eigenvalue equation for F difficult in the following sense. C from diagonalizing F is always unitary, and the subsequent density matrix D is therefore identity. Since D is identical to the density matrix of RHF, one obtains an RHF solution immediately after one single diagonalization of F. A direct energy minimization⁹⁶ or the use of the HF projector⁸⁵ is necessary to obtain broken-symmetry HF states.

2. The RHF energy does not depend on the basis set size. However, broken-symmetry solutions depend on the basis set size. It is important to converge their energies with respect to M when discussing their existence.

3. Both Coulomb and kinetic energies are minimized in an RHF solution. In particular, the Coulomb energy is always zero in RHF.

4. The R to U symmetry breaking is driven by the lowering of exchange energy which compensates the increase in the Coulomb and kinetic energies.

We are interested in the paramagnetic phase of the UEG. As r_s increases, the ferromagnetic (i.e., spin-polarized) phase becomes the ground state.⁹⁹ The GHF solution can appear in the transition between these two phases at quite high r_s values. Other than this transition, a genuine GHF solution does not appear, and therefore, we study only the UHF solutions for the purpose of this study.

We discuss the UHF solutions in the 14-electron UEG model. In Fig. 1(a), we show the basis set convergence behavior of UHF. Unlike RHF, the UHF energy does depend on the basis set size. For the 14-electron UEG model, it is sufficient to converge the UHF energies over r_s ≤ 10.0 with M = 925. Based on Fig. 1(a), we see that the energy lowering from RHF to UHF starts to appear for r_s > 3.5. The critical Wigner-Seitz radius for the 14-electron model is $rsc∈(3.5,4.0]$⁠. This is more obvious from looking at the ⟨Ŝ²⟩ value of the UHF solutions for M = 943 as a function of r_s, as shown in Fig. 1(b). Nonzero ⟨Ŝ²⟩ indicates the appearance of a UHF solution. It is clear that the RHF solution becomes unstable above r_s = 3.5. The emergent strong correlation as increasing r_s is most obvious from looking at the momentum distribution (MD) [Fig. 1(c)] and natural orbital occupation numbers (NOONs) [Fig. 1(d)]. The MD is defined as the diagonal elements of a one-particle reduced density matrix,

The MD for the UEG model was thoroughly studied in Ref. 100. NOONs are the eigenvalues of a one-particle density matrix which is closely related to the MD. Both MD and NOONs show the increasing number of open-shell electrons as increasing the r_s values. The open-shell electrons appear for r_s > 3.5 which is consistent with Figs. 1(a) and 1(b).

We also studied larger UEG models, 54-electron and 114-electron, using UHF. As expected, there are many more local minima than the 14-electron model. This makes locating the global minimum even more challenging. These multiple minima lead to ambiguity in the subsequent ph-AFQMC calculations since there are many UHF solutions that can be used as trial wavefunctions. Thorough studies of the UHF solutions in the UEG model were presented in Ref. 85. For the purpose of this work, instead of trying to locate the global minimum, we investigated the critical Wigner-Seitz radius to determine whether one can employ UHF trial wavefunctions for simulating dense UEG models. We summarize the range for the critical radius for each UEG model studied in this work in Table I.

Given Table I, it is necessary to employ RHF trial wavefunctions for the 14-electron model at r_s ≤ 3.5, the 54-electron model at r_s ≤ 4.5, and the 114-electron model at r_s ≤ 2.5. We will study the 14-electron model at r_s = 0.5, 1.0, 2.0, 3.0, 5.0 and investigate the utility of various trial wavefunctions including RHF and UHF at r_s = 5.0. We will focus on using RHF trial wavefunctions for the 54- and 114-electron models at r_s = 0.5, 1.0, 2.0, where there is no R to U instability. For the 54-electron model at r_s = 5.0, there is a UHF solution with ⟨Ŝ²⟩ = 0.94 which exhibits marginal symmetry breaking.

We begin by studying the 14-electron UEG which was studied in detail by Shepherd and co-workers in Ref. 62. This small benchmark system is helpful as it is accessible to most quantum chemistry methods, whilst still exhibiting the typical challenges one faces when simulating real solids, namely, basis set incompleteness error, and strong correlation (when r_s is large). In addition, it has of late emerged as a standard benchmark system for the UEG.⁶⁶ 

The basis set convergence of wavefunction based quantum chemistry methods for the UEG has been explored a number of times by various methods,^61,101,102 and we will only briefly comment on it here. For our purposes, it is sufficient to note that the convergence to the CBS limit is slowest at high densities (low r_s), and thus, it is sufficient to converge the basis set error here. This slower convergence can be understood simply because the electron-electron cusp is more pronounced at high densities (the electron are more likely to coalesce). This is seen in Fig. 2 and Table II where on the order of 2000 PWs are necessary to converge the total energy to within 1 mHartree in absolute energy (not per electron). Similar to previous studies, we observe a more or less linear relationship between E_c and 1/M for M greater than 925.^61,66 The linear extrapolation to the CBS limit when using PWs was thoroughly studied and understood by Shepherd and co-workers in Ref. 101. We point out that the use of transcorrelated approaches in AFQMC could greatly accelerate the convergence to the CBS limit.¹⁰³ 

FCIQMC with the initiator approximation (i-FCIQMC) is a formally exact approach as long as there is no initiator bias. Therefore, comparing ph-AFQMC with i-FCIQMC is a good way to assess the accuracy of ph-AFQMC. We see that ph-AFQMC agrees well with i-FCIQMC within the error bar up to M = 179, and it starts to deviate from i-FCIQMC beyond that. However, i-FCIQMC numbers for the larger M values should be taken cautiously as it has been noted elsewhere that the bias from the initiator approximation was not completely removed⁶⁶ and that the i-FCIQMC results for r_s = 0.5 may be too low by approximately 1 mE_h.

Comparing ph-AFQMC and CC methods is perhaps more relevant for the purpose of this paper. The UEG model r_s = 0.5 is relatively weakly correlated, and thus, CC methods on top of an RHF reference work very well. Neufeld and Thom found that for r_s = 0.5, CCSDT is enough to converge the correlation energy with respect to the excitation levels.⁶⁶ Therefore, the CCSDT numbers in Table II should be considered to be exact for a given basis set. As it is clear from Table II, the CCSD correlation energies are all above those of ph-AFQMC and ph-AFQMC agrees with CCSDT up to sub-millihartree.

These results are particularly encouraging for the following reasons. This dense UEG model may be analogous to a weakly correlated molecular system, in the sense that for a finite number of electrons, it is relatively well described by HF theory. In such a system, CCSDT [or CCSD(T)] should be more or less exact. Even if their absolute energies were not exact, the relative energies such as barrier heights and interaction energies should be close to exact. The results here suggest that ph-AFQMC is a potentially powerful tool to handle such weakly correlated systems. For the rest of this section, we will assess the accuracy of ph-AFQMC for higher r_s where there can be a good mixture of weak and strong correlation. Furthermore, the quality of the RHF trial wavefunction will start to degrade, so we will show how this affects ph-AFQMC.

In Table III, we present the comparison of ph-AFQMC, i-FCIQMC, RCCSD, RCCSDT, and RCCSDTQ for selected basis sets. All of our ph-AFQMC is available in the supplementary material. At r_s = 1.0, ph-AFQMC agrees with i-FCIQMC within the error bar of each result when M = 1189. Small basis set (M = 179) results suggest that from RCCSDT to RCCSDTQ, only small correlation energy is gained. Therefore, we consider RCCSDT to be near-exact for larger basis sets. Near the CBS limit (M = 2109), we found that the ph-AFQMC energy is 15–16 mE_h lower than RCCSD and is within the error bar of RCCSDT.

At r_s = 2.0, ph-AFQMC agrees with i-FCIQMC within each error bar. However, RCCSDT struggles to obtain quantitatively accurate results for M = 925. RCCSD is about 36 mE_h above, and RCCSDT is about 5 mE_h above the i-FCIQMC (and ph-AFQMC) correlation energies. Like usual strongly correlated systems, the effect of quadruples is not negligible here, and it accounts for about 2 mE_h correlation energy in a small basis (M = 81). As shown in Table III, ph-AFQMC provides a lower correlation energy than even RCCSDTQ in the M = 81 basis set. Since neither ph-AFQMC nor RCCSDTQ is variational, such correlation energy comparisons should be taken with caution. Nevertheless, since ph-AFQMC and i-FCIQMC agree for M = 925, we expect ph-AFQMC to be accurate (i.e., near-exact) for M = 81. This result highlights the utility of RHF+ph-AFQMC. Namely, it can provide quantitatively accurate results when the role of quadruples is not negligible and yet still small enough for the RHF trial wavefunction to behave well.

Although at r_s = 0.5, 1.0, 2.0, ph-AFQMC provides more or less exact correlation energies as the density is lowered further, we find that the RHF trial wavefunctions performance degrades significantly. Not only does the ph-AFQMC correlation energy become above RCCSDTQ by about 2 mE_h at r_s = 3.0, but the stability of the simulations suffers noticeably. Nevertheless, at r_s = 3.0, we observe that ph-AFQMC is comparable to RCCSDT and is able to reach the CBS limit reliably.

However, r_s = 5.0 is much more difficult to handle with an RHF trial wavefunction. This is typically evidenced by an increase in the number of rare event population fluctuation. These rare events are well understood and arise due a divergent importance function which occurs when ⟨Ψ_T|ϕ⟩ approaches zero. Although these rare events can be effectively controlled by the use of bounds on the local (and/or hybrid energy),¹⁰⁶ they nevertheless signify a worsening in the quality of trial wavefunction for a fixed system size. To demonstrate this, in Fig. 3, we plot the convergence of the ph-AFQMC energy with projection time for a range of densities with M = 93 as well as an estimate for the overlap ∑_n$wn$|⟨Ψ_T|ϕ_n⟩|/∑_n$wn$⁠. We see from Fig. 3(a) that as r_s increases, the projection time necessary to converge to the ground state increases as well as the frequency of rare events. This is correlated with a decrease in the magnitude in the overlap as is seen from Fig. 3(b).

Indeed, at r_s = 5.0, we found that the ph-AFQMC energy did not converge monotonically with increasing basis set size past M = 389. Rather, the ph-AFQMC correlation energy decreases in magnitude with increasing basis set size. This signals a complete breakdown of the phaseless constraint with this trial wavefunction. We note that a similar effect can be observed in i-FCIQMC when the initiator error is not fully converged for increased basis set sizes, where one finds that the correlation energy begins to plateau as a function of basis set. This suggests that an improved trial wavefunction is necessary to attain sensible results for this system.

It is noteworthy to point out that this unusual behavior of ph-AFQMC energy with respect to the basis set size could indicate the “nonvariational” failure of ph-AFQMC. ph-AFQMC is formally nonvariational in the sense that a variational energy estimator of a given ph-AFQMC wavefunction can be above the mixed energy estimator in Eq. (15).¹⁴ Similarly, CC methods are also formally nonvariational due to their projective nature. With an RHF reference, it has shown catastrophic nonvariationality for strongly correlated systems.^6,107,108 It is possible that RCCSD (and even RCCSDT) is also exhibiting nonvariationality for this r_s value. This can be confirmed with more sophisticated CC methods.^109,110 The investigation of the nonvariationality of ph-AFQMC and CC methods in the context of strong correlation will be an interesting subject for future study.

As discussed in Sec. IV A, above r_s = 3.5, UHF solutions appear and they can be often powerful trial wavefunctions for ph-AFQMC.⁷⁸ As in Ref. 78, we employ the spin-projection technique to remove the spin-contamination completely. That is, we use the RHF wavefunction as the initial wavefunction while using the UHF wavefunction for the constraints. As shown in Sec. II A 3, because both RHF and UHF are size-consistent, the resulting UHF+ph-AFQMC performed with an RHF initial wavefunction must be also size-consistent. In Table IV, we present the correlation energies from UHF+ph-AFQMC. Given its substantial symmetry breaking at r_s = 5.0 shown in Fig. 1, it is surprising that it does not provide much improvement over RHF+ph-AFQMC. The UHF+ph-AFQMC energies are about 1 mH lower than the RHF+ph-AFQMC energies, and they are still far away from more accurate i-FCIQMC benchmark energies.

To investigate the ph-AFQMC results at r_s = 5.0 further, we explored the use of nonorthogonal multi-Slater determinant (NOMSD) expansions as trial wavefunctions generated using a version of the projected Hartree-Fock (PHF) algorithm developed by Scuseria and co-workers.^111–113 Interested readers are referred to Ref. 32 for further details. In Fig. 4, we find an initial rapid decrease in the error of the ph-AFQMC correlation energy and correspondingly a reduction in the ph-AFQMC statistical variance in the local energy estimator. The statistical variance in the estimator should not be confused with the wavefunction variance in variational MC (VMC). Therefore, it is possible to observe some improvement in the statistical error bar, while the energy estimator does not improve noticeably as observed before.³² This long tail in the convergence of the ph-AFQMC energy is indicative that the system is strongly correlated. We note that the FCI space for M = 57 contains on the order of 10¹⁶ determinants, and therefore, the improvement in the trial wavefunction via determinantal expansions is eventually limited by the exponential wall. We find similar behavior for larger basis sets.

Table IV summarizes our ph-AFQMC results using a 10 determinant expansion. Note that the M = 93 and M = 389 values are roughly within error bars of each other, despite the i-FCIQMC correlation energy decreasing by approximately 10 mE_h. We found that for even larger basis sets, the RHF+ph-AFQMC correlation energies lay above those at M = 93. Nevertheless, we see that by improving the trial wavefunction, the ph-AFQMC correlation energies begin to slowly approach those of i-FCIQMC values. The slow convergence is evidence of the limitation of using multi-Slater determinant trial wavefunctions in strongly correlated systems.

As explained in Sec. II A 3, ph-AFQMC is size-consistent and thus can reliably reach the thermodynamic limit using larger super cells along combined with finite size corrections and twist averaging.^{67,68,114–117} As the total energy of the UEG model in the thermodynamic limit is already well understood for solid state densities,¹¹⁸ here, we instead just study finite-sized UEG models and compare with other available methods when applicable.

Following the sequence of “magic numbers” in the UEG model, we study larger supercells (54 electrons and 114 electrons) with ph-AFQMC for r_s = 0.5, 1.0, 2.0. In the 14-electron UEG model, we obtained energies with error bars of the order of 1 mE_h. This is important for molecular applications where we aim for energy differences between two finite systems. On the other hand, the cost of achieving the same statistical error for larger systems adds an extra $O(N)$ to the computational cost of ph-AFQMC. This extra cost for sampling may be avoided by the correlated sampling technique,^119,120 but here, we instead compare the total energy per electron. This metric is well-suited for ab initio solids (or extended systems), in general.

High-quality DMC numbers are available for the 54-electron UEG model, and we compare ph-AFQMC against this. For the 114-electron UEG model, there are only variational MC (VMC) results available, so we will compare against these.

We found that from M = 1419 to M = 2109, the change in E_tot/N at r_s = 0.5 is only of the order of 0.1 mE_h. We therefore do not perform the CBS extrapolation for the comparison between ph-AFQMC and DMC. The reported that ph-AFQMC numbers are obtained from M = 1419 which enables a direct comparison between ph-AFQMC and i-FCIQMC at r_s = 0.5 and r_s = 1.0. The ph-AFQMC results for M = 2109 at r_s = 0.5, 1.0, 2.0 are available in the supplementary material.

In Table V, we summarize the comparison between ph-AFQMC, i-FCIQMC, and DMC for the 54-electron UEG model. At r_s = 0.5, ph-AFQMC and i-FCIQMC agree with each other within the error bar. Shepherd and co-workers found that the DMC-BF energy is somewhat higher than i-FCIQMC and suggested that the fixed-node error with the backflow (BF) trial wavefunction may not be small. Indeed, we reach the same conclusion with ph-AFQMC. As explained in Sec. I, the difference between DMC and AFQMC is mainly the discretization (or basis set) we work with. It is interesting that the fixed-node error in FN-DMC can be nonnegligible even with more sophisticated trial wavefunctions such as Slater-Jastrow (SJ) and BF. It is encouraging that we can achieve near-exact accuracy with ph-AFQMC at r_s = 0.5 with this simplest possible RHF trial wavefunction.

At r_s = 1.0, ph-AFQMC is in a better agreement (the difference is about 0.2 mE_h/e) with FN-DMC-BF than i-FCIQMC is. In this case, i-FCIQMC suffers from the initiator bias and results into about 1 mE_h/e above the FN-DMC-BF energy. In fact, the ph-AFQMC energy is lower than that of FN-DMC-BF, which may indicate nonnegligible fixed-node errors even in FN-DMC-BF. Since ph-AFQMC is not variational while FN-DMC-BF is, more careful calibration is highly desirable to see if there are indeed fixed-node errors in FN-DMC-BF. ph-AFQMC agrees better with FN-DMC-BF than does FN-DMC-SJ by 1 mE_h/e similarly to the r_s = 0.5 case.

No i-FCIQMC results are available at r_s = 2.0 due to the severity of the sign problem, so we instead must compare only to FN-DMC. We find that the ph-AFQMC energy is 1.1 mE_h/e below the FN-DMC-SJ energy and 0.3 mE_h/e below the FN-DMC-BF energy. With a larger basis set M = 2109, the ph-AFQMC energy lies 0.4 mE_h/e below the FN-DMC-BF energy, as shown in the supplementary material. Further increasing r_s to 5.0, we observe that ph-AFQMC is no longer comparable to FN-DMC-SJ and FN-DMC-BF as expected. As the use of UHF and NOMSD trials was found to be ineffective in the 14-electron model at r_s = 5.0, we did not perform such calculations here.

In summary, for the 54-electron UEG model at r_s = 0.5, 1.0, 2.0, we observe that ph-AFQMC can obtain nearly exact E_tot/N. In particular, its accuracy is comparable to other state-of-the-art methods such as i-FCIQMC and FN-DMC-BF. The general conclusions are similar to the 14-electron UEG model: ph-AFQMC is particularly well-suited for r_s values smaller than 5 where there exists moderate strong correlation. It is encouraging that ph-AFQMC achieved these highly accurate results using the simplest trial wavefunction, RHF.

Encouraged by the near-exact accuracy of ph-AFQMC for low r_s values in the 14- and 54-electron UEG models, we used ph-AFQMC to provide benchmark numbers for the 114-electron UEG model for future method development. The 114-electron UEG model has been relatively less explored. For determinant-based algorithms like i-FCIQMC, the sign problem is likely to preclude its application except for very high densities. On the other hand, for ph-AFQMC, this does not pose a significant challenge especially when considering r_s ≤ 2.0.

At r_s = 0.5, the total energy per electron changes by 0.5 mE_h/e when increasing M from 1419 to 2109. For higher r_s values, we expect this energy change to be smaller. We will present the ph-AFQMC energies at r_s = 0.5, 1.0, 2.0 all obtained with M = 2109. We expect that our ph-AFQMC energies reported here have the basis set incompleteness error of the order of 0.5 mE_h/e per electron. Therefore, the numbers reported here may be considered as an upper bound to the ph-AFQMC energies at the CBS limit.

The ph-AFQMC results are presented in Table VI. The only data available in literature are r_s = 1.0 with a VMC approach with a Slater-Jastrow wavefunction. The VMC energy is 0.60395(25) E_h/e¹²² which is at least 5 mE_h/e higher than our ph-AFQMC energies. Comparing ph-AFQMC energies in Tables V and VI, we note that the finite-size effect is still very large. Namely, the energy per electron is far from the convergence with respect to the system size. It will be interesting to investigate finite-size effects with ph-AFQMC in more realistic systems in the future. Although further comparisons are not possible due to the lack of benchmark data, we believe that the ph-AFQMC numbers in Table VI are close to the exact energies and the correlation energy error is smaller than 1 mE_h per electron given the results for the 54-electron model.

In this paper, we examined the performance of phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC) with the spin-restricted Hartree-Fock (RHF) trial wavefunction (i.e., RHF+ph-AFQMC) on the uniform electron gas (UEG) problem. We considered the 14-electron, 54-electron, and 114-electron UEG models. Through these studies, we found the following conclusions:

1. In the 14-electron case, we compared RHF+ph-AFQMC with spin-restricted coupled-cluster (RCC) methods. Compared to RCC with singles and doubles (RCCSD) and CC with singles, doubles, and triples (RCCSDT), RHF+ph-AFQMC performs better than RCCSDT and similarly to or slightly worse than RCCSDTQ for r_s ≤ 3.0.

2. For the 14-electron problem at r_s = 5.0 where CCSDT is inadequate, RHF+ph-AFQMC exhibits rare fluctuations in the energy estimator, which makes the phaseless approximation difficult to use. We found that using a small multideterminant trial wavefunction as well as a spin-unrestricted HF wavefunction is effective in stabilizing the simulations but still ineffective in obtaining highly accurate results in such cases.

3. In the case of the 54-electron UEG model, the comparison with initiator full configuration interaction QMC (i-FCIQMC) and fixed-node diffusion MC (FN-DMC) suggested that RHF+ph-AFQMC is a promising tool for simulating dense solids. Such connections between dense solids and the UEG model were previously made in Refs. 123–125. RHF+ph-AFQMC confirmed that the fixed-node error in FN-DMC for r_s = 0.5 as noted before in an i-FCIQMC study. Moreover, RHF+ph-AFQMC revealed that the initiator bias in i-FCIQMC for r_s = 1.0 is large (about 1 mE_h per electron) and the fixed-node error in FN-DMC with a back flow trial wavefunction (FN-DMC-BF) may not be negligible (0.3 mE_h per electron). A similar trend was observed in the case of r_s = 2.0. Finally, r_s = 5.0 was found to be challenging for RHF+ph-AFQMC to tackle, and the ph-AFQMC correlation energy was simply inadequate compared to FN-DMC-BF for this case. Overall, RHF+ph-AFQMC was found to be as accurate as or potentially more accurate than FN-DMC-BF wavefunction for r_s up to 2.0.

4. We produced RHF+ph-AFQMC energies of the 114-electron problem for r_s ≤ 2.0 where not many benchmark data are available. Given its performance for the 54-electron case, we expect the RHF+ph-AFQMC correlation energy error to be less than 1 mE_h per electron.

It is the central message of this paper that even with the simplest trial wavefunction (RHF), ph-AFQMC is a powerful tool for simulating molecules and solids where there is no noticeable strong correlation between electrons. In particular, its scope lies between CCSD and CCSDT. Given its low scaling (⁠$O(N3)−O(N4)$⁠), RHF+ph-AFQMC remains a promising tool.

The future study should include a more extensive benchmark of RHF+ph-AFQMC on more chemically relevant systems such as the W4-11¹²⁶ set as well as designing better and yet compact trial wavefunctions for AFQMC. Using dynamically correlated orbitals such as those from orbital-optimized Møller-Plesset perturbation theory can be an economical way to go beyond HF trial wavefunctions.^49,94,95 Some essential symmetry breaking in the HF trial wavefunction can potentially improve the performance of ph-AFQMC greatly such as using complex, restricted HF orbitals.^127,128 Finally, the finite-temperature extension of ph-AFQMC has been well established.^129–131 The assessment of ph-AFQMC for the warm dense UEG model,¹³² which has been the subject of intense research of late,^76,133–136 is currently work in progress.

The supplementary material includes HF and ph-AFQMC energies reported in this work, as well as timestep and population control bias tests.

The work of J.L. was partly supported by the CCMS summer internship in 2018 at the Lawrence Livermore National Lab. J.L. thanks Martin Head-Gordon and Soojin Lee for consistent encouragement. We would like to thank Carlos Jimenez-Hoyos for providing us access to the PHF code used to generate NOMSD expansions. The PHF code used in this work was developed in the Scuseria group at Rice University.^111–113 This work was performed under the auspices of the U.S. Department of Energy (DOE) by LLNL under Contract No. DE-AC52-07NA27344. Funding support was from the U.S. DOE, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and Center for Predictive Simulation of Functional Materials (CPSFM). Computer time was provided by the Livermore Computing Facilities.