Twist-averaged Boundary Conditions in Continuum Quantum Monte Carlo

12 Pages·0.37 MB·English

by  C. Lin

Twist-averaged Boundary Conditions in Continuum Quantum Monte Carlo
C. Lin, F.-H. Zong and D. M. Ceperley
Dept. of Physics and NCSA, University of Illinois at Urbana-Champaign, Urbana, IL 61801
We develop and test Quantum Monte Carlo algorithms which use a“twist” or a phase in the
1
0 wave function for fermions in periodic boundary conditions. For metallic systems, averaging over
0 thetwistresultsinfasterconvergencetothethermodynamiclimitthanperiodicboundaryconditions
2 for properties involving the kineticenergy with the same computational complexity. We determine
exponentsfortherateofconvergencetothethermodynamiclimit forthecomponentsoftheenergy
n
ofcoulomb systems. Weshowresultswithtwist averagedvariational MonteCarloonfreeparticles,
a
J theStonermodelandtheelectrongasusingHartree-Fock,Slater-Jastrow, three-bodyandbackflow
wavefunction. We also discuss the use of twist averaging in the grand canonical ensemble, and
2
2 numerical methods to accomplish thetwist averaging.

] PACS Numbers: 02.70.-c, 82.20.Wt, 71.15.-m
h
c
Almost all quantum Monte Carlo (QMC)calculations can be restricted to be in the range:
e
m inperiodicboundaryconditionshaveassumedthatphase
of the wavefunction returns to the same value if a parti- −π <θ ≤π. (2)
- i
t clegoesaroundtheperiodicboundariesandreturnstoits
a
original position. However, with these boundary condi-
t Forsystemswitharealpotential(e.g. nomagneticfield),
s tions,delocalizedfermionsystemsconvergeslowlytothe
. onecanfurtherrestrictthetwisttobeintherange[0,π].
t thermodynamic limit because of shell effects in the fill-
a For a degenerate Fermi liquid, finite-size shell effects
m ing of single particle states. In this paper we explore an
aremuchreducedif the twistangle is averagedover. We
alternative boundary condition: one can allow particles
- call this twist averaged boundary conditions (TABC)5.
d to pick up a phase when they wrap around the periodic
This is particularly important in computing properties
n boundaries,
that are sensitive to the single particle energies such
o
c Ψ(r +Lx,r ,···)=eiθxΨ(r ,r ,···). (1) as the kinetic energy and the magnetic susceptibility.
1 2 1 2
[ By reducing shell effects, much more accurate estima-

1 The boundary cobndition θ = 0 is called periodic bound- tions of the thermodynamic limit of these properties can
v aryconditions(PBC),θ =π anti-periodicboundarycon- be made. What makes this even more important is
9 ditions (ABC) and the general condition with θ 6= 0, that the most accurate quantum methods have compu-
3 twisted boundary conditions (TBC)1. tational demands which increase rapidly with the num-
3 In periodic boundary conditions, the Hamiltonian is ber of fermions. Examples of such methods are exact
1 invariantwith respect to translating any particle around diagonalization6 (exponentialincreaseinCPUtime with
0
the periodic boundaries. According to Bloch’s theorem, N), variational Monte Carlo7(VMC) with wavefunctions
1
0 this implies that any solution can be characterized by having backflow and three-body terms8,9 (increases as
/ a given twist angle. The twist angle also has a physi- N4), and transient-estimate and released-node Diffusion
at cal origin: consider a toroidal geometry. One can either Monte Carlo methods10 (exponential increase with N).
m rotate the torus2 and go into rotating coordinates, or Methodswhichcanextrapolatemorerapidlytothether-
- add a magnetic flux3 to the center of the torus. The modynamic limit are crucial in obtaining high accuracy.
d physicalpropertieswill be unchanged. Inboth cases one Twist averaging is especially advantageous for stochas-
n can transform away the perturbation by applying TBC tic methods (i.e. QMC) because the averaging does not
o
with the twist angle given by θ = mR2ω/h for rotation necessarily slow down the evaluation of averages,except
c
: and θ = eφ/(c¯h) for magnetic flux. A torus is topologi- forthe necessityofdoingcomplexratherthanrealarith-
v cally equivalentto periodic boundaryconditions, sothat metic.
i
X a non-zero twist will be allowed in periodic boundaries. The use of twisted boundary conditions is common-
r The twist is a degreeof freedom, or boundary condition, place for the solution of the band structure problem for
a that can be varied to enable a finite system to approach a periodic solid. Band structure methods begin by as-
the thermodynamic limit more quickly or to make de- suming the wavefunction factors into single particle or-
tailed studies of the properties of the quantum state. bitals characterized by a lattice momentum. Then in
If the periodic boundaries are used in all three direc- orderto calculate properties of an infinite periodic solid,
tions, each dimension can have an independent twist4. properties must be averagedby integrating overthe first
Hence, in 3D, the twist is a three component vector, Brillouin zone. Baldereschi11 pointed out that in an in-
θ with i = {1,2,3}. The free energy and hence all sulator, in integrating over the Brilliouin zone, one can
i
equilibriumpropertiesare (triply)periodic3 in the twist: with high accuracy replace the integral with a “special
F(θ +2π)=F(θ ) so that each component of the twist k-point.” This was generalized to a grid of k-points12.
i i
1
Twisted boundary conditions has been discussed in con- sume the simulation cell is a cube (or square in 2D) of
nection with polarization of insulators13; we do not con- side L. To satisfy the twisted boundary conditions, the
siderthathere. Theuseoftwistedboundaryconditionsis wave vectors obey:
common in the analysis of lattice models6,14. Gammel15
showed using perturbation arguments for certain lattice kn =(2πn+θ)/L (3)
modelswhyitwillconvergefastertothethermodynamic
where n is an integer vector. These states have energy
limit and applied it to calculating optical properties.
E =(h¯2/2m)k2. The ground state in the canonicalen-
Gros16 studied size effects in the Hubbard model with n n
semble consists ofthe N lowestenergystates; the many-
exact diagonalization and showed TABC gives exact re-
body wavefunction is a determinant of those states. In
sultsinthegrandcanonicalensemblefornon-interacting
thissectionwewillignorespin,sinceforanon-interacting
systems.
system,spinmodifiestheresultsonlybydoublingthede-
Though twisted boundary conditions have a long his-
generacy of each level. Figure 1 shows the occupation of
tory within quantum physics, their use in quantum
statesfor 13spin-lessfermions in2Dfor θ =0 i. e. with
Monte Carlo (QMC) has been limited. In continuum
periodic boundary conditions, and also with a non-zero
QMC, Rajagopal17 used special k-points to reduce fi-
twist. The occupied states lie within a circle centered at
nitesizeerrorsforcalculationsofinsulatorswithinVMC.
the origin with radius ≈k =2(πρ)1/2.
Their use in Diffusion Monte Carlo (DMC) results were F
restricted to PBC and ABC in order to work with real
wavefunctions. Kralik et al.18 used generalized bound-
ary conditions to compute the momentum distribution
with VMC for silicon (a semiconductor) and Filippi and
Ceperley19 forlithium (a metal). This wasdone inorder
toenlargethenumberofmomentumvectorstoprobethe
behaviornearthefermisurface. TheuseofTABCwithin
QMC has not been further developed.
The applications considered here are for systems with
a fermi surface, i.e. metals. We begin by discussing
the method for non-interacting (NI) fermions. Fermi
liquid theory asserts that the spectrum of states of in-
teracting liquids are intimately related to those of the
non-interacting fermion systems, hence a detailed anal-
ysis for the NI system carries over to strongly interact-
ing fermi liquids. We then discuss interacting systems
inthe Hartree-Fockapproximation: the electrongasand
the Stoner model. In the Stoner model, we show how
TABC can be used to determine a polarization phase
transition at non-zero temperature. Results for TABC
are given for the interacting electron gas using a pair
product and backflow wavefunction in 3d. The electron
gas system has been previously treated with an extrap-
olation method based on Fermi liquid theory. We show
thatTABCgivesthesameresultsinthe thermodynamic
limit and verify the applicability of the NI analysis, in FIG.1. Momentum distibution for 13 spinless fermions in
particular to examine how the energy depends on the a 2D square with side L= 2π. The top panel shows the oc-
twist of a given system size. We then present VMC re- cupied states (closed symbols) and empty states (open sym-
sults of the polarization energy of the electron gas using bols) with zero twist (circles, PBC) and a twist equal to
the new method. In future publications we will study 2π(0.3,0.15) (triangles). The circle shows the infinitesystem
the low density properties of the electron gas using this fermi surface. The bottom panel shows the occupied states
technique. The Appendix discusses some details arising with TABC. The colored regions show the occupied region
in the implementation of TABC in QMC. forthelowest level(middlesquare),thethirdlevel,uptothe
outermost 13th level.
Figure 2 shows the relative error in energy versus the
I. NON-INTERACTING FERMIONS
number of fermions with PBC. The energy converges
slowly to the exactresult. One sees “cusps” in the curve
In a non-interacting homogenous system with PBC, atcertainvaluesofN. Theseoccuratclosed-shellvalues
the single particle states are plane waves: exp(ikr)η(σ) ofN,e.g. thestatedepictedinFig. 1forPBCisaclosed
where η is a spin function. For simplicity, we always as- shellsincestatesrelatedbysymmetryareeitherallfilled
2
or all empty. For large N the curve is “quasi-random”, onaplanebisectingandperpendiculartothelinejoining
with an envelope decaying algebraically as N−ν. the originwithanintegervector;preciselythe Lauecon-
We find numerically that the exponent of the decay of dition for the Bragg planes20 of the reciprocal lattice of
therelativeerroroftheenergyisapproximatelyν =1.33 the supercell. In fig. 3 is shown the dependence of total
in 2D and ν =1 in 3D (see Table I). To characterize the energyonthe twistanglefor afixednumber ofparticles.
approach to the thermodynamic limit, we introduce two OneseescuspsasthefilledstatescrosstheBraggplanes.
different measures. Defining the relative scaled error δ Thedependenceissimilartothebandenergyofasimple
N
as: metal. Later, we will discuss this band structure for an
interacting system. The bandwidth E (the spread of
E =E (1+N−νδ ), (4) BW
N ∞ N energy values fig. 3) is defined as:
we define a=max|δ |, b=hδ i and c=h(δ −b)2i1/2.
N N N
Table I shows estimates of these coefficients and ex- E2 =(2π)−d dθ(E(θ)−E )2 (5)
ponents obtained numerically by examining values of BW Z ∞
10≤N ≤104.
dependsonthenumberofparticlesandscalesasE ∝
BW
N−ν where the exponent is the same as describes the
convergence of the KE in PBC.
There are several alternative procedures by which the
twist angle can be varied: i) one can average the twist
over all possible values, ii) the twist can become a dy-
namicalvariableandiii) specialvalues ofthe twistcould
be used. Of these approaches, none are right or wrong
in general; which method approaches the thermody-
namic limit faster depends on the order of the phase
in question, whether fermi liquid, ferromagnetic or anti-
ferromagnetic. However, to compute a variety of prop-
erties for a metallic systems, we find the TABC best re-
duces size effects.
n O d ν a b c
1 T 2 1.33 4.5 0.37 1.77
8 T 2 1.5 0.47 0.27 0.093
1 N 2 0.67 2.18 0 0.65
16 N 2 0.75 0.71 0 0.52
1 V 2 1 0.50 -0.35 0.069
8 V 2 1 0.38 -0.367 0.0058
1 T 3 1 2.4 0.25 1.0
8 T 3 1.33 0.50 0.292 0.065
16 T 3 1.33 0.35 0.21 0.06
FIG.2. Relative error of theenergy versusnumberof par- 32 T 3 1.33 0.35 0.19 0.06
ticleswithPBC(△)andTABC(2)in2Dand3D.Thepoints 1 N 3 0.55 2.97 0 1.00
shownareonlythosewheretherelativeerrorhasalocalmax- 16 N 3 0.67 0.83 0 0.63
imum. Curves are shown only for N <100. 1 V 3 0.67 0.742 -0.549 0.072
16 V 3 0.67 0.587 -0.582 0.0043
Now consider twisting the boundary conditions, i. e.
TABLEI. Coefficientsoftheasymptoticdecayoftheerror
usinganonzerophase. Thisdisplacesthesetofk-vectors
intherelativeNIenergy. ν istheexponentofthedecay. The
as shown in the top panel of Fig. 1. Aside from a set of exponentshavebeendeterminedfromthenumericaldataand
specialtwists havingzero measure,the energylevels will are accurate to about 0.02. The amplitude was determined
no longer be degenerate. (When we sort the states to numerically by examining the values for 10 ≤ N ≤ 10000.
decide the filling, all states will have a different energy.) The coefficients are defined as a ≡ max(|δ |), b = hδ i,
N N
This is because inversion symmetry and rotational sym- c=h[δ −b]2i1/2. nisthenumberofphaseanglesusedforthe
N
metry through 90◦ is broken. This breaking of symme- summation in each dimension: θ = 2πi/n for i = 1,...,n.
i
tries and absence of degeneracies is a crucial difference n=1 corresponds to PBC; d is the dimensionality. O is the
between TBC and PBC. property: T, the kinetic energy; V, the Hartree-Fock poten-
Atcriticalvalues ofthe twist,whena filledandempty tial energy of the electron gas; N, the numberof particles in
state have the same energy, the occupation of the states theTA-GCE method.
changes. The condition for the degeneracy is that k lie
3
II. TWIST AVERAGING
The twist average of a property A is defined by:
b
π
hAi=(2π)−d dθhψ(R,θ)|A|ψ(R,θ)i (6)
Z
−π
b b
whereitisassumedthatthewavefunctionψ(R,θ)isnor-
malized for each θ. The momentum distribution, n(k) is
a key property to calculate for delocalized quantum sys-
tems. A discontinuity in n(k) at the Fermi surface is
responsible for the validity of Fermi Liquid Theory for
metals. The kinetic energy is the second moment of the
momentum distribution:
T =(h¯2/2m) dkk2n(k). (7)
N Z
Let us analyze the momentum distribution for NI
fermions in the canonical ensemble. For any given twist
θ, the N lowest energy states from those given by Eq.
3 are occupied. But any value of k can only be reached
FIG. 3. Dependence of the energy of NI unpolarized
by a unique combination of (n,θ) if θ is restricted by
fermions on the twist angle for N=54 in 3D. The solid line
Eq. 2. This proves that the averaged momentum dis-
shows theenergy along the(100) direction, dotted line along
tribution is a constant for states which can be reached the(110) direction anddashedlinealongthe(111) direction.
bysomecombinationof(n,θ)andzerootherwise. Hence Thecurvesarepiecewisequadratic,withacuspwhentheoc-
within TABC, the set of filled states comprises a “solid cupation of thestates changes. Werefer to ther.m.s. spread
volume” bounded by a Fermi surface. In contrast, for a of energy valuesas thebandwidth.
single twist value, the momentum distribution is a point
set. Thetotalvolumeink-spaceinsidethe Fermisurface
Figure 2 shows the convergence of the error of the ki-
is precisely (2π)dρ, just as it is in the thermodynamic
netic energy within TABC versus the number of parti-
limit,sotheconstantisdeterminedbythenormalization
cles. Numerical estimates of the relative error are given
condition:
in Table I. One sees a dramatic improvementin the con-
vergence with respect to PBC. The exponents governing
dkn(k)=1. (8) the decay rate are larger and the errors are up to two
Z
ordersofmagnitudesmallerforsystemsizesinthe range
N ≈ 100. Note that the TABC kinetic energy must ap-
As mentioned above, the Fermi surface is a subset of proachthe exact energyfrom above. This is because the
the Bragg planes. For N particles the occupied states shape of a givenvolume with the smallestmoment ofin-
comprise the union of the first N Brillouin zones20. The ertiaisasphere,sothatthedistortedshapeshowninfig.
(N +1)th electron will go in the (N +1)th zone, an area 1 has a higher energy. Also shown in Table I is depen-
formed by planes surrounding the Nth zone. Figure 1 dence of the error on the number of twist values in the
showsthemomentumdistributionof13spinlessfermions average. One needs from 16 to 32 values of θ along each
in 2D using TABC. axis to achieve the full reduction in size effects (better
than a percent accuracy in the relative error of the size
In 1D, TABC gives the exact momentum distribution
effects). In the Appendix are discussed the relative mer-
because the normalization condition determines every-
its of performing the average on a grid versus sampling
thing. In higher dimensions the fermi surface is not per-
the twist values from a uniform distribution.
fectly circular (spherical) as shown in fig. 1. However,
one can see that n(k) is much closer to a disk than the Letusnowexaminehowthepotentialenergyconverges
momentum distribution obtained with PBC. A perfect withPBCandTABC.Thiswillgiveussomeideaofhow
fermi surface (no finite size corrections) in any dimen- two particle correlations are affected by the boundary
sion,canbe obtainedbyallowingvariationsinthe parti- conditions since the potentialenergy is a particularinte-
clenumber andworkinginthe grandcanonicalensemble gral over the pair correlation function. The calculation
as we discuss below. performed below is particularly simple for a power law
4
potential, v(r)=r−n. In particular, we examine the po- TABC. For all values of N and twists the potential en-
tentialenergy ofanelectrongas(n=1)computed using ergy of the finite systems approaches that of the infinite
the NIwavefunction(Hartree-Fockapproximation). The systems from below. Twist averagingserves to make the
NItrialfunctionisvalidforhighdensitywhenthekinetic decay more regular but does not reduce its overall mag-
energydominatesthepotentialenergy. Thepotentialen- nitude which is determined by a charge interacting with
ergy (using the Ewald image potential) is conveniently the correlation hole in the nearby supercell. Similar ef-
evaluated in Fourier space as a sum over the structure fects are expected for other two particle quantities. The
factor: smootherconvergenceobtainedwithTABCshouldallow
for more accurate extrapolation to N →∞.
Nρ
V = v (S −1)+Nv (9)
k k M
2
X
k
III. GRAND CANONICAL ENSEMBLE
wherev istheMadelungenergyofachargeinteracting
M
withitselfandv istheFouriertransformoftheinterpar-
k
ticlepotential. Fora1/r potential,vk =2π(d−1)/kd−1. The use of TABC in the grand canonical ensemble
The values of k in the sum are k = 2πn/L where n is (GCE) gives the exact single particle occupations for NI
anintegervector. FortheNIwavefunction,thestructure particles,asshownbyGros16withintheHubbardmodel.
factoratwave-vectorqisproportionaltotheprobability Suppose(α,N)isalabelofthequantumstates,bothfor
that after we have displaced a filled state by q we are the number of particles and for other quantum numbers
still in a filled state. suchas the momentum and let E (θ) be the energy of
α,N
this state. Then the probability of a given state in the
1
Sq =1− Nh δ(k−k′−q)i (10) GCEisproportionaltoexp(−β(Eα,N(θ)−Nµ))whereµ
kX<k′ is the chemical potential. In the ground state, β → ∞,
theoccupiedmany-bodystatewillbetheoneminimizing
where the sum is overoccupied states and the averageis
E (θ)−Nµ. Thus N can depend on θ.
α,N
over twisted boundary conditions.
We now show that for a NI system, TABC within the
GCE gives exact single particle properties; i.e. there are
no finite size effects. Suppose the single particle energy
levels are e . Then the probability of occupying the N
k
states e ,...,e is exp(− N β[e (θ)−µ]). In the oc-
1 n k=1 k
cupation number nk basisP, this probability distribution
factorizes as:
b
n e−β[ek(θ)−µ]+(1−n ) (11)
k k
Yh i
k
b b
so the probability of state k being occupied is precisely
the fermi distribution law n = (exp(β(e −µ))+1)−1.
k k
As the twist angle is varied over its range, each momen-
tum state of the infinite system occurs precisely once.
Hence the averagedoccupationnumberis preciselywhat
it would be in the thermodynamic limit. The distortion
ofthe fermisurfacesobservedinthe lowerpaneloffig. 1
is a consequence of using the canonical ensemble.
The momentum distributionand hence the kinetic en-
ergy will be exactly equalto their infinite system values.
Other properties may have finite size corrections, only
the single particle propertiesare guaranteedto be exact.
We callthis procedurethe twistaveragegrandcanonical
ensemble (TA-GCE).
Withthis procedureone doesnothavea fixednumber
FIG. 4. Relative error in the evaluation of the potential ofparticlessinceforagiventwistandfermiwavevector,
energy for an electron gas using the Hartree-Fock wavefunc- thenumberofoccupiedstateswillvary. Thefluctuations
tionforN spinlesselectrons. ThesolidlineshowsTABC,the in the number of particles is closely related to a famous
dashed line PBC. problem in analytic number theory, “Gauss’ circle prob-
lem”, to determine the number of lattice points inside a
Shown in fig. 4 is the convergence of the potential circleofareaAasitsradiustendstoinfinity21. AsGauss
energy versus the number of particles using PBC and posed the problem, the center of the circle was fixed on
5
a lattice site while in TA-GCE the circle is placed ran- (3π2ρ)5/3
E/N =λ [(1+ζ)5/3+(1−ζ)5/3+b(1−ζ2)]
domly (the shifted circle problem). It has been shown22 10π2
that in 2D one obtains the following fluctuation in the
(13)
particle number for TA-GCE:
where b = 5gρ(3π2ρ)−2/3 and λ= h¯2/2m. For b <1.111
lim h(N −(LkF)2/4π)2i1/2 ∝N1/4. (12) the system6λhas an unpolarized ground state and for
N→∞
b > 1.3228 the ground state is ferromagnetic. For in-
termediatecouplings,the groundstate hasapartialspin
Our numerical estimates for the convergence are shown
polarizationatzerotemperature,similartothe behavior
in Table I. Note that the Gauss circle problem differs
suggested26,27 for the electron gas at low density.
from the convergence of the energy in the canonical en-
We performed a MC simulation of the 3D NI sys-
semble since the energyis the averagesecondmomentof
tem, within the occupation representation. The state
thelowestN points,notthenumberofpointsinacircle.
variables of the system consist of occupation numbers
The preceding discussion refers to NI systems, how-
(both spin and wavevector) and the twist angles. We
ever, fermi liquid theory asserts that the low lying ex-
makeMonteCarlomovesconsistingofflippingspinsand
cited states of an interacting system are in one-to-one
moving the spatial occupation while keeping the parti-
correspondence with the NI states. Hence, as argued by
cle number fixed. If the new state is not already occu-
Gross16andGammel15,TA-GCEislikelytoreducefinite
pied, the move was accepted with probability equal to
size effects substantially for interacting systems.
exp(−β(e (θ) − e (θ))). Fig. 5 shows the conver-
Foraninteractingmany-bodysystem,onedifficultyin new old
genceofmagnetizationdistributionversusthenumberof
using the GCE is the need to optimize the wavefunction
fermions. We found that one could determine the phase
at each twist value; in particular to pick out which or-
transition within a few percent accuracy in the coupling
bitalsshouldbeoccupied. Foranisotropicsystemhaving
constantb using only 54 fermions. Such accuracy within
a sphericalfermi surface, the order of filling the states is
PBCwouldrequirethousandsoffermions becauseofthe
simple. For a metal with a non-spherical Fermi surface,
strong shell effects.
the usual procedure is to determine the filling of single
particle states according to a mean-field theory such as
the Kohn-Sham method in density functional theory. If
oneusesthesameprocedurewithinQMC,theFermisur-
face will be substantially unchanged. Hence, it would be 1
better to try other fillings, choosing the one which mini-
mizes E −Nµ.
N
ThereareotherproblemsinusingtheGCEforcharged 0.8
systems in periodic boundary conditions. Usually the
positive compensating charge, either a uniform back-
ground or a fixed array of charged nuclei, is at a fixed 0.6
density. But if the number of electrons fluctuates, the
periodic cell can have a net charge, which causes prob- m
lems in calculating the long-range potential. Although
0.4
verypromising,wedonotconsiderthe TA-GCEmethod
further inthis paper. The followingexamples areforthe
canonical ensemble.
0.2
IV. THE STONER MODEL
0
1 1.5 2 2.5
b
FIG.5. Magnetization asafunctionofbforthe3DStoner
To test the utility of TABC for determining a phase model. Circlesaretheexactresults,diamondsfor54electrons
transition, we simulated the Stoner model, an analyti- with PBC, squaresfor 54 electrons with TABCand triangles
callysoluablemodel23,24forstronglycorrelatedfermions, for N=200 with TABC. The temperature was T = 0.224E
F
particularly related to itinerant magnetism25. The and a 83 grid of twist values was used.
Stoner model differs from NI fermions by the addition
of a contact repulsive potential: gδ(r ). The en-
i<j ij
ergy is evaluated within the meaPn field (Hartree-Fock) V. ELECTRON GAS
approximation using the NI wavefunction. For a spin
(1/2) system, the potential energy is E =E +gn n . We now present results for the 3-D electron gas as
NI ↑ ↓
In 3D the energy at zero temperature in the thermody- a test of TABC on a correlated many-body continuum
namic limit is: system. The electron gas is a very important model
6
in condensed matter physics being the basis for the Fig. 6 shows the convergence of the VMC energy ver-
density functional theory method of electronic structure susthe numberofparticleswithin TABC andPBC.One
computations.28 The phase diagram of the electron gas canseetheconvergencetothethermodynamiclimitthat
at low density is still not resolved27. The wavefunctions was found within the NI system is also evident within
we use are the most accurate known for a continuum VMC using the SJ wavefunction. In addition to the
many-body system with optimized 2-body, 3-body and Slater-Jastrow trial functions, we also have used opti-
back flow terms. Variational MC (VMC) and Diffusion mizedbackflow-threebody functions(BF3B)9 whichgive
(DMC) are QMC methods appropriate to zero tempera- a more accurate description of the low-density electron
ture, and Path Integral(PIMC) to T >0. In this paper, gas. (They pick up about 3/4 of the remaining corre-
we will only discuss VMC and DMC. PIMC will be dis- lation energy9 and break certain symmetries of the SJ
cussed in future publications. wavefunctions.)
In VMC, one assumes an analytic form for a trial In the Diffusion Monte Carlo (DMC) method, one
function Ψ (R) where R symbolizes the 3N coordi- startswithatrialfunctionandusesexp(−tH)toproject
T
nates. Then one samples |Ψ (R)|2 using a random out the ground state using a branching random walk.
T
walk7. An upper bound estimate to the exact ground Fermi statistics pose a significant problem to the pro-
state energy is the average of the local energy E (R) = jection method, since exact methods such as transient
L
Ψ (R)−1HΨ (R) over the random walk. An accurate estimate or release-node QMC suffer a exponential loss
T T
trial wavefunction is obtained from the NI wavefunction of efficiency for large numbers of particles. For this rea-
by multiplying by paircorrelationterms. The pairprod- sontheapproximatefixed-nodemethodisnormallyused.
uct or Slater-Jastrowwave function is: Usingthefixed-nodemethod,oneobtainsthebestupper
bound to the energy consistent with an assumed sign of
Ψ2(R)=Det(eikirj)e− i<ju(rij). (14) the wave function. Both the FN method and the exact
P
transient estimate can be generalized to treat complex-
valued trial functions.31 These methods are called the
fixed-phaseandreleasedphaseQMC.Wehavetriedboth
of these approaches using TABC.32
−0.1055
3D r=10
s
VI. FERMI LIQUID THEORY AND TABC
−0.1060
Fermi liquid theory (FLT) for metallic systems allows
both a method to extrapolate to the thermodynamic
)
yd limit and a way of understanding the twist dependence
R−0.1065
( of the QMC results. According to Landau, the low-lying
E
excitations of an interacting system are in close relation
to those of the NI system.
−0.1070
E =E0+Z dkδnke(k)+Z dkdk′δnkδnk′f(k,k′)+...
−0.1075
0.000 0.005 0.010 0.015 0.020 (15)
−1
N
FIG.6. Energy versus the number of electrons for the 3D where δnk is the deviation of the quasi-particle occupa-
electron gas using SJ wavefunction. The circles and connect- tion from the ground state, and e(k) and f(k,k′) are
ing line are TABC with 103 twist values. The filled △’s are 1 and 2 quasiparticle energy functionals. For simplicity
using PBC. The arrow shows theextrapolation toan infinite we have not indicated dependence on spin. The energy
system made usingthe PBC calculations using Eq. 16. functionalsareusuallyfurtherexpandedaboutthe fermi
surface in spherical harmonics and applied to calculate
One can determine the correlation factor u(r) either properties in the thermodynamic limit. Here we discuss
with an analytic argument or by a minimizing of the en- how to apply FLT when the “excitation” is caused by
ergy and/or the variance of an assumed form. We use a the boundary conditions. For example, we consider the
parameter-free analytic form so that the systematic size momentum distribution shown in fig. 1. The change in
and twist effects are not masked by noise in the trial the energy caused by the non-circular shape should also
function itself. For the electron gas very accurate ana- be given by Eq. 15.
lytic form for u(r) based on the energy within the RPA We cananalyze this dependence by comparingthe en-
approximation29 has as low an energy as those with op- ergies of the non-interacting system within a given twist
timized parameters. We used optimized Ewald sums30 with the interacting system. This is done in Fig. 7. One
bothforthepotentialandforthecorrelationfactorsoas sees a linear relation between the two energies, confirm-
to have the correct long wavelength behavior. ingthatforthesewavefunctions,Fermiliquidtheoryisan
7
appropriate description. The slope is the effective mass
of the quasi-particles. Then since the NI infinite energy
is known this gives a way of determining the interacting
system energy.
Previouscalculationsontheelectrongashaveusedan-
otherapplicationofFLT,the extrapolationmethod29,33.
In this method, one calculates accurate groundstate en-
ergiesforasequenceofparticlenumbers,anddetermines
the effective mass, potential correction and infinite sys-
tem energy by fitting these energies to the relation:
E =E +(m/m∗)∆T +ǫN−ν (16)
N ∞ N
where∆T isthedeviationoftheNIkineticenergyfrom
N
theinfinitesystemandν istheexponentforthepotential
energygiveninTableI.Figure6showsthattheestimate
FIG. 7. Plot of the VMC energy versus the deviation of
of the infinite system energy obtained using TABC and
the NI energy from E∞. Each phase angle is plotted sepa-
the extrapolation method agree within errors. This is
rately. Simulations are for N = 54 unpolarized 3D electron
reassuring,but expected, since both are based on FLT.
gas at r = 50. Since the excitations are linearly related,
s
Fermi liquid theory describes the phase angle dependence.
However, applied in practice, there are significant ad- UpperpointsaretheSlater-Jastrowwavefunction,lowerones
vantages to the TABC method. The fitting method re- are done with the backflow-3-body trial function. Effective
quireswellconvergedrunsforatleastthreedifferentval- masses (the slope) are respectively 1 and 0.61 for the two
ues of N. This can be difficult, for example when the trial functions.
unit cell is large. For example, suppose one is doing a
simulation of bcc hydrogen where one is limited to val-
4e-05
ues of N equal to 2K3 = {2,16,54,128,250...}. The
last 3 values are reasonable to simulate, however, if the 3D r=40 TABC=103
s
unit cell contained10 electrons,for example, the extrap-
olation, would be very difficult without further approx-
2e-05
imations. There is also the problem that convergence
and trial functions for large N may differ from that for
small N. For example, one typically determines m/m∗
0)
and ǫ within VMC and then applies the extrapolation ζ=
Ε( 0e+00
using the DMC energies. There are other problems with −
ζ)
extrapolation. For example, suppose we have a partially E(
spin polarized system with say n spin up particles and
1
n down particles. Using PBC there are no necessarily
2
larger,closedshell systems with precisely the same ratio -2e-05
ζ =(n −n )/(n +n ).
1 2 1 2
Within TABC, there is the possibility of obtaining re-
-4e-05
sults in the thermodynamic limit, using only quantities 0 0.2 0.4 0.6 0.8 1
ζ
computed for a single value of N. Figures 7-8 demon- FIG.8. Energy versus spin polarization for the 3D elec-
strate that the kinetic contribution is well controlled. tron gas at r =40 and N = 54. The dashed line is with the
s
Other methods34 are needed to make the potential cor- Slater-Jastow function, the solid line with backflow-3-body
rection also for a fixed N. Use of TABC should make trial function. Calculations are done using VMC with 103
this more accurate as well, as demonstrated in fig. 4. twist values.
8
VII. OTHER PROCEDURES USING TWISTED The special k-point method is not appropriate for a
BOUNDARY CONDITIONS metal because of the discontinuity in properties at cer-
tain twist values. One cannot replace the average by
Dynamical twist method. As an alternative to a single point because the Fermi surface is not given in
TABC, one can let the twist angle be a dynamical vari- advancebysymmetry,andcanchangebetweentheinter-
able (DTBC) and be determined self consistently. By acting and non-interacting wave functions. In addition,
dynamical is meant that the probability of a given twist itis not expectedthat the same twistvalues appropriate
is proportional to its free energy: fortheNIenergy,willbeappropriateforotherquantities
such as the potential energy, or spin susceptibility. It is
P(θ)∝exp(−βF(θ))= exp(−βE (θ)) (17) better to have a method which can give a spectrum of
α
Xα properties correctly, rather than only the kinetic energy.
As we discuss in the Appendix, TABC does not impose
where E (θ) is the energy of the state α with twist θ.
α anexcessivecomputationalburden,andistobepreferred
Thisdistributionoftwistanglescouldbeattainedwithin
over using only special twist values.
QMCbyenforcingdetailedbalanceonmovesofthetwist
angle during a random walk. The expectation of an op-
erator A will be:
VIII. CONCLUSIONS
1
hAi = dθe−βF(θ)hΨ (θ)|A(θ)|Ψ (θ)i. (18)
DTBC Z Z α α
Note that there is a significant difference in the ef-
At zero temperature, a special set of twist angles, those
ficiency for stochastic (QMC) methods versus explicit
whichhavethe lowestenergy,willbe singledout,sothat
methods such as density functional theory or exact di-
the ground state energy will be E = min (E(θ)).
DTBC θ agonalization in regards to the TABC. In explicit meth-
Clearly E ≤ E , in contrast to the TABC which
DTBC ∞ ods,computationsforeachtwistrequireanequalamount
gives an upper bound. In general,the dynamicalenergy,
ofcomputer time so that averagingoverN twist values
G
E converges as slowly to the thermodynamic limit as
DT willtakeroughlyN timesaslongasasingletwistvalue.
G
does the PBC energy: the exponents and fluctuations
Onecanuseinversionandrotationalsymmetrytoreduce
are the same. The dynamical twist method is not an
this, so that for a grid of 163 points, TABC will require
improvement over PBC for approaching the thermody-
only165twists. Ontheotherhand,inQMCallthetwists
namic limit for a metallic system at temperatures much
reduce the statistical error of the average. The twists
lower than the fermi temperatures.
are simply three more degrees of freedom on top of the
Althoughthedynamictwistmethodisnotsatisfactory
3N coordinate variables to be averaged over. However,
foraFermiliquidatlowtemperatures,forcertainlattice
one must also take into account startup costs associated
models such as an antiferromagnetic Heisenberg model,
witheachtwist,suchasre-optimizingthetrialwavefunc-
it can be a definite improvement over PBC and TABC.
tion or equilibrating the random walk. Neglecting these
This is because one can only establish a defect-free N´eel
startup costs, there is no loss of statistical efficiency in
orderingiftheunitcelliscommensuratewiththebound-
performing TABC so that the gain in reducing the sys-
aryconditions. FortheHeisenbergmodelonatriangular
tematic error is free. This is examined in more detail in
lattice,thegroundstatehasagiventwistperlatticespac-
the Appendix.
ing. UsingDTBC allowsonetoestablishthistwistvalue
Therearemanyexampleswheretwistaveragingcanef-
automatically, without imposing it in advance. This is
fect considerable improvement over the use of PBC. We
equivalenttotheclassicalvariablecellmethodwherethe
have been able to perform quite accurate calculations of
dimensions and aspect ratio of the supercell of a crystal
the polarization energy of the 2 and 3D electron gas32
become dynamical variables, so that one can determine
and of liquid 3He using backflow wavefunctions with on
the most stable crystal lattice structure35 instead of ex-
the order of 100 fermions. As pointed out by Ortiz et
amining each crystal structure explicitly.
al.27calculationsonsuchsmallsystemsinPBChavecon-
siderable systematic errors. One property that could be
Special points For eachvalue of N there exists a setof computed more accurately with TABC is the estimation
twist values for which E(θ) = E . One can determine of fermi liquid parameters by calculating particle-hole
∞
thesespecialtwistvaluesfortheNIsystemandthenper- excitation36. TABC can reduce the shell effects which
formsimulationsatonlyoneofthosetwistvaluesforthe caused much difficulty in that calculation. A related ex-
interacting system, thereby getting rid of single particle ample is in computation of the charge response of the
size effects. This is similartothe specialk-pointmethod electron gas37 where a considerable effort was made to
of Baldereschi11 for insulators where a single k-point is cancel out effects of the PBC. We are presently study-
determined by symmetry, thus allowing one to replace ingtheelectrongasconfinedtoaslab38 todeterminethe
an integral over the Brillouin zone with evaluation at a work function and surface energy of a metallic surface.
single k-point. This method was used within QMC by The filled states consist of a set of disks, each of which
Rajagopal et al.17. will have a certain occupation number. By doing twist
9
averagingwehaveshownreducedsizeeffectswithrespect earlier, even a factor of 2.5 in computer time is worth-
to PBC. whileifoneisabletoapproachthe thermodynamiclimit
Experimental systems are at a non-zero temperature. quicker,since QMC methods scale asNν with1≤ν ≤4
For NI systems one occupies the states with probability or, in the case of exact fermion methods, as exp(γN). If
givenbytheFermi-Diracdistribution. BecausetheFermi TABCsavesgoingtolargerN,theadditionaltimedoing
function at non-zero temperature is a continuous func- complex arithmetic is well justified.
tion,theconvergencetothethermodynamiclimitwillbe There are severalalternativesfor performing the twist
much faster, even with PBC. Howeverin practice, one is averaging:
interestedinelectronicsystemsclosetothegroundstate;
1. Evaluateas w E usinga griddefinedbypoints
therelevantquantityisthethermaldeBrogliewavelength i i i
of the electron: h¯/(m k T)1/2 ≈32˚A at T =300K. Be- θi with weigPhts wi (with iwi = 1) in the region
e b specified by Eq. 2. P
causethislengthisusuallylargerthanthesimulationcell
in QMC, the localization of the density matrix does not 2. Sample the twist during the QMC random walk
helpatreducingfermionfinitesizeatthesetemperatures. and take the average.
Using a non-zero temperature just to achieve faster con-
vergence to the thermodynamic limit is not practically 3. Acombinationofthetwoapproaches: workingona
useful. However, TABC is extendable to finite temper- gridthatisaugmentedwithrandomdisplacements.
ature PIMC simulations using the fixed-phase method31
As we discuss below, all three methods are satisfactory;
and will reduce size effects at low temperature.
there is no fundamental difference in efficiency. The
TheTABCmethodislikelytobevaluableforallQMC
choice of whether to sample or use a grid is primarily
calculationsinsystemswithafermisurface. Thecalcula- based on convenience and programming considerations
tions on the electron gas demonstrate that even though
and only secondly on efficiency. We note that all meth-
it may be a little slower per step of the random walk,
ods are easy to parallelize.
it is better to do TABC than a larger system with peri-
odicboundaryconditionsbecauseTABCconvergesmuch
Grid averaging. First,wemustaddressthequestionof
fastertothe thermodynamiclimit. Theoverallefficiency
which grid and integration rule to use. Since all prop-
ofanynumericalmethodisultimatelyjudgedbythecom-
erties are periodic with respect to the twist, the grid
puter time needed to reduce systematic and statistical
shouldbe a Bravaislattice with equallyweightedpoints.
error below a given value. TABC is effective in reduc-
One must keep in mind that the properties, though pe-
ing the systematic errors and thus improve the overall
riodic and continuous, have discontinuous derivatives at
efficiency.
the Bragg planes. Unless grid points can be located on
This research was supported by NSF DMR-98-02373
theseplanes(whichisdifficulttoachieveinpractice),the
and the Department of Physics at the University of Illi-
integrationerrorwillgoasǫ ∝∆θ2 ∝N−2/D whereN
nois Urbana-Champaign. We acknowledge useful dis- G G G
is the number of grid points. Numerically, we find that
cussions with R. M. Martin and G. Bauer. Computa-
163 grids are needed for an accuracy of 10−3 (see table
tional resources were provided by the NCSA. We thank
I). This slowly convergent, systematic error is the main
H. Edelsbrunner for references concerning Gauss’ circle
drawback of the grid integration method. For an insu-
problem.
lator with a large enough gap to excitations, properties
would be analytic for all θ since the occupation of single
particlestateswillnotchangeasafunctionoftwistangle
APPENDIX
and the grid error would converge exponentially fast.
Once a grid is chosen, one can use symmetry (e.g. in-
Here we discuss numerical details of implementing version and rotation through 90o) to reduce the num-
twist averaged boundary conditions. Many of the ber of grid points and give them a weight (wi with
changes caused by twisted boundary conditions arise NG w = 1) proportional to their multiplicity. It
i=1 i
fromthe need to havecomplex wavefunctions. Although iPs easy to show that the optimal amount of computer
the wavefunction and energies in special cases are real time at eachgridpoint shouldbe chosenproportionalto
(e.g. PBC and ABC), complex functions are needed for w /ζ(θ )1/2 where the MC efficiency at θ is defined as
i i i
general twist angles. There is a factor of roughly 2.5 in ζ(θ ) = 1/(var(E(θ ))cputime). We have found on the
i i
additional CPU time to do the arithmetic to evaluate calculations of the electron gas described earlier, that
the determinant and its derivatives. The actual impact this efficiencyis independent ofthe twistangleexceptat
on the total speed is smaller than this because the cal- the special PBC and ABC points where real functions
culations of two-particle quantities such as the potential can be used. Even though we have symmetry and can
energy and correlation factors are still done with real integrate over a reduced set of twist values, we must in-
arithmetic; the actual penalty of working with non-zero tegrate longer at high multiplicity points since they con-
twist depends on the number of particles and the type tribute more to the average. Hence, the symmetry does
ofthe trialwavefunction. However,aswe havediscussed notsignificantlyreducetheneededamountofCPUtime.
10