The recently developed density matrix quantum Monte Carlo (DMQMC) algorithm stochastically samples the N-body thermal density matrix and hence provides access to exact properties of many-particle quantum systems at arbitrary temperatures. We demonstrate that moving to the interaction picture provides substantial benefits when applying DMQMC to interacting fermions. In this first study, we focus on a system of much recent interest: the uniform electron gas in the warm dense regime. The basis set incompleteness error at finite temperature is investigated and extrapolated via a simple Monte Carlo sampling procedure. Finally, we provide benchmark calculations for a four-electron system, comparing our results to previous work where possible.

The overwhelming majority of electronic structure studies of matter have been conducted at zero temperature. This state of affairs has been justified as typically one is interested in the low-energy properties of condensed matter systems or the room-temperature properties of chemical systems. Due to recent experimental advances, however, there has been renewed interest in the thermodynamic properties of electron plasmas such as those found on the pathways to inertial confinement fusion,1 in the interiors of Jupiter and other gas giants,2 at the surfaces of solids after laser irradiation,3 and in plasmonic catalysis.4 

Of fundamental importance to the theoretical understanding of these systems is the uniform electron gas (UEG), which has been pivotal to the development of modern quantum mechanical approaches to the low-temperature chemistry and physics of molecules and solids.5,6 At finite temperatures, the UEG can be described in terms of the density parameter, rs, and the degeneracy temperature, Θ = T/TF, where TF is the Fermi temperature. When both rs and Θ are of order one the system is said to be in the warm dense regime, with quantum, thermal, and interaction effects all being important. It is in this region where analytical methods, such as Ref. 7, are least useful and computational approaches such as quantum Monte Carlo are most important.

In a pioneering study, Brown et al.8 provided the first accurate data for the UEG in the warm dense regime using the restricted path integral Monte Carlo (RPIMC) method,9 from which the first accurate parameterizations of finite-temperature density functionals have been produced.10,11 Recent configuration path integral Monte Carlo (CPIMC) results12 have called into question the validity of the restricted path approximation used in Ref. 8, with significant disagreement between the two methods at high densities and low temperatures. Simulations using a third method such as density matrix quantum Monte Carlo (DMQMC) would help to resolve this discrepancy.13,41

The exact thermodynamic properties of the UEG can be determined from the (unnormalized) N-particle density matrix

ρ ˆ = e β H ˆ ,
(1)

where β = 1/kBT. A direct evaluation of Eq. (1) requires all eigenvalues and eigenstates of H ˆ to be known and is an impossible task in all but the simplest of systems. The infinite basis set limit can be approached by only including determinants that can be constructed using a finite basis set of M plane waves, reducing the problem to the diagonalization of an M N × M N matrix. Even this problem is only tractable for very small M and N.

Recently, Booth et al. have shown, through the development of the full configuration interaction QMC (FCIQMC) method, that FCI quality results can be obtained at zero temperature with no prior knowledge of the nodal structure of the wavefunction.14 FCIQMC also often offers a substantially reduced memory cost compared to conventional FCI calculations. This has most dramatically been seen in the case of the UEG, where a space of O ( 1 0 108 ) Slater determinants was successfully sampled using the initiator adaptation of FCIQMC.15,16 Subsequently, three of us used these ideas to develop the finite-temperature analogue of FCIQMC: DMQMC. This was applied to the Heisenberg model as a proof of concept17 but has not previously been used to study more realistic systems.

Here, we show how DMQMC can be applied to fermionic systems, starting with the UEG, thus opening the door to providing accurate, unbiased thermodynamic results for problems of chemical interest. We note that CPIMC18 and Krylov-projected FCIQMC19 will likely be complementary approaches to both DMQMC and PIMC in the treatment of real systems.

In Section II, we outline the DMQMC method and show how moving to the interaction picture provides substantial improvements in statistical accuracy when treating weakly correlated systems. In Section III, we discuss basis-set extrapolation at finite temperatures in detail. In Section IV, we present benchmark results for a four-electron system across the relevant parameter space, comparing, where possible, to previous results. Finally, in Section V, we outline the limitations and future prospects of DMQMC.

In this section, we briefly outline the DMQMC algorithm; a more complete description is available in Ref. 17. DMQMC is applicable to any Hamiltonian, but here we focus on the specific example of the UEG. We then explain how to sample the density matrix in the interaction picture, show that this overcomes sampling issues found when treating weakly correlated systems, and introduce a simple Monte Carlo scheme for sampling non-interacting density matrices in the canonical ensemble. Hartree atomic units are used throughout.

The unnormalized density matrix in Eq. (1) obeys the symmetrized Bloch equation

d ρ ˆ d β = 1 2 ( H ˆ ρ ˆ + ρ ˆ H ˆ ) ,
(2)

which can be solved using a simple Euler update scheme

ρ ˆ ( β + Δ β ) = ρ ˆ ( β ) Δ β 2 ( H ˆ ρ ˆ ( β ) + ρ ˆ ( β ) H ˆ ) + O ( Δ β 2 ) .
(3)

In DMQMC, Eq. (3) is solved stochastically by evolving a population of signed “psi-particles,” or “psips,” in a discrete operator space made of tensor products of Slater determinants. To this end, we rewrite Eq. (3) in matrix form,

ρ i j ( β + Δ β ) = ρ i j ( β ) Δ β 2 k ( H i k S δ i k ) ρ k j ρ i k ( H k j S δ k j )
(4)
= ρ i j ( β ) + Δ β 2 k ( T i k ρ k j + ρ i k T k j ) ,
(5)

where ρ i j = D i | ρ ˆ | D j , |Di〉 is a Slater determinant in the finite but large basis set, S is a variable shift introduced to control the psip population,14,20 and the last line defines the update matrix Tij = − (Hijij).

The rules for evolving the psips, which resemble those used in FCIQMC,14 follow from Eq. (5).

  1. Psips can spawn from a density matrix element ρik to ρij with probability ps(ikij) = Δβ|Tkj|/2, with sign(ρij) = sign(ρik) × sign(Tkj); a similar spawning process takes place from ρkj to ρij.

  2. Psips on the density matrix element ρij clone/die, whereby their population is increased or decreased, with probability p d ( i j ) = Δ β 2 | T i i + T j j | . The population is increased if sign(Tii + Tjj) × sign(ρij) > 0 and decreased otherwise.

Additionally, we annihilate psips of opposite sign on the same density matrix element to improve the efficiency of the algorithm and overcome the fermion sign problem.14,21 We note that, unlike PIMC methods where the quality of averages depends on the average sign of the sampled paths,8,18 in FCIQMC and DMQMC, we require a system specific and basis set dependent critical psip population to obtain correct low temperature and ground state estimates.14,21–23

The simplest starting point for a simulation is at β = 0, where the density matrix is the identity and can be sampled by occupying diagonal density matrix elements with uniform probability. A simulation then consists of propagating the initial distribution of psips with the rules described above to a desired value of β. Estimates for thermodynamic quantities can be found by averaging over many such simulations, a single one of which we call a “β-loop.”

In Fig. 1(a), we see that a direct application of DMQMC to the dense UEG can result in estimates for the internal or total energy that are too high in the intermediate temperature range. This can be understood by noting that, at rs = 1, the ground state of a few-electron UEG system is well described by a single (Hartree-Fock) determinant, |D0〉. The probability of initially selecting this determinant, however, is M N 1 , which rapidly approaches zero as M increases. If, by chance, the Hartree-Fock determinant or another low-energy determinant is sampled at β = 0, the population of psips arising from that low-energy determinant will dominate the simulation, but most simulations miss the low-energy part of the Hilbert space altogether. As shown in Fig. 1(b), this sampling problem reduces as the number of β-loops (or the population of psips per β-loop) increases, thus increasing the chance of sampling the low-energy space; however, this approach soon becomes impractical. Fig. 1(a) shows that, by moving to the interaction picture, we can effectively solve this sampling issue and regain FCI-quality thermodynamic averages.

FIG. 1.

Panel (a) shows the total energy calculated for a seven-electron spin-polarized electron gas at rs = 1 with M = 33 plane waves in the total momentum K = 0 subspace using DMQMC and an initial psip population Np = 104. Increasing the number of β-loops, Nβ, from 102 (squares) to 103 (crosses) results in a more accurate answer being reproduced. We also see that the error bars do not reflect the true errors for Nβ = 102 in the intermediate β regime. Panel (b) shows that the average occupation on the Hartree-Fock density matrix element (|D0〉〈D0|) is under-represented in the intermediate temperature range. Also plotted in (a) are the IP-DMQMC results (circles) for the total energy calculated using Np = 103 and only 10 β-loops (see Section II B).

FIG. 1.

Panel (a) shows the total energy calculated for a seven-electron spin-polarized electron gas at rs = 1 with M = 33 plane waves in the total momentum K = 0 subspace using DMQMC and an initial psip population Np = 104. Increasing the number of β-loops, Nβ, from 102 (squares) to 103 (crosses) results in a more accurate answer being reproduced. We also see that the error bars do not reflect the true errors for Nβ = 102 in the intermediate β regime. Panel (b) shows that the average occupation on the Hartree-Fock density matrix element (|D0〉〈D0|) is under-represented in the intermediate temperature range. Also plotted in (a) are the IP-DMQMC results (circles) for the total energy calculated using Np = 103 and only 10 β-loops (see Section II B).

Close modal

There are two sampling issues present when treating real systems; the distribution of weight in the density matrix changes rapidly as a function of β, and important determinants are rarely present in our initial configurations. Feynman originally pointed out that if we can write H ˆ = H ˆ 0 + V ˆ , where V ˆ is small compared to H ˆ 0 , then the quantity e β H ˆ 0 ρ ˆ will be a slowly varying function of β.24 This does not solve the issue of selecting important determinants, so we define an auxiliary matrix

f ˆ ( τ ) = e ( β τ ) H ˆ 0 e τ H ˆ ,
(6)

which has the properties

f ˆ ( τ = 0 ) = e β H ˆ 0 ,
(7)
f ˆ ( τ = β ) = e β H ˆ = ρ ˆ ( β ) .
(8)

From Eq. (7) we see that, by working with the operator f ˆ , we can start the simulation from e β H ˆ 0 instead of the identity. For most weakly correlated systems, this should provide a good first approximation to the distribution of weight in the fully interacting density matrix.

Differentiating Eq. (6) with respect to τ, we find

d f ˆ d τ = H ˆ 0 f ˆ f ˆ H ˆ
(9)
= V ˆ I ( τ β ) f ˆ ,
(10)

where we have used the usual definition of an operator in the interaction picture,

A ˆ I ( τ ) = e τ H ˆ 0 A ˆ e τ H ˆ 0 .
(11)

In practice, the exponential factors appearing in Eq. (10) are time-consuming to evaluate and we prefer to work with Eq. (9). Since we choose H ˆ 0 to be diagonal in our determinantal basis set of Hartree-Fock eigenstates, the only modification to the original DMQMC algorithm is that p d ( i j ) = Δ τ | H i i 0 H j j | . This scheme, which we dub interaction-picture DMQMC (IP-DMQMC), has the added benefit that there is typically no death along the diagonal as long as H i i 0 H i i . This overcomes a related issue in DMQMC simulations of large systems, whereby the weight of walkers on the diagonal decays nearly to zero with β; this was previously remedied with the use of importance sampling.17 

Eq. (8) shows that IP-DMQMC only samples the correct distribution at τ = β so that, unlike in DMQMC, where the whole temperature range is sampled in one simulation, separate simulations are required for each β value. As in DMQMC, estimates for observables require averaging over multiple independent simulations.

Referring back to Fig. 1 we see that working in the interaction picture effectively eliminates this sampling issue, with the correct total energy being reproduced using small numbers of psips and β-loops.

The choice of H ˆ 0 is somewhat arbitrary, but it should allow for an efficient sampling of f ˆ ( τ = 0 ) and this is most easily achieved if H ˆ 0 is non-interacting. In principle, any initial density matrix can be sampled using the Metropolis algorithm,25 but we have found this approach problematic due to the long equilibration times required at low temperatures and in large basis sets. An alternative method, which is free from such issues, is to use knowledge of the grand canonical density matrix corresponding to H ˆ 0 and sample this in such a way that the desired, canonical, distribution is reached.

Consider the grand canonical density matrix

ρ ˆ GC 0 = e β ( H ˆ 0 μ N ˆ ) ,
(12)

where H ˆ 0 = i ε i c ˆ i c ˆ i is some non-interacting Hamiltonian whose single-particle eigenvalues εi are known, N ˆ is the number operator, and μ is the chemical potential, which can be determined from the implicit relationship

N ˆ GC = i 1 e β ( ε i μ ) + 1
(13)
= i p i ,
(14)

where we have identified the usual Fermi-factor, pi. The grand canonical partition function, ZGC, can be written as

Z GC = N { n i N } e β i ( ε i μ ) n i ,
(15)

where { n i N } denotes a set of occupation numbers such that i n i = N and ni ∈ {0, 1} for fermions.

The probability of selecting a particular set { n ̄ i N } is

P GC ( { n ̄ i N } ) = 1 Z GC n i { n ̄ i N } e β ( ε i μ ) n i .
(16)

However, we wish to generate determinants in the canonical ensemble where the correct probability is

P C ( { n ̄ i N } ) = 1 Z C n i { n ̄ i N } e β ( ε i μ ) n i ,
(17)

where

Z C = { n i N } e β i ( ε i μ ) n i .
(18)

We see that P C ( { n i N } ) P GC ( { n i N } ) . Thus, by independently occupying orbitals with probability pi and then discarding those configurations with N ˆ N , we attain the correct proportionality factor ZGC/ZC. Only about one in N of the configurations sampled has the right value of N, but the sampling process is so fast that very little computer time is required for any system size in the reach of many-body simulation methods. The chemical potential can be obtained by numerically inverting Eq. (13) in the appropriate finite basis set, a procedure we carry out using Scipy.26 A demonstration of the above procedure is given in Fig. 2, where we see that H ˆ is indeed a slowly varying function of τ and that the correct estimate is reproduced at τ = β.

FIG. 2.

Variation of H ˆ with τ using H ˆ 0 = k ε k c ˆ k c ˆ k , i.e., the free-electron expression with ε k = 1 2 k 2 . The grand canonical procedure described in Section II C was used. The system shown is a four-electron, spin polarized gas at rs = 1, M = 81, and β = 1. For these results, we used approximately 103 psips and averaged over 100 simulations. The dashed line represents the exact FCI result, which IP-DMQMC reproduces at τ = β as expected. For comparison, H ˆ = 50 . 751 ( 4 ) hartree at β = 0.

FIG. 2.

Variation of H ˆ with τ using H ˆ 0 = k ε k c ˆ k c ˆ k , i.e., the free-electron expression with ε k = 1 2 k 2 . The grand canonical procedure described in Section II C was used. The system shown is a four-electron, spin polarized gas at rs = 1, M = 81, and β = 1. For these results, we used approximately 103 psips and averaged over 100 simulations. The dashed line represents the exact FCI result, which IP-DMQMC reproduces at τ = β as expected. For comparison, H ˆ = 50 . 751 ( 4 ) hartree at β = 0.

Close modal

Finally, we note that any diagonal density matrix can be obtained by reweighting the configurations which result from the above sampling procedure as

P new ( { n ̄ i N } ) = P old ( { n ̄ i N } ) e β ( E new E old ) ,
(19)

where Enew and Eold are the new and old total energies of a given configuration { n ̄ i N } , respectively.

To treat the UEG using DMQMC, we need to work in a finite basis set of M plane waves; thermodynamic quantities will therefore need to be extrapolated to the M → ∞ limit. At T = 0, it has been found16,27,28 that the correlation energy, Ec = EEHF, for unpolarized systems converges like M−1. Recent CPIMC results29 were also obtained using an M−1 extrapolation, although in principle this relationship only holds for an unpolarized system; we have found that, for polarized systems, extrapolating with M−5/3 results in a better fit,30 consistent with similar observations by other authors.31,32 Based on the analysis presented here, Schoof et al.12 have used the M−5/3 extrapolation with CPIMC. In addition, we find that the convergence of the total energy is strongly dependent on temperature.

At T > 0, there is a competition between the convergence of the kinetic and potential energies with M. To investigate this further we focus on a two-electron spin-polarized system, which can be solved exactly using diagonalization in large basis sets. In Fig. 3, we see this competition between energy scales: the total energy initially increases rapidly with basis-set size before appearing to saturate. As the size of the basis set is further increased, a slight reduction in the total energy is observed, with the residual error apparently proportional to M−5/3 (see inset of Fig. 3). In this regime, the convergence of the total energy is dominated by exchange and correlation effects.

FIG. 3.

Behavior of the FCI total energy with basis-set size for an N = 2, rs = 1 spin-polarized system at Θ = 0.5. Here, we see the competition between the exponential convergence of the total energy at low M (main plot) and the M−5/3 behavior for high M (inset).

FIG. 3.

Behavior of the FCI total energy with basis-set size for an N = 2, rs = 1 spin-polarized system at Θ = 0.5. Here, we see the competition between the exponential convergence of the total energy at low M (main plot) and the M−5/3 behavior for high M (inset).

Close modal

The initial increase of the total energy with respect to M at non-zero temperatures can be understood by looking at the non-interacting total energy as a function of basis-set size, which is most easily analyzed in the grand canonical ensemble. The non-interacting basis-set error is

Δ E 0 ( M ) = E 0 ( ) E 0 ( M )
(20)
= k > k c ε k p k ,
(21)

where kc is the plane-wave cutoff. For ε c = 1 2 k c 2 1 , this can be approximated as

Δ E 0 ε c ε 3 / 2 e β ( ε μ ) d ε ,
(22)

where we have used pkeβkμ) for εkμ. Hence,

Δ E 0 ( M ) 0 ( ε c + x ) 3 / 2 e β ( ε c + x μ ) d x
(23)
ε c 3 / 2 e β ( ε c μ ) 0 e β x d x ,
(24)

if εcβ−1 so that ( ε c + x ) 3 / 2 ε c 3 / 2 everywhere eβx is significant. It then follows that the leading-order correction is

Δ E 0 ( M ) β 1 ε c 3 / 2 e β ( ε c μ ) .
(25)

From Eq. (25), we see that the kinetic energy begins to converge exponentially once εcβ−1 or M V ( Θ r s 2 ) 3 / 2 N Θ 3 / 2 , where V is the simulation-cell volume. In practice, we find that for large Θ the kinetic energy and hence the total energy converge quite slowly. This is an issue for DMQMC simulations as the cost of a calculation increases with basis-set size.

We can mitigate some of these issues by instead extrapolating the temperature-dependent correlation energy, Ec(β, M) = E(β, M) − EHF(β, M). The infinite basis-set total energy can then be reconstructed as E(β, M = ∞) = EHF(β, M = ∞) + Ec(β, M = ∞). We calculate the “Hartree-Fock” energy, E HF = H ˆ HF , using the density matrix33,42

ρ ˆ HF = i e β E i HF | D i D i | ,
(26)

where E i HF = D i | H ˆ | D i and the sum runs over all determinants in the basis set. EHF(β) can be found using the sampling procedure outlined in Section II C, i.e.,

E HF = 1 Z HF i E i HF e β E i HF
(27)
= 1 Z HF i E i HF e β ( E i HF E i 0 ) e β E i 0
(28)
= i E i HF w ( i ) p ( i ) i w ( i ) p ( i ) ,
(29)

where w ( i ) = e β ( E i HF E i 0 ) and p ( i ) = Z 0 1 e β E i 0 . Thus, by generating determinants as described in Section II C and reweighting them using w(i), we can instead sample ρ ˆ HF and, as a result, estimate EHF as desired. In Fig. 4 we show the convergence of EHF(β, M) as a function of basis set for a four-electron, spin-polarized system at rs = 1 and Θ = 4. Note the large basis-set sizes required to converge the total energy to within statistical error bars. Fig. 4 also shows various other “non-interacting” or “mean-field” energy estimates as functions of M. Any of these could in principle be subtracted from E(β, M) to define a correlation energy, but the quantity defined by subtracting EHF(β, M) extrapolates most smoothly to the infinite M limit. The non-interacting grand canonical energy, H ˆ 0 GC , is significantly larger than the canonical estimates.

FIG. 4.

Exponential convergence of the Hartree-Fock total energy with basis-set size for N = 4, rs = 1, and Θ = 4. Note that no Madelung contribution is included for the Hartree-Fock estimates in this figure. The infinite basis-set limit for EHF is estimated as 70.792(1) hartree. For comparison, we plot H ˆ 0 and H ˆ 0 0 , where the trace is now with respect to the non-interacting density matrix. Also plotted is the total energy calculated in the grand canonical ensemble using a finite basis set, i.e., H ˆ 0 GC = k k c ε k p k .

FIG. 4.

Exponential convergence of the Hartree-Fock total energy with basis-set size for N = 4, rs = 1, and Θ = 4. Note that no Madelung contribution is included for the Hartree-Fock estimates in this figure. The infinite basis-set limit for EHF is estimated as 70.792(1) hartree. For comparison, we plot H ˆ 0 and H ˆ 0 0 , where the trace is now with respect to the non-interacting density matrix. Also plotted is the total energy calculated in the grand canonical ensemble using a finite basis set, i.e., H ˆ 0 GC = k k c ε k p k .

Close modal

Fig. 5 shows how Ec(β, M) depends on M at a number of different temperatures. For small basis sets, Ec shows a power-law decay with M, but this ceases for large enough M and the energies begin to increase again. We believe that the increase is due to kinetic effects that are not captured in the non-interacting expression we subtract. The value of M at which the correlation energy begins to increase again corresponds to the onset of the power-law convergence of the total energy observed in Fig. 3. This non-variational behavior of the internal energy with respect to M is not surprising as, at finite temperatures, it is the free energy that satisfies a variational principle.

FIG. 5.

Comparison of the convergence of Ec with basis-set size calculated using exact diagonalization for a two-electron spin-polarized system at rs = 1 for different values of Θ. The Θ = 2 data have been shifted down by 0.0002 hartree for visibility. The dependence of Ec on Θ is non-monotonic. The converged-basis-set correlation energies at Θ = 8 and Θ = 2 are very similar, but the correlation energy at Θ = 4 is more negative.

FIG. 5.

Comparison of the convergence of Ec with basis-set size calculated using exact diagonalization for a two-electron spin-polarized system at rs = 1 for different values of Θ. The Θ = 2 data have been shifted down by 0.0002 hartree for visibility. The dependence of Ec on Θ is non-monotonic. The converged-basis-set correlation energies at Θ = 8 and Θ = 2 are very similar, but the correlation energy at Θ = 4 is more negative.

Close modal

We can estimate Ec(β, M = ∞) by taking the value calculated with the largest basis set after the minimum is reached. This will in general over estimate Ec (it will be too negative) but the remaining discrepancy is typically smaller than the stochastic error bar. The systematic errors left after extrapolating Ec(β, M) in this manner are well within chemical accuracy (∼kBT) and can be orders of magnitude smaller than those introduced by a direct extrapolation of the total energy.

All calculations discussed in this paper were performed using the HANDE code.34 Unless otherwise stated, the QMC calculations used real amplitudes35–37 to sample the density matrix, which improves stochastic efficiency compared to integer weights. The full data set is available in the supplementary material.38 

Using the procedures outlined above, we are now in a position to provide exact benchmarks for the UEG in small simulation cells across the relevant parameter space. In this first study we focus on the four-electron spin-polarized system, which is the smallest non-trivial system and one for which there already exist benchmark calculations.29 All energies contain the Madelung contribution39 where appropriate. As a first step, we compare our four-electron IP-DMQMC results to FCI results in small basis sets and see perfect agreement across the whole temperature range (Fig. 6).

FIG. 6.

Comparison of IP-DMQMC results (markers) with FCI results (dashed-lines) for the UEG with N = 4 and rs = 1 in two different basis-sets. The inset shows the low T behavior where we see increasing the basis-set size serves to decrease the total energy in contrast to the high T behavior, where the opposite occurs. These calculations used integer rather than real weights and were run using approximately 1000 psips and we averaged over 100 simulations.

FIG. 6.

Comparison of IP-DMQMC results (markers) with FCI results (dashed-lines) for the UEG with N = 4 and rs = 1 in two different basis-sets. The inset shows the low T behavior where we see increasing the basis-set size serves to decrease the total energy in contrast to the high T behavior, where the opposite occurs. These calculations used integer rather than real weights and were run using approximately 1000 psips and we averaged over 100 simulations.

Close modal

We have extended these results to basis sets far beyond the reach of conventional full diagonalization procedures; the largest space sampled here contains approximately 1022 density matrix elements. We used the initialization procedure outlined in Section II C and the free-electron Hamiltonian for H ˆ 0 for rs ≤ 1; for rs > 1 we found it advantageous to use Hartree-Fock density matrix defined in Eq. (26). The calculations were initialized with 103–107 psips and the results averaged over 100–5000 simulations, each using a different random number seed. Time steps Δτ ranging from 0.01/EF to 0.001/EF were used, with a smaller time step required at lower rs; the values chosen were small enough that we could resolve no time step error within the statistical errors. Each (rs, Θ, M) calculation was typically run for 2 h on 48 cores with a total computational cost of approximately 80 000 core hours. The separate calculations of EHF required 9000 core hours.

In Figs. 7 and 8, we show the convergence of the IP-DMQMC results with basis set at low and high temperatures, respectively. We note that the behavior found in the two-electron system is also found in the four-electron system; in particular the non-trivial dependence of the correlation energy on M is reproduced. We find that a direct extrapolation of the total energy with respect to M is best for Θ ≤ 0.25 as it is here where kinetic effects are minimal and there is a clear trend in the total energy for the basis set sizes considered. The procedure outlined in Section III is best suited for temperatures above this, becoming increasingly useful for Θ ≥ 2 as more highly excited states become accessible, requiring prohibitively large basis sets for a direct extrapolation to be possible. In between these two regimes, both methods produce statistically identical results.38 Fig. 9 summarizes our results and shows perfect agreement with the available CPIMC data from Ref. 29. Further results at higher temperatures and other rs values are available in tabular form in the supplementary material and again agree with the available CPIMC results.38 

FIG. 7.

Total energy of a four-electron spin-polarized system at rs = 1 and Θ = 0.0625, showing a convergence with M−5/3. The dashed line represents an extrapolation to the infinite-basis-set limit carried out using a least-squares fit as implemented in Scipy.26 At this low temperature, a direct extrapolation of the total energy works better than the correlation-energy extrapolation technique discussed in Section III.

FIG. 7.

Total energy of a four-electron spin-polarized system at rs = 1 and Θ = 0.0625, showing a convergence with M−5/3. The dashed line represents an extrapolation to the infinite-basis-set limit carried out using a least-squares fit as implemented in Scipy.26 At this low temperature, a direct extrapolation of the total energy works better than the correlation-energy extrapolation technique discussed in Section III.

Close modal
FIG. 8.

Convergence of Ec with basis-set size for the four-electron system at rs = 5 calculated using IP-DMQMC. Note the similar behaviour to that found in the two-electron case in Fig. 5. At the high temperatures considered here, the correlation-energy extrapolation technique introduced in Section III works much better than the total-energy extrapolation illustrated in Fig. 7.

FIG. 8.

Convergence of Ec with basis-set size for the four-electron system at rs = 5 calculated using IP-DMQMC. Note the similar behaviour to that found in the two-electron case in Fig. 5. At the high temperatures considered here, the correlation-energy extrapolation technique introduced in Section III works much better than the total-energy extrapolation illustrated in Fig. 7.

Close modal
FIG. 9.

(a) Extrapolated total energies per particle for the four-electron system at rs = 1, 2, 4 showing exact agreement with the CPIMC results of Ref. 29. Dashed lines are meant as guides to the eye. (b) Relative deviation, (EDMQMCECPIMC)/EDMQMC, as a function of temperature showing statistically identical results for rs = 1 (circles), rs = 2 (diamonds), and rs = 4 (squares). Further results at higher temperatures and other rs values are available in tabular form in the supplementary material38 and again agree with the available CPIMC results.

FIG. 9.

(a) Extrapolated total energies per particle for the four-electron system at rs = 1, 2, 4 showing exact agreement with the CPIMC results of Ref. 29. Dashed lines are meant as guides to the eye. (b) Relative deviation, (EDMQMCECPIMC)/EDMQMC, as a function of temperature showing statistically identical results for rs = 1 (circles), rs = 2 (diamonds), and rs = 4 (squares). Further results at higher temperatures and other rs values are available in tabular form in the supplementary material38 and again agree with the available CPIMC results.

Close modal

In this paper, we have demonstrated how DMQMC can be applied to realistic systems. By moving to the interaction picture, we have removed sampling issues found when treating weakly correlated systems with large basis sets.

We have examined in detail the convergence of the total and correlation energies with respect to basis-set size M and temperature using a system accessible to FCI calculations. We found that, in general, these quantities exhibit a non-trivial dependence on M attributable to the competing energy scales present. By developing a simple Monte Carlo sampling scheme, we showed that it is possible to reduce the error in extrapolating these quantities to the complete basis-set limit by at least an order of magnitude at high temperatures. We believe this analysis and developments will be useful when treating molecular systems at finite temperatures. In addition, our approach to calculating the “Hartree-Fock” energy should be a useful first approximation when providing accurate benchmark calculations for systems away from the thermodynamic limit and in analyzing single-particle finite-size effects for the UEG at non-zero temperatures.

Using these developments, we have reproduced the four-electron CPIMC benchmarks of Ref. 29 and provided results at higher temperatures. We hope that these small-system results will aid the analysis of the apparent discrepancies between other QMC methods for larger system sizes40 and serve as benchmarks for other QMC methods based in configuration space.

Whilst the results presented here are for much smaller systems than those accessible by RPIMC and CPIMC, DMQMC provides access to exact finite-temperature data for a given basis set. The main limitation on the system size is the critical population (determined by the plateau height17,21) required to sample the density matrix. There are several grounds for optimism. The sign problem is much weaker at higher temperatures, implying that larger systems will be accessible albeit over a restricted temperature range. Further, in our previous study17 we found that the plateau height in DMQMC, which is a measure of the strength of the sign problem, was roughly the square of that in FCIQMC, which would suggest a rather limited utility of our method. Remarkably, however, we have found that, for the UEG at various rs values and for some other models, the plateau height is only a small multiple of the FCIQMC plateau height. For example, for N = 4, rs = 5, and M = 1045, the FCIQMC plateau height is roughly 8000 psips in comparison to 90 000 in IP-DMQMC at Θ = 0.0625. Finally, there is evidence that methodological developments will increase the scope of systems that can be treated with DMQMC. Given that the initiator approximation15 enabled FCIQMC to be applied successfully to large UEG systems across a range of densities,16,27 we are confident that, using similar approximations and developments in importance sampling, DMQMC will become a competitive method in treating degenerate Fermi systems.

F.D.M. is funded by an Imperial College Ph.D. Scholarship. N.S.B. acknowledges Trinity College, Cambridge for funding. J.J.S. thanks the Royal Commission for the Exhibition of 1851 for a Research Fellowship. J.S.S. and W.M.C.F. acknowledge the research environment provided by the Thomas Young Centre under Grant No. TYC-101. W.M.C.F. received support from the UK Engineering and Physical Sciences Research Council under Grant No. EP/K038141/1. Computing facilities were provided by the High Performance Computing Service of Imperial College London, by the Swiss National Supercomputing Centre (CSCS) under Project ID No. s523, and by ARCHER, the UK National Supercomputing Service, under EPSRC Grant No. EP/K038141/1 and via a RAP award.

1.
M.
Koenig
,
A.
Benuzzi-Mounaix
,
A.
Ravasio
,
T.
Vinci
,
N.
Ozaki
,
S.
Lepape
,
D.
Batani
,
G.
Huser
,
T.
Hall
,
D.
Hicks
,
A.
MacKinnon
,
P.
Patel
,
H. S.
Park
,
T.
Boehly
,
M.
Borghesi
,
S.
Kar
, and
L.
Romagnani
,
Plasma Phys. Controlled Fusion
47
,
B441
(
2005
).
2.
J. J.
Fortney
,
S. H.
Glenzer
,
M.
Koenig
,
B.
Militzer
,
D.
Saumon
, and
D.
Valencia
,
Phys. Plasmas
16
,
041003
(
2009
).
3.
R.
Ernstorfer
,
M.
Harb
,
C. T.
Hebeisen
,
G.
Sciaini
,
T.
Dartigalongue
, and
R. J. D.
Miller
,
Science
323
,
1033
(
2009
).
4.
S.
Mukherjee
,
F.
Libisch
,
N.
Large
,
O.
Neumann
,
L. V.
Brown
,
J.
Cheng
,
J. B.
Lassiter
,
E. A.
Carter
,
P.
Nordlander
, and
N. J.
Halas
,
Nano Lett.
13
,
240
(
2013
).
5.
D. M.
Ceperley
and
B. J.
Alder
,
Phys. Rev. Lett.
45
,
566
(
1980
).
6.
J. P.
Perdew
and
A.
Zunger
,
Phys. Rev. B
23
,
5048
(
1981
).
7.
F.
Perrot
and
M. W. C.
Dharma-wardana
,
Phys. Rev. A
30
,
2619
(
1984
).
8.
E. W.
Brown
,
B. K.
Clark
,
J. L.
DuBois
, and
D. M.
Ceperley
,
Phys. Rev. Lett.
110
,
146405
(
2013
).
9.
D. M.
Ceperley
,
Rev. Mod. Phys.
67
,
279
(
1995
).
10.
E. W.
Brown
,
J. L.
DuBois
,
M.
Holzmann
, and
D. M.
Ceperley
,
Phys. Rev. B
88
,
081102
(
2013
).
11.
V. V.
Karasiev
,
T.
Sjostrom
,
J.
Dufty
, and
S. B.
Trickey
,
Phys. Rev. Lett.
112
,
076403
(
2014
).
12.
T.
Schoof
,
S.
Groth
,
J.
Vorberger
, and
M.
Bonitz
, “
Ab initio thermodynamic results for the degenerate electron gas at finite temperature
,” preprint arXiv:1502.04616 (
2015
).
13.

There has also been disagreement reported at high densities between RPIMC and direct path integral Monte Carlo,41 although the origin of these discrepancies remains unclear.

14.
G. H.
Booth
,
A. J. W.
Thom
, and
A.
Alavi
,
J. Chem. Phys.
131
,
054106
(
2009
).
15.
D.
Cleland
,
G. H.
Booth
, and
A.
Alavi
,
J. Chem. Phys.
132
,
041103
(
2010
).
16.
J. J.
Shepherd
,
G.
Booth
,
A.
Grüneis
, and
A.
Alavi
,
Phys. Rev. B
85
,
081103
(
2012
).
17.
N. S.
Blunt
,
T. W.
Rogers
,
J. S.
Spencer
, and
W. M. C.
Foulkes
,
Phys. Rev. B
89
,
245124
(
2014
).
18.
T.
Schoof
,
M.
Bonitz
,
A.
Filinov
,
D.
Hochstuhl
, and
J. W.
Dufty
,
Contrib. Plasma Phys.
51
,
687
(
2011
).
19.
N. S.
Blunt
,
A.
Alavi
, and
G. H.
Booth
, “
Krylov-projected quantum Monte Carlo
,” preprint arXiv:1409.2420 (
2014
).
20.
C. J.
Umrigar
,
M. P.
Nightingale
, and
K. J.
Runge
,
J. Chem. Phys.
99
,
2865
(
1993
).
21.
J. S.
Spencer
,
N. S.
Blunt
, and
W. M. C.
Foulkes
,
J. Chem. Phys.
136
,
054110
(
2012
).
22.
M. H.
Kolodrubetz
,
J. S.
Spencer
,
B. K.
Clark
, and
W. M. C.
Foulkes
,
J. Chem. Phys.
138
,
024110
(
2013
).
23.
J. J.
Shepherd
,
G. E.
Scuseria
, and
J. S.
Spencer
,
Phys. Rev. B
90
,
155130
(
2014
).
24.
R. P.
Feynman
,
Statistical Mechanics: A Set of Lectures
(
Westview Press
,
1998
).
25.
N.
Metropolis
,
A. W.
Rosenbluth
,
M. N.
Rosenbluth
,
A. H.
Teller
, and
E.
Teller
,
J. Chem. Phys.
21
,
1087
(
1953
).
26.
E.
Jones
,
T.
Oliphant
,
P.
Peterson
 et al, SciPy, a Python-based ecosystem of open-source software for mathematics, science, and engineering, 2001, see http://www.scipy.org/citing.html.
27.
J. J.
Shepherd
,
G. H.
Booth
, and
A.
Alavi
,
J. Chem. Phys.
136
,
244101
(
2012
).
28.
J. J.
Shepherd
,
A.
Grüneis
,
G. H.
Booth
,
G.
Kresse
, and
A.
Alavi
,
Phys. Rev. B
86
,
035111
(
2012
).
29.
T.
Schoof
,
S.
Groth
, and
M.
Bonitz
,
Contrib. Plasma Phys.
55
,
136
(
2015
).
30.

This power law can be derived by analyzing the MP2 correlation energy of a polarized system in a manner analogous to the treatment of the unpolarized electron gas in Ref. 28.

31.
J.
Klimeš
,
M.
Kaltak
, and
G.
Kresse
,
Phys. Rev. B
90
,
075125
(
2014
).
32.
A.
Gulans
,
J. Chem. Phys.
141
,
164127
(
2014
).
33.

This differs from the usual approach to finite-temperature Hartree-Fock theory,42 in which the density matrix is obtained by minimizing the free energy over a choice of all non-interacting trial density matrices, with the single-particle eigenvalues and (in general) orbitals allowed to depend on temperature. Here, in place of a sum of single-particle eigenvalues, we use the many-particle total energies of specific determinants in our Boltzmann factors. Fortunately, since our “Hartree-Fock” density matrix is used only to improve the efficiency of the IP-DMQMC algorithm, the simulations yield exact results whether or not the standard definition is used.

34.
See http://www.hande.org.uk/ for the directions on how to acquire the source code.
35.
F. R.
Petruzielo
,
A. A.
Holmes
,
H. J.
Changlani
,
M. P.
Nightingale
, and
C. J.
Umrigar
,
Phys. Rev. Lett.
109
,
230201
(
2012
).
36.
C.
Overy
,
G. H.
Booth
,
N. S.
Blunt
,
J. J.
Shepherd
,
D.
Cleland
, and
A.
Alavi
,
J. Chem. Phys.
141
,
244117
(
2014
).
37.
N. S.
Blunt
,
S. D.
Smart
,
J. A. F.
Kersten
,
J. S.
Spencer
,
G. H.
Booth
, and
A.
Alavi
,
J. Chem. Phys.
142
,
184107
(
2015
).
38.
See supplementary material at http://dx.doi.org/10.1063/1.4927434 for the full QMC data set as well as results at further temperatures and rs values.
39.
L. M.
Fraser
,
W. M. C.
Foulkes
,
G.
Rajagopal
,
R. J.
Needs
,
S. D.
Kenny
, and
A. J.
Williamson
,
Phys. Rev. B
53
,
1814
(
1996
).
40.

A full comparison will be presented elsewhere.

41.
V. S.
Filinov
,
V. E.
Fortov
,
M.
Bonitz
, and
Z.
Moldabekov
,
Phys. Rev. E
91
,
033108
(
2015
).

Supplementary Material