Ab Initio Quantum Monte Carlo Simulations of the Uniform Electron Gas without Fixed Nodes II: Unpolarized Case

0.78 MB·

by  T. Dornheim

Ab Initio Quantum Monte Carlo Simulations of the Uniform Electron Gas
without Fixed Nodes II: Unpolarized Case
T. Dornheim1,†, S. Groth1,†, T. Schoof1, C. Hann1,2, and M. Bonitz1∗
†These authors contributed equally to this work.
1Institut fu¨r Theoretische Physik und Astrophysik,
Christian-Albrechts-Universit¨at zu Kiel, D-24098 Kiel, Germany
2Department of Physics, Duke University, Durham, North Carolina 27708, USA
(Dated: January 20, 2016)
In a recent publication [S. Groth et al., PRB (2016)], we have shown that the combination of two
novel complementary quantum Monte Carlo approaches, namely configuration path integral Monte
Carlo (CPIMC) [T. Schoof et al., PRL 115, 130402 (2015)] and permutation blocking path integral
6 Monte Carlo (PB-PIMC) [T. Dornheim et al., NJP 17, 073017 (2015)], allows for the accurate
1 computation of thermodynamic properties of the spin-polarized uniform electron gas (UEG) over a
0 widerangeoftemperaturesanddensitieswithoutthefixed-nodeapproximation. Inthepresentwork,
2 we extend this concept to the unpolarized case, which requires non-trivial enhancements that we
n describeindetail. WecompareournewsimulationresultswithrecentrestrictedpathintegralMonte
a Carlo data [E. Brown et al., PRL 110, 146405 (2013)] for different energy contributions and pair
J distributionfunctionsandfind,fortheexchangecorrelationenergy,overallbetteragreementthanfor
9 the spin-polarized case, while the separate kinetic and potential contributions substantially deviate.
1
PACSnumbers: 05.30-d,05.30.Fk,71.10.Ca,02.70.Ss

]
l
e
- I. INTRODUCTION beeninvestigatedinourpreviousworksareofrelevancefor
r the description of e.g. ferromagnetic materials or strongly
t
s magnetized systems, they constitute a rather special case,
. Quantum Monte Carlo (QMC) simulations of fermions
t since most naturally occuring plasmas are predominantly
a are of paramount importantance to describe manifold
unpolarized. Therefore, in the present work we modify
m aspects of nature. In particular, recent experimental
both our implementations of PB-PIMC and CPIMC to
- progress with highly compressed matter [1–3] such as
d plasmas in laser fusion experiments [4–9] and solids after simulate the unpolarized UEG. So far only a single data
n set for a small system (N =14 electrons, one isotherm)
laserirradiation[10], butalsotheneedforanappropriate
o could be obtained in our previous work [25] because the
descriptionofcompactstarsandplanetcores[11–13], has
c paramagnetic case turns out to be substantially more
[ leadtoahighdemandforaccuratesimulationsofelectrons
difficult than the ferromagnetic one. Therefore, we have
inthewarmdensematter(WDM)regime. Unfortunately,
1 developed novel nontrivial enhancements of our CPIMC
theapplicationofallQMCmethodstofermionsisseverely
v algorithm that are discussed in detail. With these im-
hampered by the fermion sign problem (FSP) [14, 15]. A
7 provements, we are able to present accurate results for
7 popularapproachtocircumventthisissueistherestricted
different energies for the commonly used case of N =66
9 path integral Monte Carlo (RPIMC) [16] method, which,
unpolarized electrons over a broad range of parameters.
4 however, is afflicted with an uncontrollable error due
0 the fixed node approximation[17–20]. Therefore, until
.
1 recently, thequalityof theonly availableQMC resultsfor Sincemanydetailsofourapproachhavebeenpresented
0 the uniform electron gas (UEG) in the WDM regime [21]
in our paper I [22], in the remainder of this paper we
6 has remained unclear.
restrictourselvestoabrief,butselfcontainedintroduction
1
To address this issue, in a recent publication (paper to CPIMC and PB-PIMC and focus on the differences
:
v I, Ref. [22]) we have combined two novel complemen- arising from their application to the unpolarized UEG,
Xi tary approaches: our configuration path integral Monte comparedtothepolarizedcase. InsectionII,weintroduce
Carlo(CPIMC)method[23–25]excelsathightomedium the model Hamiltonian, both in coordinate space (IIA)
r
a density and arbitrary temperature, while our permuta- andsecondquantization(IIB)and, subsequently, provide
tion blocking path integral Monte Carlo (PB-PIMC) ap- a brief introduction to the employed QMC approaches
proach [26, 27] significantly extends standard fermionic (Sec. III), namely PB-PIMC (IIIA) and CPIMC (IIIB).
PIMC[28,29]towardslowertemperatureandhigherden- Finally, in Sec. IV, we present combined results from
sity. Surprisingly, it has been found that existing RPIMC both methods for the exchange correlation, kinetic, and
results are inaccurate even at high temperatures. potential energy (IVA) as well as the pair distribution
However,althoughthespin-polarizedsystemsthathave function (IVB). Further, we compare our data to those
from RPIMC [21], where available. While we find better
agreement than for the spin-polarized case [22, 27], there
nevertheless appear significant deviations towards lower
∗Electronicaddress: [email protected] temperature.
2
II. HAMILTONIAN OF THE UNIFORM B. Hamiltonian in second quantization
ELECTRON GAS
In second quantization with respect to spin-orbitals of
Theuniformelectrongas(“Jellium”)isamodelsystem plane waves, (cid:104)rσ |kiσi(cid:105)= L31/2eiki·rδσ,σi with ki = 2Lπmi,
of Coulomb interacting electrons in a neutralizing homo- m Z3 and σ , , the model Hamiltonian, Eq. (1),
i i
geneous background. As such, it explicitly allows one to take∈s the form ∈{↑ ↓}
study effects due to the correlation and exchange of the
(cid:88) (cid:88)
electrons, whereas those due to the positive ions are ne- Hˆ = k2aˆ†aˆ +2 w− aˆ†aˆ†aˆ aˆ +Ne2ξ, (4)
i i i ijkl i j l k
glected. Furthermore, the widespread density functional i i<j,k<l
theory (DFT) crucially depends on ab initio results for i(cid:54)=k,j(cid:54)=l
the exchange correlation energy of the uniform electron
with the antisymmetrized two-electron integrals, w− =
gas(UEG),hithertoatzerotemperature[30]. However,it ijkl
w w , where
is widely agreed that the appropriate treatment of matter ijkl− ijlk
under extreme conditions requires to go beyond ground
4πe2
state DFT, which, in turn, needs accurate results for the w = δ δ δ , (5)
ijkl L3(k k )2 ki+kj,kk+kl σi,σk σj,σl
finite temperature UEG. While the electron gas itself is i− k
defined as an infinite macroscopic system, QMC simula-
and the Kronecker deltas ensuring both momentum and
tions are possible only for a finite number of particles N.
spin conservation. The first (second) term in the Hamil-
Hence, we always assume periodic boundary conditions
tonian Eq. (4) describes the kinetic (interaction) energy.
andincludetheinteractionoftheN electronsinthemain The operator aˆ† (aˆ ) creates (annihilates) a particle in
simulationcellwithalltheirimagesviaEwaldsummation i i
the spin-orbital k σ .
and defer any additional finite-size corrections [31–33] to | i i(cid:105)
a future publication.
III. FERMIONIC QUANTUM MONTE CARLO
WITHOUT FIXED NODES
A. Coordinate representation of the Hamiltonian
Throughout the entire work, we consider the canonical
Following Refs. [27, 31], we express the Hamiltonian
ensemble, i.e., the volume V, particle number N and
(we measure energies in Rydberg and distances in units
inverse temperature β =1/k T are fixed. In equilibrium
B
of the Bohr radius a ) for N = N +N unpolarized
0 ↑ ↓ statistical mechanics, all thermodynamic quantities can
electrons in coordinate space as
be derived from the partition function
N N N
Hˆ = (cid:88) 2+(cid:88)(cid:88)e2Ψ(r ,r )+Ne2ξ , (1) Z =Trρˆ, (6)
− ∇i i j
i=1 i=1 j(cid:54)=i which is of central importance for any QMC formula-
tion and defined as the trace over the canonical density
withthewell-knownMadelungconstantξandtheperiodic
operator
Ewald pair interaction
1 (cid:88) e−π2G2/κ2e2πiG(r−s) ρˆ=e−βHˆ . (7)
Ψ(r,s) =
V G(cid:54)=0 πG2 The expectation value of an arbitrary operator Aˆ is given
π (cid:88)erfc(κr s+R) by
+ | − | . (2)
− κ2V R |r−s+R| Aˆ = Tr(Aˆρˆ) = 1Tr(Aˆρˆ). (8)
(cid:104) (cid:105) Trρˆ Z
Here R = n L and G = n /L denote the real and re-
1 2
ciprocal space lattice vectors, respectively, with the box However, for an appropriate description of fermions,
length L, volume V = L3 and the usual Ewald param- Eqs. (6) and (8) must be extended either by antisym-
eter κ. Furthermore, PB-PIMC simulations require the metrizing ρˆ ρˆ− or the trace itself [23], Tr Tr−.
evaluation of all forces within the system, where the force Therefore, it→holds →
between two electrons i and j is given by
Z =Trρˆ− =Tr−ρˆ. (9)
(cid:18) (cid:19)
F = 2 (cid:88) G sin[2πG(r r )]e−π2G2/κ2 (3)
ij V G2 i− j While defining the trace in Eq. (9) as either expression
G(cid:54)=0
does not change the well-defined thermodynamic expec-
(cid:18) (cid:19)
+ (cid:88)ri−rj +R erfc(κα)+ 2καe−κ2α2 , tation values, it does lead to rather different formulations
α3 √π of the same problem. The combination of antisymmetriz-
R
ing the density matrix and evaluating the trace in co-
with the definition α= r r +R. ordinate space is the first step towards both standard
i j
| − |
3
PIMC and PB-PIMC, cf. Sec. IIIA, but also RPIMC. All matrix, with (cid:15) = β/P, and approximate each of the P
these approaches share the fact that they are efficient factors at a P times higher temperature by the fourth-
when fermionic quantum exchange does not yet domi- order factorization [40, 41]
nate a systm, but they will become increasingly costly
towards low temperature and high density. Switching to e−(cid:15)Hˆ e−v1(cid:15)Wˆa1e−t1(cid:15)Kˆe−v2(cid:15)Wˆ1−2a1
≈
second quantization and carrying out the trace in anti- e−t1(cid:15)Kˆe−v1(cid:15)Wˆa1e−2t0(cid:15)Kˆ . (11)
symmetrized Fock space, on the other hand, is the basic ×
idea behind our CPIMC method, cf. Sec. IIIB, and, in It should be noted that Eq. (11) allows for sufficient
a different way, behind the likewise novel density matrix accuracy, even for small P. The Wˆ operators in Eq.
QMCmethod[34]. Thelatterapproachhasrecentlybeen
(11) denote a modified potential that combines the usual
applied to the the case of N = 4 spin-polarized elec- potential energy Vˆ with double commutator terms of the
trons [35], where complete agreement with our CPIMC
form
results[24]wasreported. TheseQMCapproachestendto
excel at high density, i.e., weak nonideality, and become (cid:126)2 (cid:88)N
eventually unfeasible towards stronger coupling strength. [[Vˆ,Kˆ],Vˆ]= F 2 , F = V(R), (12)
i i i
m | | −∇
Therefore, it is a natural strategy to combine different i=1
representations at complementary parameter ranges as
and, therefore, require the evaluation of all forces within
this does effectively allow to circumvent the numerical
the system, cf. Eq. (3). The final result for the PB-PIMC
shortcomings with which every single fermionic QMC
partition function is given by
method is necessarily afflicted [22, 27].
(cid:90)
1
Z = dX
(N !N !)3P
↑ ↓
A. Permutation blocking PIMC
P(cid:89)−1(cid:16)e−(cid:15)V˜αe−(cid:15)3u0(cid:126)m2F˜αDα,↑Dα,↓(cid:17), (13)
1. Basic idea
α=0
Inthissection,wewillbrieflyintroduceourpermutation with V˜ and F˜ containing all contributions of the poten-
α α
blocking PIMC approach. A more detailed description of tial energy and the forces, respectively. The exchange-
themethoditselfanditsapplicationtothespin-polarized diffusion functions are defined as
UEG can be found in Refs. [26, 27].
The basic idea behind PB-PIMC is essentially equal Dα,↑ = det(ρα,↑)det(ραA,↑)det(ραB,↑) (14)
to standard PIMC in coordinate space, e.g., Ref. [29], D = det(ρ )det(ρ )det(ρ )
α,↓ α,↓ αA,↓ αB,↓
but, in addition, combines two well-known concepts: 1)
antisymmetric imaginary time propagators, i.e., determi- and contain the determinants of the diffusion matrices
nants [36–38], and 2) a fourth-order factorization of the
density matrix [39–42]. Furthermore, since this leads to a ρα,↑(i,j)=λ−t13(cid:15)(cid:88)e−λ2tπ1(cid:15)(rα,↑,j−rαA,↑,i+nL)2 , (15)
significantly more complicated configuration space with-
n
out any fixed paths, one of us has developed an efficient
set of Metropolis Monte Carlo [43] updates that utilize with λ =(cid:112)2π(cid:15)t (cid:126)2/m being the thermal wavelength
t1(cid:15) 1
the temporary construction of artificial trajectories [26]. of a single “time slice”.
As mentioned above, we evaluate the trace within the In contrast to standard PIMC, where each permuta-
canonicalpartitionfunctionforN =N +N unpolarized tion cycle has to be explicitly sampled, we combine both
↑ ↓
electrons in coordinate representation positively and negatively signed configuration weights
in the determinants both for the spin-up and spin-down
1 (cid:88) (cid:88) electrons. This leads to a cancellation of many terms
Z = sgn(σ ) sgn(σ )
↑ ↓
N !N ! and, consequently, a significantly increased average sign
↑ ↓
σ↑∈SN↑σ↓∈SN↓ in our Monte Carlo simulations. Yet, this “permutation
(cid:90)
× dR (cid:104)R|e−βHˆ |πˆσ↑πˆσ↓R(cid:105) , (10) bmloeackninign”teirs-poanrlyticelffeecdtiisvteanwchee,ni.λe.t,1(cid:15)wihsecnomthpearreabalreetoboththe
large diagonal and off-diagonal elements in the diffusion
with πˆσ↑,↓ being the exchange operator that corresponds matrices. Withanincreasingnumberofhigh-temperature
tSoNa↑,↓pawrittihcualassroeclieamteedntsiσg↑n,↓sfgrno(mσ↑t,↓h)eapnedrm↑u(t↓a)tiodnengortoiunpg foancltyobrsutPa,sλint1g(cid:15)ledleacrrgeeaseelesmaenndt,ineveeancthuraollwy,owfthheenρtαh,↑e,rethies
spin-up (spin-down) electrons. However, since the kinetic average sign converges towards that of standard PIMC.
andpotentialcontributionstotheHamiltonian,Kˆ andVˆ, For this reason, it is crucial to combine the determinants
do not commute, the low-temperature matrix elements from the antisymmetric propagators with a higher order
of ρˆare not known. To overcome this issue, we use the factorization of the density matrix, cf. Eq. (11). It is only
common group property ρˆ(β)=(cid:81)P−1ρˆ((cid:15)) of the density thiscombinationwhichallowsforsufficientaccuracywith
α=0
4
P=2 a) P=2
-5.14 P=3 0.003 P=3
P=4 P=4
CPIMC P=5
0.002
-5.148
y V|
V/R -5.156 ΔV/| 0.001
0
-5.164
-0.001
-5.172
0 0.05 0.1 0.15 0.2 0.25 0.1 1 10
t0 rs
0.0015
Figure 1: Influence of the relative interslice spacing t on the
0 b)
convergence–ThepotentialenergyfromPB-PIMCsimulations
0.001
ofN =4unpolarizedelectronsatθ=0.5andr =1isplotted
s
versus t for the fixed choice a =0.33.
0 1
0.0005
K
as few as two or three propagators while, at the same K/ 0
Δ
time, the benefit of the blocking within the determinants
-0.0005
is maximized. Furthermore, we note that electrons with
different spin-projections do not exchange at all. There-
-0.001
fore, PB-PIMC simulations of the unpolarized UEG with
N = N +N do suffer from a significantly less severe
↑ ↓ -0.0015
sign problem than for N =2N spin-polarized electrons.
↑ 0.1 1 10
rs
2. Application to the unpolarized UEG
Figure 2: Density dependence of the relative time step error
The accuracy of our PB-PIMC simulations crucially from PB-PIMC with a =0.33 and t =0.14 – The relative
1 0
depends on the systematic error due to the employed differences between PB-PIMC results with P =2,3,4,5 and
higher order factorization [26, 27]. Thus, we begin the reference data from CPIMC are plotted versus rs for the
investigation of the unpolarized electron gas with the potential energy (a) and the kinetic energy (b).
analysis of the convergence behavior with respect to the
two free parameters from Eq. (11), namely a (weighting
0
thecontributionsoftheforcesondifferenttimeslices)and from CPIMC. The statistical uncertainty is mainly due
t (controllingtherelativeinterslicespacing). InFig.1,we to PB-PIMC, except for r =4 where the CPIMC error
0 s
seta =0.33fixed,whichcorrespondstoequallyweighted bar predominates. For the kinetic energy, even for P =3
0
forces on all slices, and plot the potential energy for there are no clear systematic deviations from the exact
P =2,3,4overtheentiret -rangeforabenchmarksystem result over the entire r -range. Only with two propaga-
0 s
of N = 4 unpolarized electrons at r = 1 and θ = 0.5. tors, our results for K appear to be slightly too large
s
Evidently,forallt valuesV convergesmonotonicallyfrom for r (0.5,1,2), although this trend hardly exceeds
0 s
above towards the exact result, which has been obtained ∆K/K∈=5 10−4. For the potential energy, the factor-
×
withCPIMC.Theoptimumvaluefort islocatedaround ization error behaves quite differently. For r 1, even
0 s
≥
t = 0.14, where all three PB-PIMC values are within with two propagators the accuracy is better than 0.1%,
0
single error bars with the black line. For completeness, while towards higher density (r < 1), the convergence
s
we mention that this particular set of the optimum free significantly deteriorates. In particular, at r =0.25 even
s
parameters for the energy is consistent with the previous with P =5 there is a deviation of ∆V/V 0.1%. This
| |≈
findings for different systems [26, 27, 41]. observation is in striking contrast to our previous inves-
A natural follow-up question is how the convergence tigation of the polarized UEG, where the relative error
with P behaves with respect to the density parameter in both K and V decreased towards r 0. The reason
s
→
r . In Fig. 2, we show results for the relative error of the for this trend lies in the presence of two different particle
s
potential (∆V/V , panel a) and kinetic energy (∆K/K, species which do not exchange with each other, namely
| |
panel b), where the reference values are again obtained N spin-up and N spin-down electrons. Even at high
↑ ↓
5
density,twoelectronsfromthesamespeciesareeffectively
1
separated by their overlapping kinetic density matrices
that cancel in the determinants, which is nothing else
than the Pauli blocking. Yet, a spin-up and a spin-down
0.75
electron do not experience such a repulsion and, at weak
coupling (small r ), can be separated by much smaller
s
distances r from each other. With decreasing r the force g 0.5
terms in Eq. (11) that scale as F(r) 1/r2 will eventu-
∝ a)
ally exceed the Coulomb potential V(r) 1/r, i.e., the
higher order correction predominates. T∝his trend must 0.25 P=2
P=3
be compensated by an increasing number of propagators
P=4
P. Hence, the fermionic nature of the electrons that P=5
manifests as the Pauli blocking significantly enhances the 0
performance of our factorization scheme, which means 0 0.2 0.4 0.6 0.8 1
that the simulation of unpolarized systems is increasingly r/rs
hampered towards high density. In addition to the Monte
Carlo inherent sign problem, this is a further reason to
1
combine PB-PIMC with CPIMC, since the latter excels
just in this regime.
In our recent analysis of PB-PIMC for electrons in 0.75
a 2D harmonic trap [26], it was found that, while the
combination a = 0.33 and t = 0.14 (parameter set
0 0
a) is favorable for a fast convergence of the energy, it g 0.5
does not perform so well for other properties like, in
that case, the density profile. To address this issue, we
0.25
again simulate a benchmark system of N =4 unpolarized
electrons and compute the pair distribution function g(r),
bb))
see, e.g. Ref. [44] for a comprehensive discussion. In
0
Fig. 3, we show results for the above combination of free
0 0.2 0.4 0.6 0.8 1
parameters (a) and P =2,3,4,5. Panel a) displays the
data for the inter-species distribution function g . We r/rs
↑↓
note that, for the infinite UEG, this quantity approaches
unity at large distances, but the small simulation box for
Figure 3: Convergence of the pair distribution function for
N =4restrictsustothedepictedr-range. Allfourcurves
N =4 unpolarized electrons at θ=1 and r =4 – Shown are
deviate from each other for r (cid:46)0.5, which indicates that PB-PIMCresultsfortheinter-(g ,panelas)andintra-species
↑↓
g↑↓ isnotyetconvergedevenforP =5atsmalldistances, (g↑↑, panel b) distribution function for different numbers of
andareequalotherwise. Thisisagainaclearindicationof propagators P and the fixed free parameters a = 0.33 and
1
the shortcomings of our fourth-order factorization, which t0 =0.14.
overestimatestheCoulombrepulsionatshortranges. The
intra-species distribution function g = g , which is
↑↑ ↓↓
shown in panel b), does not exhibit such a clear trend
(b) is the choice a = 0, which means that the forces
0
since only the green curve that corresponds to P = 2
are only taken into account on intermediate time slices.
can be distinguished from the rest. This is, of course,
Due to the diagonality of the pair distribution function
expected and a consequence of the Pauli blocking as
in coordinate space, it is measured exclusively on the
explained above.
main slices, for whose distribution the force terms do
Evidently, our propagator with the employed choice not directly enter. For this reason, the inter-species pair
of free parameters (a) does not allow for an accurate distribution function is not as drastically affected by the
description of the Coulomb repulsion at short distances. divergence of the F(r) 1/r2 terms at small r and the
∝
Tounderstandthisissue,werepeatthesimulationswitha convergenceofthisquantityissignificantlyimproved. For
differentcombinationa =0andt =0.04(parameterset completeness, in panel b) we again show results for g ,
0 0 ↑↑
b), which has already proven to be superior to parameter which, for parameter set (b), are almost converged even
set (a) for the radial density in the 2D harmonic trap. fortwopropagators. Itisimportanttonotethatwhilethe
The results are shown in Fig. 4 for different numbers of description of the Coulomb repulsion at very short ranges
propagators. The data with P = 2 are nearly equal to is particularly challenging, this does not predominate
the results from parameters (a) and P =5. The data for in larger systems since the average number of particles
P = 4 and P = 5 almost coincide and are significantly within distance r [r˜,r˜+∆r˜) increases as N(r˜) r˜2.
∈ ∝
increased with respect to the other curves. The main For N =66 unpolarized electrons, which is the standard
reason for the improved convergence of parameter set system size within this work, these effects are by far not
6
1 1
0.75
0.1
g 0.5
>'
aa)) s
<
0.25 a1=0.33, P=5
a1=0.0, P=2 0.01
a1=0.0, P=4 θ=0.75
a1=0.0, P=5 θ=1.0
0 θ=2.0 PIMC (θ=1)
0 0.2 0.4 0.6 0.8 1 θ=4.0
r/rs θ=8.0
0.001
0.1 1 10
1
r
s
0.75 Figure 5: Average sign for PB-PIMC simulations of N =66
unpolarizedelectronsatdifferenttemperatures–AllPB-PIMC
data have been obtained for P = 2 with a = 0.33 and
0
g 0.5 t0 =0.14 and the standard PIMC data (red curve) have been
taken from Ref. [21].
0.25
degeneracy. With increasing r , coupling effects become
bb)) s
more important, which leads to a stronger separation of
0
0 0.2 0.4 0.6 0.8 1 the electrons. Thus, there is less overlap of the kinetic
densitymatricesandthedeterminantsbecomeexclusively
r/rs
positive. For θ =1, the average sign already significantly
deviates from unity at r =40 and exhibits a more severe
s
Figure4: onvergenceofthepairdistributionfunctionforN =4 decrease towards smaller rs. Nevertheless, it attains a
(cid:48)
unpolarizedelectronsatθ=1andrs =4–Shownisthesame finite value (cid:104)s(cid:105) ≈ 0.01 even at high density rs = 0.1,
information as in Fig. 3, but for a different combination of which means that simulations are more involved but still
free parameters, i.e., a =0 and t =0.04. manageable over the entire coupling range. This is in
0 0
starkcontrasttostandardPIMCwithoutthepermutation
blocking (red circles), for which the sign exhibits a sharp
as important and, for the same combination of r and θ drop and simulations become unfeasible below r 5.
s s ≈
as in Fig. 4, both the inter- and intra-species distribution Finally, the green curve corresponds to θ = 0.75 where
function are of much higher quality, cf. Fig. 12. PB-PIMCiscapabletoprovideaccurateresultsforr 3.
s
≥
Uptothispoint,onlydataforsmallbenchmarksystems
with N = 4 electrons have been presented. To obtain
meaningful results for the UEG, we simulate N = 66
unpolarized electrons, which is a commonly used model B. Configuration PIMC
system since it corresponds to a closed momentum shell
and, therefore, is well suited as a starting point for an
1. Basic idea
extrapolation to the thermodynamic limit (finite size
corrections). In Fig. 5, the average sign, cf. Eq. (20), is
plotted versus the density parameter r for five different In this section, the main aspects of our CPIMC ap-
s
temperatures. Forθ =2,4,8, s (cid:48) isalmostequaltounity proachareexplained. AdetailedderivationoftheCPIMC
for r = 40 and decreases ju(cid:104)st(cid:105) a trifle towards higher expansionofthepartitionfunctionandtheutilizedMonte
s
density, until it saturates at r 0.5. Consequently, Carlo steps for the polarized UEG can be found in
s
simulationsarepossibleovertheen∼tiredensityrangewith Refs. [22, 24].
relatively small computational effort. The slight increase For CPIMC, instead of evaluating the trace of the
(cid:48)
of s around r [1,10] is a nonideality effect: At high partition function Eq. (6) in coordinate representa-
s
(cid:104) (cid:105) ∈
density, the system is approximately ideal and the Fermi tion, we switch to second quantization and perform the
temperature θ is an appropriate measure for quantum trace with anti-symmetrized N particle states (Slater-
F
−
7
determinants) n(2) = 001100110...
|{ }i | i
n = n ,n ,... , (16) 8
|{ }(cid:105) | 1 2 (cid:105) ↑
7
w0it,h1 n)iofbtehiengi-tthhespfienr-moribointaicl okcicσuip,awtihoenrenwumecbheoros(enith∈e i56 s1 =(2,5,0,3) ↑↓spinp
{order}ing of orbitals such that e|ven(cid:105)(odd) orbital numbers al ↓ro
have spin-up (spin-down) σ = ( ). In this represen- bit4 ↑jec
tation, fermionic anti-symmetry↑is↓automatically taken or3 ↓tio
2 n
into account via the anti-commutation relations of the ↑σ
1 i
creation and annihilation operators, and thus, an explicit ↓
0
anti-symmetrizationofthedensityoperatorisnotneeded. ↑
The expansion of the partition function is based on the 0 τ1 τ2 τ3 β
concept of continuous time QMC, e.g., Refs. [45, 46], imaginary time τ
where the Hamiltonian is split into a diagonal and off-
diagonalpartHˆ =Dˆ+Yˆ withrespecttothechosenbasis.
Figure 6: Typical closed path of N = 4 unpolarized parti-
Summing up the entire perturbation series of the density
cles in Slater determinant (Fock) space. The state with four
operator e−βHˆ in terms of Yˆ finally yields occupied orbitals |k ↑(cid:105),|k ↓(cid:105),|k ↓(cid:105),|k ↑(cid:105) undergoes a two-
0 1 3 6
particleexcitations attimeτ replacingtheoccupiedorbitals
1 1
(cid:88)∞ (cid:88) (cid:88) (cid:90)β (cid:90)β (cid:90)β |k0 ↑(cid:105),|k3 ↓(cid:105) by |k2 ↑(cid:105),|k5 ↓(cid:105). Two further excitations occur
Z = dτ1 dτ2... dτK (17) at τ2 and τ3. The states at the “imaginary times” τ = 0
and τ =β coincide. In addition, the total spin projection is
KK=(cid:54)=01, {n}s1...sK−1 0 τ1 τK−1 conserved at any time. All possible paths contribute to the
(−1)Ke−i(cid:80)K=0D{n(i)}(τi+1−τi)(cid:89)K Y{n(i)},{n(i−1)}(si), partition function Z, Eq. (17).
i=1
with the Fock space matrix elements of the diagonal and to the modulus weights, and s measures the sign of each
off-diagonal operator path. Therefore, s (cid:48) is the average sign of all sampled
(cid:104) (cid:105)
(cid:88) (cid:88) paths during the MC simulation. It is straightforward
D = k2n(i)+ w− n(i)n(i) , (18)
{n(i)} l l lklk l k to show that the relative statistical error of observables
l l<k computed according to Eq. (20) is inversely proportional
Y{n(i)},{n(i−1)}(si)=ws−i(−1)αsi . (19) teoxpthecetaavteiornagveasliugens.cAansabceoonbsetaqiuneendceif,itnhepraavcetricaeg,erseilgianblies
Here, s = (pqrs) defines the four occupation numbers larger than about 10−4.
i
in which n(i) and n(i−1) differ, where it is p<q and
{ } { }
r < s. In this notation, the exponent of the fermionic
phase factor is given by
α =α(i) =(cid:88)q−1n(i−1)+(cid:88)s−1n(i) . 2. Application to the unpolarized UEG
si pqrs l l
l=p l=r
The difference between CPIMC simulations of the po-
Duetothetrace, eachsummandinEq.(17)fulfills n = larizedandunpolarizedUEGentersbasicallyintwoways.
{ }
n(0) = n(K) and hence can be interpreted as a β- First,inadditiontotheparticlenumberN,thetotalspin
{ } { }
periodic path in Fock space. An example of such a path projection in the summation over the starting determi-
for the case of an unpolarized UEG is depicted in Fig. 6. nant n(0) inEq.(17)hastobefixed,i.e.,thenumberof
{ }
The starting determinant n at τ = 0 undergoes K spin-up N and spin-down electrons N . Thus, if a whole
{ } ↑ ↓
excitations of type si at time τi, which we refer to as occupied orbital is excited during the MC procedure (for
”kinks”. The weight of each path is computed according details see Ref. [24]), it can only be excited to an orbital
to the second line of Eq. (17), which can be both positive with the same spin projection. For example, orbital 6 in
andnegative. SincetheMetropolisalgorithm[43]canonly Fig. 6 could only be excited to orbital 8 or some higher
be applied to strictly positive weights, we have to take unoccupied orbital with spin up (not pictured). More-
the modulus of the weights in our MC procedure and over, when adding a kink or changing two kinks via some
compute expectation values according to two-particle excitation, it is most effective to include spin
conservation in the choice of the four involved orbitals,
Os (cid:48)
O = (cid:104) (cid:105) , (20) sinceallotherproposedexcitationswouldberejecteddue
(cid:104) (cid:105) s (cid:48) to a vanishing weight.
(cid:104) (cid:105)
whereOisthecorrespondingMonteCarloestimatorofthe For the second aspect, we have to explicitly consider
observable, (cid:48) denotestheexpectationvaluewithrespect the modulus weight of some kink s = (pqrs), which is
i
(cid:104)·(cid:105)
8
100 5.26
a) N =4 5.24 CPIMC
N =14 extrapolation
10−1 N =66 E 5.22 PB-PIMC
s0hi10−2 upnolp.ol. 55..1280 a)
10−3 100 b)
s0i10−1
10−4 h
104
10−2
103 b) 15 c)
10
102 K0i
h 5
K0i 101
h 0
0 0.1 0.2 0.3 0.4 0.5
100
1/κ
10−1
Figure 8: Convergence of a) the internal energy, b) the
10−2 average sign and c) the average number of kinks with respect
10−1 100 101 102 to the kink potential parameter κ of N = 66 unpolarized
rs electronsatrs =2andθ=4inNB =88946spinorbitals. The
potenital parameter δ has been fixed to one. The blue (green)
line show a horizontal (linear) fit to the last converged points.
Figure7: Averagesigna)andaveragenumberofkinksb)ofdi-
The asymptotic value (black point) in the limit 1/κ → 0 is
rectCPIMC,plottedversusthedensityparameterforthreedif-
enclosed between the blue and green lines and, within error
ferentparticlenumbersN =4,14,66inN =2109,4169,5575
B bars, coincides with the PB-PIMC result (orange points).
plane wave basis functions, respectively, at θ=1. Shown are
the results from the simultation of the polarized (circles) and
unpolarized(dots)UEG,wherefortheunpolarizedcase2·N
B
spin-orbitals have been used.
We address this issue in Fig. 7, where we plot the
given by the modulus of Eq. (19)
average sign a) and the average number of kinks b) for
thepolarized(circles)andunpolarized(dots)UEGofN =
Y (s ) =
| {n(i)},{n(i−1)} i | 4,14and66electronsatθ =1. Comingfromsmallvalues
(cid:12) (cid:12)
(cid:12) 1 1 (cid:12) of r , the average number of kinks grows linearly with r .
(cid:12)(cid:12)(kp kr)2δσp,σrδσq,σs − (kp ks)2δσp,σsδσq,σr(cid:12)(cid:12) Depsending on the particle number, at some critical valuse
− −
4πe2 of rs, it starts growing exponentially, until it eventually
· L3 δkp+kq,kr+ks , (21) turns again into a linear dependency. The onset of the
exponential growth is connected to a drop of the average
wherewehaveusedthedefinitionoftheanti-symmetrized sign due to the combinatorial growth of potential sign
two-electron integrals from Sec. IIB. If all of the involved changes in the sampled paths with increasing number of
spin-orbitalshavethesamespinprojection,theKronecker kinks. Thisbehaviorbecomesmoreextremethelargerthe
deltas due to the spin obviously equal one, and the two- particle number, both for the polarized and unpolarized
electron integrals are efficiently blocked, i.e., in most system, so that for N = 66 electrons (blue lines), the
(momentum conserving) cases it is |wp−qrs|<|wpqrs| and averagenumberofkinkssuddenlyincreasesfromlessthan
w− < w . However, if the involved orbitals have about two to a couple of hundred, which corresponds to
|diffpeqrrse|ntsp|inpqpsrro|jections,oneofthetwotermsinEq.(21) a drop of the average sign from almost one to below 10−3.
is always zero and w− = w or w− = w . However, for the unpolarized system, the critical value of
| pqrs| | pqsr| | pqrs| | pqrs|
Hence,forotherwisefixedsystemparameters,theaverage r at which the average sign starts dropping drastically
s
weight of kinks in the unpolarized system is significantly is approximately half of that of the polarized system
larger. Since the diagonal matrix elements, cf. Eq. (18), containing the same number of electrons. In practice,
are independent of the spin, there ought to be more kinks this means that for N =66 polarized electrons at θ =1
in simulations of the unpolarized system, which in turn direct CPIMC calculations are feasible up to r 0.6,
s
∼
results in a smaller sign, because each kink enters the whereas for N =66 unpolarizd electrons direct CPIMC
partition function with three possible sign changes. is applicable only up to r 0.3.
s
∼
9
3. Auxiliary kink potential 8.55 a)
In Ref. [22], it has been shown that the use of an E 8.50
auxiliary kink potential of the form
8.45
1
Vδ,κ(K)= e−δ(κ−K+0.5)+1 (22) 100 b)
significantlyextendstheapplicabilityrangeofourCPIMC s0hi1100−−21 CCPPIIMMCC VVcc ==01.e03 09
method towards larger values of r . This is achieved by −
adding the potential to the seconds line of the partition 10−403 CPIMC Vc =0.0
function Eq. (17), i.e., multiplying the weight of each c) extrapolation
30
PB-PIMC
path with the potential. Obviously, since Vδ,κ(K) → 1 K0i 20
in the limit κ , performing CPIMC simulations for
h
→∞ 10
increasing values of κ at fixed δ always converges to the
exact result. Yet, to ensure a monotonic convergence of 0
0 0.1 0.2 0.3 0.4 0.5
the energy, it turned out that the value of δ has to be
1/κ
sufficiently small. Both for the polarized and unpolarized
system, choosing δ =1 is sufficient. In fact, the potential Figure 9: Convergence of a) the internal energy, b) the
is nothing but a smooth penalty for paths with a larger average sign and c) the average number of kinks with respect
number of kinks than κ. to the kink potential parameter κ of N = 66 unpolarized
electrons at r =0.8 and θ=1 in N =11150 spin orbitals.
In Fig. 8, we show the convergence of a) the internal s B
The potential parameter δ has been fixed to one. The three
energy (per particle), b) the average sign and c) the av-
curves correspond to CPIMC calculations where the kink
erage number of kinks with respect to the kink potential potential has been cut off at different values V , i.e., V (K)
c 1,κ
parameter κ of N = 66 unpolarized electrons at rs = 2 (cf. Eq. (22)) is set to zero if it takes values smaller than Vc.
and θ = 4. We have performed independent CPIMC Theblue(green)lineshowsahorizontal(linear)fittothelast
simulations for different κ, using integer values from 2 to converged red points. The asymptotic value (black point) in
17. While the energy almost remains constant for κ 10 thelimit1/κ→0isenclosedbetweentheblueandgreenlines
with a corresponding average sign larger than 0.1,≥the and, within error bars, coincides with the PB-PIMC result
(orange points).
average sign and number of kinks themselves clearly are
not converged. Further, the direct CPIMC algorithm
(without the kink potential) would give a couple of thou-
sand kinks with a practically vanishing sign. However, 4. Further enhancement of the kink potential
for the convergence of observables like the energy, appar-
ently, a significantly smaller number of kinks is sufficient.
It turns out that, in case of the unpolarized UEG, even
This can be explained by a near cancellation of all addi-
withtheuseofakinkpotentialwithδ =1,thesimulation
tional contributions of the sampled paths with increasing
may approach paths with an extremely large number of
number of kinks. For a detailed analysis, see Ref. [22].
kinks. This is demonstrated by the turquoise data points
Wegenerallyobserveans-shapedconvergenceofobserv- in Fig. 9 c), where the average number of kinks is shown
ables with 1/κ, where the onset of the cancellation and for N =66 unpolarized electrons at θ =1 and r =0.8.
s
near convergence are clearly indicated by the change in For example, at κ = 8, there are on average about 30
curvature. This allows for a robust extrapolation scheme kinks. However, increasing the penalty for paths with a
to the asymptotic limit 1/κ , which is explained number of kinks larger than κ, by increasing δ, is not a
→ ∞
in detail in Ref. [22]. An upper (lower) bound of the solution, since this would cause a non-monotonic conver-
asymptotic value is obtained by a horizontal (linear) fit gence, oscillating with even and odd numbers of κ, as has
to the last points after the onset of convergence. The been demonstrated in Ref. [22]. Therefore, we choose a
extrapolated result is then computed as the mean value different strategy which is justified by the fact that paths
of the lower and upper bounds with the uncertainty esti- with a very large number of kinks do not contribute to
mated as their difference. In Fig. 8, both, the horizontal physical observables, cf. Sec. IIIB3 and Ref. [22]: we cut
(blue line) and linear fit (green line) almost coincide due off the potential once it has dropped below some criti-
to the complete convergence (within statistical errors) of cal value V , thereby completely prohibiting paths where
c
thelastpoints. TheasymptoticCPIMCresult(blackdot) V (K)<V . If the cut-off value is too large, we again
1,κ c
perfectly agrees (within error bars) with the PB-PIMC recover an oscillating convergence behavior of the energy
result (orange dot). This confirms the validity of using with even and odd numbers of κ rendering an extrapo-
the kink potential also for the unpolarized UEG. lation difficult. This is shown by the purple data points
10
in Fig. 9 a), where the simulations have been performed 0.5
−
with Vc = 0.03 so that paths with a number of kinks θ=8 a)
0.6
largerthanκ+3areprohibited. Ontheotherhand, ifwe −
set V =10−9, so that paths with up to κ+20 kinks are θ=4
c 0.7
allowed, the oscillations vanish (within statistical errors) −
andwecanagainapplyourextrapolationscheme. Indeed, 0.8 θ=2
even with the additional cut-off the extrapolated value −
(black dot) coincides with that of the PB-PIMC simula- rs 0.9 θ=1
tion (orange dot) within error bars. In all simulations c·− θ=0.75
x
presentedbelowwehavecarefullyverifiedthatthecut-off E 1.0
−
θ=0.5
value is sufficiently small to guarantee converged results.
1.1
To summarize, as for the polarized UEG [22], the ac- −
cessible range of density parameters rs of our CPIMC 1.2 CPIMC
method can be extended by more than a factor two by − PB-PIMC
the use of a suitable kink potential, in simulations of the 1.3 RPIMC
−
unpolarized UEG as well. For example, at θ = 1 direct
CPIMC simulations are feasible up to rs 0.3, see Fig. 7, 0.8
whereas the kink potential allows us to∼obtain accurate − b)
0.9
energies up to rs = 0.8, as demonstrated in Fig. 9. In rs − θ=1
addition to the extrapolation scheme that has been intro- xc· −1
duced before for the spin-polarized case [22], we have cut E 1.1
−
off the potential at a sufficiently small value to prevent 1.2
the simulation paths from approaching extremely large −
0.1 0.2 0.5 1 2 5 10
numbersofkinks. WeexpectthisenhancementofCPIMC
r
s
to be useful for arbitrary systems. In particular, it will
allow us to further extend our previous results for the
Figure 10: Exchange-correlation energy E times r of the
polarized UEG to larger r -values. xc s
s
unpolarized N =66 particle UEG over the density parameter
r for different temperatures. In graphic a), only the best
s
resultsfromCPIMC(dots)orPB-PIMC(crosses)calculations
IV. COMBINED CPIMC AND PB-PIMC
are shown, cf. table I in the appendix. In addition, RPIMC
RESULTS
resultsbyBrownetal.[21,50]areplottedforcomparison(lines
with light colors and open circles). Graphic b) also shows
A. Exchange correlation energy PB-PIMC data for r <1 at θ=1.
s
The exchange-correlation energy per particle, E , of
xc
the uniform electrons gas is of central importance for is subtracted, which, obviously, means that the compara-
the construction of density functionals and, therefore, tively small remainder is afflicted with a larger statistical
has been the subject of numerous previous studies, e.g., uncertainty.
Refs. [21, 22, 25, 47–49]. It is defined as the difference To illustrate the overlap between PB-PIMC and
betweenthetotalenergyofthecorrelatedsystemandthe CPIMC, we show all available data points for θ = 1
ideal energy U0, for both methods in panel b). This is the lowest temper-
ature for which this is possible and, therefore, the most
E =E U . (23)
xc − 0 difficult example, because the systematic propagator er-
In Fig. 10 a), we show results for this quantity for six ror from PB-PIMC at small rs is most significant here.
different temperatures in dependence on the density pa- Evidently, both data sets are in excellent agreement with
rameter r . All data are also available in Tab. I in the each other and the deviations are well within the error
s
appendix. In order to fully exploit the complementary bars. Although we do expect that the deterioration of
natureofourtwoapproaches, wealwayspresentthemost the convergence of the PB-PIMC factorization scheme for
accurate data from either CPIMC (dots) or PB-PIMC small rs, cf. Fig. 2, should become less severe for larger
(crosses). This allows us to cover the entire density range systems, anysystematictrendismaskedbythesignprob-
forθ 1,sincehere,thetwomethodsallowforanoverlap lem anyway and cannot clearly be resolved for the given
with≥respect to r . For completeness, we mention that statistical uncertainty.
s
the apparently larger statistical uncertainty for θ =8 in Let us now consider temperatuers below θ = 1. For
comparison to lower temperature is not a peculiar mani- θ = 0.75, CPIMC is applicable only for r 0.7, while
s
≤
festation of the FSP, but, instead, an artifact due to the PB-PIMC delivers accurate results for r 3. Thus, the
s
≥
definition(23). Athightemperature, thesystembecomes intermediate regime remains, without further improve-
increasingly ideal and, therefore, the total energy E ap- ments, out of reach and, for θ = 0.5, PB-PIMC is not
proaches U . To obtain E at θ = 8, a large part of E applicableforN =66unpolarizedelectronsinthisdensity
0 xc