The sign problem and population dynamics in the full configuration interaction quantum Monte Carlo method

0.34 MB·English

by  J. S. Spencer

The sign problem and population dynamics in the full configuration
interaction quantum Monte Carlo method
J.S. Spencer,1,2 N.S. Blunt,2 and W.M.C. Foulkes2
1)Department of Materials, Imperial College London, Exhibition Road, London, SW7 2AZ,
U.K.
2)Department of Physics, Imperial College London, Exhibition Road, London, SW7 2AZ,
U.K.
TherecentlyproposedfullconfigurationinteractionquantumMonteCarlomethodallowsaccesstoessentially
2 exact ground-state energies of systems of interacting fermions substantially larger than previously tractable
1 withoutknowledgeofthe nodalstructureofthe ground-statewavefunction. We investigatethenatureofthe
0 signprobleminthismethodandhowitsseveritydependsonthesystemstudied. Weexplainhowcancelation
2 of the positive and negative particles sampling the wave function ensures convergence to a stochastic repre-

n sentationofthemany-fermiongroundstateandaccountsforthecharacteristicpopulationdynamicsobserved
a in simulations.
J

8
2 I. INTRODUCTION culation using the same basis. The memory required by
FCIQMC simulations still scales factorially with system

] size,buttheexponentappearstobesubstantiallysmaller
h Oneofthemajorgoalsofelectronicstructuremethods
than for FCI simulations. Moreover, unlike the iterative
p is to produce accurate ground-state energies and prop-
p- erties of many-electron systems.1 Quantum chemistry daligaogrointhamliatisiorneasdchilyempeasrarleleqluiziraebdlefaonrdFcCaIn, rtuhne eFffiCcIiQenMtlCy
m provides a hierarchy of ab initio methods2 based upon
on thousands of cores. Alavi and co-workers have used
Hartree-Fock, of which coupled-cluster singles and dou-
o FCIQMC to reproduce essentially every molecular FCI
bleswithperturbativetriples(CCSD(T)) isthemostac-
c calculationeverdoneandhaveobtainedground-stateen-
. curate applicable to medium-sized molecules.3 Quantum
s ergies for systems with Hilbert spaces many orders of
c Monte Carlo methods such as Green’s function Monte magnitude larger than the largest FCI calculations.13–15
i Carlo4,5 (GFMC), the closely related diffusion Monte
s We believe that the FCIQMC method will become in-
y Carlo6–8 (DMC), and auxiliary-field quantum Monte
creasingly important in the electronic structure commu-
h Carlo9(AFQMC)produceaccurateresultsviaastochas-
nity, especially if improved algorithms or substantially
p tic sampling of the many-electron wave function, but
cheaper approximations can be found. Our motivation
[ none of these methods is exact: GFMC and DMC sim-
for exploring the behavior of the method is to provide
3 ulations of all but the smallest systems converge to the insight into possible improvements.
v physicallyirrelevantmany-bosongroundstateunlessthe
9 fixed-node approximationis made;6,7 whilstAFQMC re- An FCI calculation based on iterative diagonalization
7 quires the phaseless approximation10 in order to avoid (using, for example, the Davidson method) requires the
4 an exponential growth in noise except in certain special storageofatleasttwovectors,eachoflengthequaltothe
5 cases.11 size of the Hilbert space. An FCIQMC simulation using
. thesamebasisrequiresthestorageofthelabelsofthede-
0 The full configuration interaction (FCI) method12
1 terminants occupied by a population of stochastic walk-
1 casts the Schr¨odinger equation as a matrix eigenvalue ers (which, following Anderson,6 we call “psi-particles”
1 problem, in which the requirement that the many- orpsips) scatteredoverthe same Hilbertspace. In order
: electronwavefunctionbeanti-symmetricwithrespectto for FCIQMC to be more efficient than FCI, the num-
v
exchange of electrons is imposed by working in a space
i ber of psips required must be a small fraction of the
X of Slater determinants formed from a finite basis set of
size of the Hilbert space. The fraction required is sys-
r single-particlewavefunctions. Thelowesteigenvalueand tem dependent and provides a measure of how “hard” it
a eigenvector of the FCI Hamiltonian matrix give the ex-
is for FCIQMC to treat that system. The hardness is
act ground-state energy and wave function of the sys-
surprisingly difficult to predict. FCIQMC is wildly suc-
tem, subject only to the error due to the finite basis set.
cessful for some systems, such as the neon atom, where
Whilst the computational cost of FCI scales factorially
it requires only 0.01% of the memory of a conven-
withsystemsize,itneverthelessrepresentsthe holygrail tionalFCI calcula∼tion.13 Evenforthe nitrogenmolecule,
of electronic structure methods.
a classic example of a strongly correlated system and
In 2009, Booth, Thom and Alavi13 introduced a new a tough test for quantum chemical methods, FCIQMC
stochastic approach in which the nodal structure of the usedonly a quarterof the memory ofthe equivalentFCI
ground-state wave function emerges spontaneously by calculation.13 However, FCIQMC struggles to describe
samplingthediscretespaceofSlaterdeterminants. Their themethanemolecule,13 forwhichHartree-Fockisavery
“full configuration interaction quantum Monte Carlo” good approximation. We show in Sec. II that FCIQMC
(FCIQMC)methodyieldsexact(i.e.FCI-quality)results also struggles when applied to the Hubbard model and
whilst requiring a fraction of the memory of an FCI cal- cannot treat systems larger than existing FCI methods
2
unless U is very small. What is it that makes a system and impose the normalization constraint using a La-
difficult? Evidentlytheanswerismorecomplicatedthan grange multiplier, the optimal coefficients c form the
0i
whether or not the system is strongly correlated. lowest eigenvector c of the matrix-eigenvalue problem
0
The aim of this paper is to understand the FCIQMC
algorithm better. Why are some systems more difficult Di Hˆ Dj c0j = Hijc0j =E0c0i, (5)
h | | i
totreatthanothers? Whatdeterminesthecharacteristic Xj Xj
population dynamics observed in FCIQMC simulations?
How does the cancelation of positive and negative psips whereE0isboththeLagrangemultiplierandtheground-
ensure convergence to the many-fermion ground state? state energy eigenvalue. This approach yields the FCI
How many psips are required to obtain correct results? wave function for the finite basis set used.
What goes wrongif the population ofpsips is too small? One way to find the ground-state eigenvector c0 of
We provide at least partial answers to all of these ques- theHamiltonianmatrixHistosolvetheimaginary-time
tions. Schr¨odinger equation:
Section II reviews the FCIQMC method and applies
dc
i
it to the Hubbard model as an example. In Sec. III = H c . (6)
ij j
dτ −
we discuss the nature of the sign problem in FCIQMC Xj
and explain the effect of canceling positive and negative
psips. SectionIVshowshowouranalysisofthesignprob- Thesolutionvectorc(τ)convergestoc0 asτ when-
→∞
lem alsoexplains the characteristicpopulationdynamics ever the starting vector c(τ=0) has a non-zero overlap
observed in FCIQMC simulations. Section V describes with c0; indeed, this is also the driving principle behind
several special cases in which the sign problem can be other projector methods such as DMC and GFMC. The
removed entirely. We offer some concluding remarks in convergence is easy to demonstrate by considering the
Sec. VI. formal solution of the imaginary-time Schr¨odinger equa-
tion:
II. FCIQMC METHOD c(τ)=e−Hτc(0). (7)
Iftheinitialvector,c(0),isexpandedintermsofthecom-
We briefly summarize the FCIQMC method here, pleteorthonormalsetofthe2MC eigenvectors, c ,of
N α
largely following the derivations given by Booth et al.,13 the Hamiltonian matrix, H, { }
withattentionpaidtodetailsrelevantlaterinthispaper.
Consider an orthonormal set of 2M single-particle c(0)= v (0)c , (8)
α α
spin-orbitals, φ1,φ2,...,φ2M . Previous FCIQMC Xα
work13–15 has u{sed a basis of H}artree-Fock spin-orbitals,
the solution of the imaginary-time Schr¨odinger equation
which is often a sensible choice but not required. One
can be written as
can construct an N-electron Slater determinant, D , by
i
selecting any N spin-orbitals (assuming N 2M):
≤ c(τ)=e−Hτ vα(0)cα = vα(0)e−Eατcα. (9)
D =D (1) Xα Xα
i i1,i2,...,iN
φ (1) φ (2) φ (N) Thesummationisdominatedbytheground-statecontri-
= √1N!(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)φii21...(1) φii21...(2) ··.··.··. φii21(...N)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12). (2) bution in the licm(τit→τ →∞)∞≈avn0d(0t)heu−sE0τc0. (10)
(cid:12)(cid:12)(cid:12)φiN(1) φiN(2) ··· φiN(N)(cid:12)(cid:12)(cid:12) The steady change in normalization due to the e−E0τ
To ensurethatalls(cid:12)uchdeterminantsareunique(cid:12)(up to a factor is awkward but can be removed by choosing the
sign), it is convenient to require that i1 <i2 < <iN. zero of energy such that E0 = 0. In practice, E0 is not
The ground-state wave function, Ψ , may be··d·efined knownuntiltheendofthesimulation,soweinsteadsolve
0
variationally as the N-electron wave function that mini-
dc
mizes the energy expectation value, i = (H Sδ )c , (11)
ij ij j
dτ − −
Xj
E[Ψ]= ΨHˆ Ψ , (3)
h | | i
where the energy shift, S, is adjusted slowly to keep the
subject to the normalizationconstraint ΨΨ =1. If we normalization more or less constant. The energy shift
h | i
restrictthe searchfor Ψ0 to the subsetofwavefunctions convergestotheground-stateeigenvalue,E0,inthelong-
thatcanbeexpandedinthebasisof2MC determinants, time limit. Note that previous work13–15 also subtracted
N
the Hartree-Fock energy, E , from the diagonal ele-
HF
Ψ = c D , (4) ments of the Hamiltonian, but this amounts merely to
0 i i
Xi a redefinition of S. To make the subsequent notation as
3
simple as possible, we define a “transition matrix”, T, boundedfrombothaboveandbelow. Italsofollowsthat
via thetimestepislimitedifhigh-energystatesareincluded
inthebasis;forexampleincreasingthebasissetincalcu-
T= (H SI). (12) lations of quantum chemical systems requires a decrease
− −
in the time step.
Note that the minimal eigenvalue of H corresponds to
Rather than considering all possible spawning events,
the maximal eigenvalue of T.
itiscomputationallyefficienttoallowa psipatD to at-
The FCIQMC algorithm proposed by Booth et al.13 i
tempttospawnonlyontoitsowndeterminantandasin-
may be summarized as follows. Consider a collection of
gleconnecteddeterminantD pertimestep. (Connected
j
markersdistributedoverthespaceofdeterminants, D .
i determinants are linked by non-zero transition matrix
{ }
The markers do not move through the Hilbert space, so
elements T .) The connected determinant is chosen
followingAnderson6 wecallthem“psiparticles”orpsips ji
stochastically according to some generation probability
instead of “walkers”. (It is easy to devise alternative
p (j i),whichmustbe non-zerowheneverT isnon-
gen ji
FCIQMCalgorithmsinwhichthe psipsdomoveandbe- ←
zero and normalized such that p (j i) = 1. To
havemuchlikewalkersinDMC,8butweusetheapproach j gen ←
maximisetheprobabilityofsuccePssfulspawning,thegen-
of Booth, Thom and Alavi here.) Each psip has both a
erationprobabilitiesshouldbeassimilarto T / T
location, i, and a “charge”, q = 1, associated with it. | ji| j| ji|
± as possible; in practice, we go for speed and set tPhe gen-
In one time step, ∆τ, we loop over the population of
eration probabilities for spawning on all connected de-
psips and allow each to “spawn” children (new psips)
terminants equal. Once the candidate spawning location
located on new determinants (which may be the same
j has been chosen, the spawning event is accepted with
as the parent’s determinant) according to the following probability T ∆τ/p (j i).13 The effective size of
ji gen
rules: | | ←
the Hilbert space can be reduced by only selecting can-
didate spawning locations of some desired symmetry.
Theprobabilitythatapsipatisuccessfullyspawns
• The population evolution of the psips over the course
a child at j is T ∆τ.
ji
| | of anFCIQMC simulationhas some commonly observed
If the spawning event is successful, the child psip features,asillustratedinFig.1;wedevelopsimplemath-
•
has charge q =sign(T )q , where q is ematical models to understand these results in Sections
child ji parent parent
the charge of the parent. III and IV. The energy shift S is initially set to a value
greater than (i.e., less negative than) the ground-state
At the end of the time step, pairs of psips of opposite
energy eigenvalue, ensuring that the largest eigenvalue
charge on the same determinant cancel each other out of T = SI H is positive and hence that total psip
(“annihilate”)andareremovedfromthe simulation. For −
population grows exponentially. Once the population
example, psips on determinant D with negative diago-
i reaches a critical size the simulation spontaneously en-
nal elements T may generate children which immedi-
ii ters a “plateau” regime, where the total psip population
ately annihilate with the parent psips. Thus, although
remains stable and during which the ground-state wave
psips on different determinants may have different signs,
function emerges. The height of the plateau relative to
all psips on any given determinant have the same sign
the size of the space of Slater determinants provides a
at the end of every time step. Any psips which remain
measureforhowhardasystemistosolveusingFCIQMC.
(whether originally “child” or “parent”) are merged and
After a system- and S-dependent waiting time, the sim-
areallowedtospawnnewpsipsinsubsequenttimesteps.
ulation exits the plateau phase and the psip population
The FCIQMC algorithm13 actually yields a stochastic
starts to grow in an exponential fashion again, albeit at
sampling of the solution of the iterative equation
a slower rate than before. At this point the distribution
ofpsipsis astochasticrepresentationofthe FCI ground-
c (τ +∆τ)= (δ +T ∆τ)c (τ), (13)
i ij ij j state wave function. Allowing the psip population to in-
Xj crease further serves only to reduce the statistical noise,
so it is convenient to hold the population constant by
which may be regarded as a finite-difference approxima-
varying the shift from this point in the simulation. The
tion to the imaginary-time Schr¨odinger equation
changesinshiftmustbecarriedoutsmoothlyandslowly
dc to avoid introducing a bias. Following Booth,13 we use
i
= Tijcj. (14) the shift-update algorithm originally proposed for DMC
dτ
Xj simulations by Umrigar16:
Assuming that S = E0, as is the case once the sim- ξ Np(τ +A∆τ)
S(τ +A∆τ)=S(τ)+ ln , (15)
ulation has converged on the ground state, the domi-
A∆τ (cid:18) N (τ) (cid:19)
p
nant eigenvectors of δ +T ∆τ and T are identical if
ij ij ij
∆τ(E E )<2,whereE isthelargesteigenvalue where A is the number of iterations between which the
max 0 max
−
oftheFCIHamiltonianmatrixH. Thismeansthatthere shift is updated, ξ is a damping parameter and N (τ) is
p
is no time-step error if ∆τ is small enough. In order for the total psip population at time τ. In the simulations
there to be a finite time step, the Hamiltonian must be presentedherewesetA=20andξ =0.1. Oncethe psip
4
populationhassettleddown,theexpectedvalueq¯ ofthe may also be written in reciprocal space,
i
sum of the charges of the psips on D is proportional to
i
theweightofD intheexactground-statewavefunction:
i
U
τl→im∞q¯i(τ)∝c0i. (16) Hˆ =Xk,σ ǫkcˆ†k,σcˆk,σ+M k1X,k2,k3cˆ†k1,↑cˆ†k2,↓cˆk3,↓cˆk1+k2−k3,↑,
(20)
Note that annihilation of psips of opposite charge does
not affect q¯i. where cˆ†k,σ = √1M rcˆ†r,σeik·r is the creation opera-
There are two main ways of obtaining an estimate of tor for a Bloch statPe of wave vector k, the sums are
the ground-state energy from an FCIQMC simulation. over the M wave vectors in the first Brillouin zone,
Theshift estimator isthemeanvalueoftheshiftS inthe ǫk =−2t di=1cos(kia)istheone-electronkineticenergy
constant-population-mode simulation after the plateau. for wave vPector k = (k1,k2,...,kd), and a is the lattice
The projected estimator is obtained by noting that, for parameter(whichissettounityfromnowon). Thecom-
any single determinant D0 with a non-zero ground-state binationofwavevectorsk1+k2 k3 is assumedto have
−
component, the quantity been reduced into the first Brillouin zone by addition of
a reciprocal lattice vector if necessary.
D Hˆ Ψ(τ)
0
E = lim h | | i (17) Whilst the FCIQMC algorithm may be used in the
τ→∞ hD0|Ψ(τ)i real-andreciprocal-spacerepresentations,here,forillus-
H q¯
j 0j j trativepurposes,wefocusonthereciprocal-spaceformu-
= (18)
P q¯ lation, as would be appropriate if U/t were small. As
0
U/t increases and so electrons become increasingly more
tends to the ground-state energy as τ . Since the localised, the real-space formulation becomes more ap-
→ ∞
Hamiltonian operator only links determinants differing propriate. Indeed, as we show in Sec. V, the real-space
bytwoorfewerexcitations,thenumberoftermsincluded basisissubstantiallymorefavourablefornon-zerovalues
inthesumislimited. Notethatitisimportanttoaverage of U/t in the one-dimensional Hubbard model.
q (τ) and q (τ) separately; the projected estimator and
j 0
its associated error can then be found by taking the ra- The imaginary-time evolution of the psip population
tiooftheaveragesandusingthecovariance,respectively. and energy estimators during a k-space FCIQMC sim-
TakingD tobethedeterminantwiththelargestoverlap ulation of the 2D 18-site half-filled Hubbard model at
0
withthe exactground-statewavefunctionminimises the U/t=4(whichisnot small)isshowninFig.1. Asmight
relative stochastic noise in the denominator of the above be expected, the k-space FCIQMC method becomes less
equation. Furthermore,suchadeterminantwilltypically efficientatfindingthegroundstateofthis systemasU/t
have single and double excitations which also have sig- increases — in fact, the psip population at the plateau
nificantcontributions to the ground-statewavefunction, increases approximately linearly with U/t for U/t > 2
and hence determinants contributing to the numerator (Fig. 2). When U/t 4, one needs at least as many
≥
will also often have a significant population. The de- psips as there are determinants in the Hilbert space and
terminant with greatest overlap may not necessarily be the memory requirements are comparable to those of it-
knowna priori (or evenbe clearly defined, as is the case erative diagonalization. As a result, unless U/t is small,
insystemsstudiedinSec.V)butinpracticetheHartree– k-spaceFCIQMCis notabletotreathalf-filledHubbard
Fock determinant is usually a good choice. The choice model systems substantially larger than those accessible
of D can be changed during the course of a simulation to the FCI method.
0
if a determinant with a particularly large population is
found. Working in a space of Slater determinants has two
We can illustrate the main features of the FCIQMC main advantages. As the basis is anti-symmetric with
method by applying it to the Hubbard Hamiltonian17 respectto exchangeoftwospin-orbitals,thereis noneed
for a d-dimensional cubic lattice, to use the fixed-node approximation to prevent collapse
to the bosonic groundstate, as in DMC andGFMC. Al-
though FCIQMC still has a sign problem (see Sec. III),
Hˆ =−thrX,r′i,σcˆ†r,σcˆr′,σ+UXr nˆr,↑nˆr,↓, (19) tahnedinisstoafbteilnityleisssnsoetvewreit.hTrehsepescetcotnodthaedvaabntoasgoeniicsstthaatet
the annihilation of psips with opposite charges proves
where cˆ†r,σ (cˆr,σ) creates (destroys) an electron of spin σ surprisingly efficient in the discrete space of Slater de-
onsite r, the number operatornˆ =cˆ cˆ counts the terminants. Walker cancelationcan in principle cure the
r,σ †r,σ r,σ
spin-σ electronsonsiter,andthe summationover r,r signproblemincontinuumDMCandGFMCsimulations
′
includes allnearest-neigborpairsoflattice sites. Ashsumi- too,21–30butsuchalgorithmsaresubstantiallymorecom-
ing a system of M sites and 2M spin-orbitals subject to plicated and less successful than the sign cancelation in
periodic boundary conditions,the HubbardHamiltonian FCIQMC.
5
coupled differential equations:
8sips/10 122...505 (a) ddnτ+i =Xj (cid:16)Ti+jn+j +Ti−jn−j (cid:17),
p (21)
1.0
No.of −01.45 (b) S(τ) ddnτ−i =Xj (cid:16)Ti+jn−j +Ti−jn+j (cid:17).
−15
E(τ) These can be decoupled by adding and subtracting:
Energy/t −−−−11119876 EFCI d(cid:0)n+id+τ n−i (cid:1) =Xj (cid:16)Ti+j +Ti−j(cid:17)(cid:16)n+j +n−j (cid:17),
−200 5 10 Uτ 15 20 d(cid:0)n+id−τ n−i (cid:1) =Xj (cid:16)Ti+j −Ti−j(cid:17)(cid:16)n+j −n−j (cid:17). (22)
As τ , n+ + n tends to the eigenvector core-
−
spondi→ng t∞o the largest eigenvalue of T+ +T , whilst
FIG. 1. Population dynamics and energy estimators for the −
2D18-sitesquarelatticehalf-filledHubbardmodelatU =4t. n+ n− tends to the eigenvector corresponding to the
Periodic boundary conditions are applied to a 45◦ tilted larg−est eigenvalue of T+ T−. We wish to find the
−
square cell with sides 3√2. The ground-state wave function latter state. However, as explained below, the largest
hasmomentumk=(0,0) andtheHilbertspaceformed from eigenvalue of T+ + T is always larger than that of
−
all determinantsofthismomentumcontains1.3 108 states. T+ T . Thus,thesignaln+ n alwaysdecays relative
Graph(a)showshowthepsippopulationevolves×withsimula- to n−++−n . Fig. 3 shows tha−t p−erforming an FCIQMC
−
tion time. The psip population initially grows exponentially
simulationwithoutannihilationdoesindeedgivethisun-
before reaching a plateau, during which the distribution of
desired state. We note that a similar analysis has been
psips settles to model the FCI wave function. The psip pop-
performedpreviously for the world-line Quantum Monte
ulation then begins to grow exponentially again, reinforcing
Carlo method31,32.
the wave function signal, at which point the shift is allowed
tovarytostabilizethenumberofpsips. Graph(b)showsthe One can provethat the largesteigenvalue of T++T−
shiftandprojectedenergyestimatorsasafunctionofsimula- is always greater than or equal to the largest eigenvalue
tiontime. TheFCIenergyisfromRef.18. FCIQMCrequires ofT+ T− asfollows. Supposethatc0 istheeigenvector
roughly the same number of psips as the size of the Hilbert corresp−onding to the largest eigenvalue, λmax, of T+
diff −
space to converge on the ground-state wave function in this T . If c is normalized and the Hamiltonian matrix is
− 0
case. real,allcomponentsofc mayalsobechosenrealandso
0
λmdiaffx = c0i(Ti+j −Ti−j)c0j. (23)
Xi,j
Nowconsiderthevector c withcomponents c . This
III. THE ORIGIN OF THE SIGN PROBLEM IN FCIQMC | 0| | 0i|
vectorisalsonormalizedandmaybeusedasatrialstate
in the (upside-down) variationalprinciple for the largest
eigenvalue, λmax, of T++T :
Booth et al.13 showed that annihilation is crucial in sum −
FCIQMC; without it the simulation never converges to λmsuamx ≥ |c0i|(Ti+j +Ti−j)|c0j|. (24)
the FCI ground state. In this section we attempt to ex- Xi,j
plain why annihilation is required and how it helps the
The right-hand side of Eq. 24 is manifestly greater than
ground state to emerge.
orequaltotheright-handsideofEq.23,soλmax λmax.
sum ≥ diff
The discrete annihilation process is difficult to model
First,considertheeffectofremovingannihilationfrom
exactly33inadifferentialequation. Wethereforeconsider
the simulationprocedure. This meansthat psips ofboth
a simpler annihilation process in which pairs of psips of
charges are permitted to reside on and propagate from
opposite sign on the same determinant annihilate each
the same determinant at the same time. We shall use
n+i and n−i to represent the populations of positive and ocothnesrtaanttaancdonthsteanfatctroarteof2κ2,iswihnetrreodκucisedassomlealyllfpoorsaitlgivee-
negative psips on determinant D . It is convenient to
i braicconvenience. Thisleadstothedifferentialequations
write the transition matrix T as T+ T , where T+
−
contains the positive transition matrix−elements, T+ = dn+
max(Tij,0), and T− contains the absolute values ofijthe dτi =Xj (cid:16)Ti+jn+j +Ti−jn−j (cid:17)−2κn+i n−i ,
Tnofe+gpaaotnsividteiTvee−leamanredenttnsh,eegrTaei−tfjoivr=eenpmosnaipx-sn(−eegvTaoitjli,vv0ee)..aTcAcholelrdpeilonepmgueltanottisothnoesf ddnτ−i =Xj (cid:16)Ti+jn−j +Ti−jn+j (cid:17)−2κn+i n−i . (25)
6
FIG. 2. Graphs (a)–(d) show the psip population dynamics for the 2D 18-site square lattice half-filled Hubbard model at
k=(0,0) for various values of U/t. Periodic boundary conditions are applied to a 45◦ tilted square cell with sides 3√2. The
dataforU =4t in(b)isidenticaltothatin Fig.1. Graph(e) showsthat thenumberof psipsat theplateau increases linearly
withU/twhenU/t 2. Thestandarderrorsinthenumbersofpsipsduringtheplateauphaseswereobtainedusingablocking
analysis19 andthest≥raightlinewasfittedusingthemethodofleastsquaresasimplementedinRef.20. Unlessshown,theerror
bars are smaller than themarkers.
Eq. 22 thus become
0.00
−−00..1005 d(cid:0)n+id+τ n−i (cid:1) =Xj (cid:16)Ti+j +Ti−j(cid:17)(cid:16)n+j +n−j (cid:17)−4κn+i n−i ,
−−00..2105 Eλ((ττ)) d(cid:0)n+id−τ n−i (cid:1) =Xj (cid:16)Ti+j −Ti−j(cid:17)(cid:16)n+j −n−j (cid:17).
−0.25 E0 (26)
λ0 ItisclearhowannihilationenablestheFCIQMCmethod
0 50 100 150 200 250 to converge upon the ground state of the Hamiltonian
τ/a.u. matrix: asthetotalpsippopulation, i(n+i +n−i ),rises,
the rate of annihilation events rises qPuadratically. This
FIG. 3. The impact of removing annihilation from the causestherateofgrowthofthen++n signaltodecrease
−
FCIQMC algorithm applied to a 100 100 randomly- until spawning of new psips is balanced by annihilation.
generated real symmetric Hamiltonian ma×trix H. E0 is the The growth of n+ n , the desired solution, is not af-
−
exactground-stateeigenvalueandE(τ)theprojectedestima- −
fected and so this signal eventually emerges. As shown
tor for H. λ and λ(τ) are the analogous quantities for the
matrix H++0H−, where H=H+ H− and all elements of inSec.IV,the ratesofgrowthofthe twosignalsactually
H+ andH− arenegative(i.e.,allel−ementsofT+ andT− are becomethe samebutduetoannihilationallpsips onthe
positive). Thesimulation convergestothelowest eigenvector same determinant have the same sign.
of H++H−,which does notcorrespond toan eigenvectorof Tosummarize,thesignprobleminFCIQMCoriginates
theHamiltonian matrix. from the in-phase combination of positive and negative
psips, which grows at a rate determined by the largest
eigenvalue of T+ +T . This eigenvalue is greater than
−
the largesteigenvalueofT+ T ,whichdetermines the
−
−
growth rate of the physical ground state of the system.
The decoupledordinarydifferentialequations(ODEs)in The severity of the sign problem depends upon the dif-
7
ference between the largest eigenvalue of T+ +T and where α=c n(0). Hence the evolution equation for p
− 0
the largest eigenvalue of T+ T — and thus on the becomes ·
−
−
prevalence of negative off-diagonalelements of T — and
dp
upontheconcentrationofpsipsrequiredtoachieveasuf- i = V p κp2+κα2e2Tmaxτc2. (29)
ficient rate ofannihilationfor the growthofthe in-phase dτ Xj ij j− i 0i
signaltobesuppressed. Adifficultsignproblemdoesnot
necessarily imply a strongly correlated system. Eq. 29 is difficult to solve exactly, but the population
Eq.21models the time dependence ofthe populations dynamics it embodies can be illustrated using a simple
ofpositiveandnegativepsipsinasimulationinwhichall one-component analogue:
annihilation is forbidden. In a realFCIQMC simulation,
however, the annihilation rate per psip does not tend to dp =V p κp2+κn2, (30)
max
zero as the psip density tends to zero. The probability dτ −
of encountering and annihilating an unrelated psip van-
where p(τ) is the total psip population at time τ, V
ishes,butaparentmaystillspawnachildoftheopposite max
is larger than T , and n(τ) = n eTmaxτ. It is common
signonitsowndeterminant,inwhichcasetheparentand max 0
to start an FCIQMC simulation with a population of
child annihilate each other. The low-density limit of the
positive psips only, in which case the initial condition is
standardFCIQMCalgorithmisbettermodeledbydefin-
p(0)=n(0)=n . As the initial psip population and an-
ing T+ and T in a different way: T+ now contains all 0
− nihilationratearesmall, p(τ) growsexponentiallyatthe
diagonal elements of T, regardless of sign, and all pos-
start of the simulation: p(τ) = n eVmaxτ. The exponen-
itive off-diagonal matrix elements; whilst T is zero on 0
− tial growth ceases when p V /κ, at which point the
thediagonalbutcontainstheabsolutevaluesofthenega- κp2 annihilationterm(whi≈chismlaaxrgerthantheκn2 term
tive off-diagonal elements. This redefinition corresponds
because V > T ) counteracts it. The psip popula-
to a change in viewpoint: instead of allowing parents to max max
tion then enters a plateau phase, remaining rougly con-
spawnchildrenoftheoppositesignontheirowndetermi-
stantuntilthen(τ)signal,whichhasbeensteadilygrow-
nant,withsubsequentannihilation,thenegativediagonal
ing like n eTmaxτ, begins to dominate. From then on the
elements of T+ remove psips from the simulation in one 0
population rises exponentially again, although now at a
step, introducing an exponential decay of the psip pop-
ulation on determinants for which T+ < 0. The psip rate determined by Tmax. The ground-state wave func-
ii tion thus spontaneously emerges from the simulation.
densities n+i and n−i remain positive at all times and The plateau begins to appear at a time τ determined
the above analysis is unchanged, but the severity of the 1
by the equation p(τ ) n eVmaxτ1 V /κ. Hence
sign problem,as measuredby the difference between the 1 ≈ 0 ≈ max
largesteigenvalues ofT++T and T+ T , is smaller
− − ln(V /κn )
− max 0
than the above analysis suggests. τ1 . (31)
≈ V
From now on it will be assumed that T+ and T are max
−
definedasexplainedinthepreviousparagraph,andhence Moreprecisely,solvingEq.30withtheκn2 termomitted
that all diagonal elements of T are zero.
− and assuming that κn /V 1 (which must be the
0 max
≪
caseiftheannihilationrateisnegligibleatthebeginning
of the simulation) shows that p(τ) reaches 95% of its
IV. POPULATION DYNAMICS plateau value at time
Separating the psip population into positive and neg- ln(19(Vmax/κn0 1))
τ = − . (32)
95%
ative contributions also explains the population dynam- Vmax
ics observed during an FCIQMC simulation. Let V =
T++T , p = n++n and n = n+ n . The decou- The end of the plateau occurs at a time τ2 determined
− − − − by the equation n(τ )=n eTmaxτ2 V /κ, and hence
pled ODEs in Eq. 26 can thus be written as 2 0 ≈ max
ln(V /κn )
dp max 0
i = V p κ(p2 n2) τ2 . (33)
dτ ij j− i − i ≈ Tmax
Xj
(27)
dn Combining Eqs. 31 and 33 yields
i
= T n .
ij j
dτ
Xj τ2 Vmax
. (34)
τ ≈ T
1 max
In the long-time limit, n(τ) becomes dominated by its
ground-state component, which grows at a rate deter- The energy shift S appears on the diagonal of T =
mined by the largest eigenvalue, T = S E , of the SI Handis incorporatedintothe diagonalelementsof
max 0
transition matrix T: − T+−;itthereforeaffectsT=T+ T− andV=T++T−
−
equally. Thisdependencecanbemadeexplicitbywriting
n(τ) αeTmaxτc , (28) T =S+T andV =S+V ,whereT andV ( T )
0 max 0 max 0 0 0 0
≈ ≥
8
Vmax/κn0=25.0;Vmax/Tmax=30.0
Vmax/κn0=12.5;Vmax/Tmax=15.0
Vmax/κn0=12.5;Vmax/Tmax=30.0
100
80
0 60
n
/
)
τ
p( 40
20
FIG. 4. Effect of the initial value of the shift on the plateau
height in an FCIQMC simulation of the 2D 18-site square
lattice half-filled Hubbard model at U = 2t and k = (0,0). 0
Periodic boundary conditions are applied to a 45◦ tilted cell 0 20 40 60 80 100
with sides 3√2. The reference determinant, D0, is formed Vmaxτ
fromoccupyingthe18lowest-energyBlochfunctions. Setting
theshiftslightlybelowH00=hD0|Hˆ|D0ireducestheplateau FIG. 5. Model population dynamics in an FCIQMC simula-
height butincreases the time spent in the plateau region.
tion. The psip population evolves according to Eq. 30 and
shows the key features of the population dynamics in real
FCIQMC calculations: an initial exponential growth phase
arethelargesteigenvaluesofTandV whenS =0. This
followedbyaplateaufollowedbyasecondslowerexponential
leads to the equation growth phase. Asexpected,thepsip population reaches95%
of its plateau value when V τ ln(19(V /κn 1)),
τ S+V S+V max 95% ≈ max 0−
2 0 = 0. (35) the height of the plateau is given by p/n0 = Vmax/κn0, and
τ1 ≈ S+T0 S−E0 the plateau ends when Vmaxτ2 ≈ ln(Vmax/κn0)Vmax/Tmax.
The vertical dashed (dotted) lines show V τ (V τ )
max 95% max 2
In most FCIQMC simulations the shift is initially held
foreachsetofparameters. Thehypergeometricfunction, F ,
0 1
at a fixed value, typically the Hartree-Fock energy, until was calculated using Ref. 20.
the desired psip population is reached. If S is reduced
and allowed to approach E from above, the simulation
0
never exits the plateau. The duration of the plateau can that V /(2T ) is not an integer. The normalization
max max
be shortened by moving the initial value of S further of u(τ) drops out of p(τ) = (κu) 1du/dτ, leaving one
−
abovetheground-stateenergyE0,butonlyatthecostof arbitrary constant to be fixed by the initial conditions.
increasingthe psippopulation(S+V0)/κattheplateau. Fig.5 shows the resulting populationdynamics for three
The effect of modifying the initial shift in an FCIQMC differentsets ofparameters,imposingthe boundarycon-
simulation is shown in Fig. 4. dition p(0)=n in all cases. Whilst this is an extremely
0
TheODEinEq.30isanexampleofaRiccatidifferen- simple model, it captures allof the features of the popu-
tial equation34 and canbe transformedfroma quadratic lation dynamics seen in actual FCIQMC calculations.
first-orderODEintoalinearsecond-orderODEusingthe LetusnowreturntoEq.29andconsiderwhathappens
substitution beyondtheendoftheplateau,wherethepsippopulation
beginstoriseagainuntilthe energyshiftisadjustedand
1 du
p(τ)= . (36) a steady state with dp /dt=0 is attained. Since the lin-
κudτ i
earVtermisnegligibleincomparisonwiththequadratic
The solution of the resultant second-order ODE is terms in this regime, Eq. 29 becomes
u(τ)=c1 0F1 ;1 Vmax ;z 0= dpi κp2+κn2, (39)
· (cid:18) − 2Tmax (cid:19) dt ≈− i i
(37)
V
+c2zVmax/2Tmax ·0F1(cid:18);1+ 2Tmmaaxx;z(cid:19), i0m. pDlyeintegrmthiantanptis≈on±nwihaicnhdthheencgerotuhnadt-nst−iat≈e 0amorplnit+iud≈e
c is positive are occupied only by positive psips and
where 0i
determinants on which c is negative are occupied only
0i
κ2p2(0)e2Tmaxτ by negativepsips. The signedpsipdensity ni =n+i −n−i
z = , (38) isproportionaltoc andtheunsigneddensityp =n++
4T2 0i i i
max
n−i to |c0i|.
c and c are constants of integration, F is a confluent We canalsouseEq.29to understandwhy the plateau
1 2 0 1
hypergeometric limit function,35 and we have assumed psippopulationsplottedinFig.2areproportionaltothe
9
Hubbard U. Assuming that the growth of the in-phase inthisexamplearespinconfigurationsratherthanSlater
signal, p, is fast enough to allow the plateau to emerge determinants.) Indeed, the sign structure of the ground-
before the ground-state signal becomes significant, the statewavefunctionofthebipartiteHeisenbergmodelhas
plateau occurs when long been known.36 Every non-zero off-diagonal matrix
element of the Heisenberg Hamiltonian flips a neighbor-
V p κ p2. (40)
ij j ≈ i ingpair ofspinsfromdown-upto up-downor vice-versa.
Xij Xi The Hamiltonian matrix element is positive if J > 0,
If U/t is large, the kinetic energy contributions to the implying that the corresponding matrix element of T is
Hamiltonian matrix (and hence T and V) can be ne- negative,sothesignsofchildrenproducedbyoff-diagonal
glected and changing U simply scales these matrices. spawningeventsalwaysdifferfromthesignsoftheirpar-
Defining V = UV and p = Up leads to the follow- ents. Initially it appears as if this ought to cause a seri-
ij i′j i ′i
ing condition for the emergence of the plateau: ous sign problem. Consider, however,a bipartite system
consistingoftwointer-penetratingsub-lattices,arranged
2
Vi′jp′j ≈κ p′i , (41) such that the pairs of spins flipped by off-diagonal ele-
Xij Xi ments of Hˆ are always on different sub-lattices. The ac-
where V and p contain no dependence upon U. The tion of any off-diagonal element of Hˆ increases the total
′ ′
valueofS ononesub-latticeby1anddecreasesthetotal
total psip population at the plateau, p = U p , z
i i i ′i
value of S on the other sub-lattice by 1. This allows all
therefore scales linearly with U. P P z
spin configurations to be assigned to one of two classes,
A or B. Starting from the classical N´eel state, which is
V. SIGN-PROBLEM-FREE SYSTEMS arbitrarilyassigned to class A, all states that differ from
the N´eel state by an even number of applications of Hˆ
are also assigned to class A, while states that differ by
If there exists a similarity transformthat maps T into
anoddnumber ofapplicationsofHˆ areassignedtoclass
V,andhencemakeseveryoff-diagonalelementofTpos-
B. If an FCIQMC simulation is initialized by placing a
itive and every off-diagonalelement of H=SI T neg-
− population of positive psips on the N´eel state, all psips
ative, then V and T have identical eigenvalues. Even
createdonconfigurationsinclassAwillbepositivewhile
without annihilation, the in-phase and out-of-phase sig-
allpsips createdonconfigurationsinclass B willbe neg-
nals, p and n, grow at the same rate and there is no
ative. Psips of different signs will never mix, no cance-
difficulty extracting the ground-state (out-of-phase) sig-
naln=n+ n . Suchsystemsaresign-problemfreeand lation is required, and there is no sign problem. In fact,
−
− byapplyingasimpleunitarytransformationinwhichthe
FCIQMC simulations of them yield correctground-state
sign of every configuration in class B is flipped, all off-
energieswitharbitrarilysmallpsippopulations,although
diagonal matrix elements of H can be made negative.37
no plateau phase is seen. Increasing the psip population
servesonlytoreducethestochasticerror. Indeed,byfor- By way of an example, the Heisenberg model for a
mulating the problem in the transformed basis, we can two-dimensional6 6squarelatticewithperiodicbound-
carry out an FCIQMC simulation in which no negative ary conditions ha×s a Hilbert space of 9.08 109 con-
psips need ever appear. figurations. Using just 1.8 106 psips w×e find the
×
Aparticularlysimpletypeofsimilaritytransformation E/N = 0.67886(2)J, which is in excellent agreement
−
does no more than change the signs of some of the basis withthevalueE/N = 0.678871(8)JobtainedbyRunge
determinants. No choice of signs is sufficient to render using Green’s function−Monte Carlo.38 Note that, unlike
alloff-diagonalelementsofT positiveinmostcases,but Runge, we did not use importance sampling. FCIQMC
ij
there are a few interesting model systems in which sim- simulations of Heisenberg models defined on other (non-
ple sign-changing transformations are effective. In some bipartite or “frustrated”) lattices display the same be-
models, for example, all psips spawned on any given de- havior as for fermionic systems (Fig. 6); in particular,
terminant D have the same sign regardless of the loca- the population plateau is a universal feature unless the
i
tion of their parent, so that positive and negative psips system is sign-problem free.
never mix. A simple sign-changing transformation that The Hubbard Hamiltonian, Eq. 19, is closely related
multiplieseverybasisdeterminantbythesignofthepsips to the Heisenberg Hamiltonian and is also sign-problem-
that occupy it then makes all off-diagonal elements of T free, although only in one dimension. Consider a 1D
positive. HubbardlatticewithN sitesandN (N )spin-up(spin-
s
The antiferromagnetic Heisenberg model defined by down)electrons,wherethesystemis↑not↓necessarilyhalf-
the Hamiltonian filledandperiodicboundaryconditionsareapplied. IfN
s
Hˆ =J Sˆr ·Sˆr′, (42) Nis evaenndaNndaNre↑baonthdeNv↓ena,rtehebnotthheordedexoirstNsasnisanoadldagaonuds
hXr,r′i tra↑nsform↓ation to that described above for the Heisen-
where J>0, Sˆ is the vector spin operator on lattice site berg model. (This is most easily seen if the orbitals are
r
r,andthesumrunsovernearestneighbors,issuchasys- ordered first by spin and then by their position in the
tem ifthe lattice is bipartite. (Note that the basis states lattice.) The ground-state energies of large 1D Hubbard
10
whereasthatwithBlochorbitalscontainsonly1.31 108
×
determinants.
VI. DISCUSSION
The aboveobservationsprovidesomedegreeofinsight
intotheFCIQMCmethodandthefactorsthatdetermine
how hard it is to apply FCIQMC to any given physical
system. We now briefly comment upon other topics re-
lated to FCIQMC before summarizing our work.
TheinitiatorapproximationtoFCIQMC(i-FCIQMC),
FIG.6. Thepopulation dynamicsofan FCIQMCsimulation wherebyspawningeventsontopreviouslyunoccupiedde-
of the antiferromagnetic Heisenberg model for a 5 5 trian- terminantsareonlyallowediftheparentdeterminanthas
×
gularlattice. Noperiodicboundaryconditionswereimposed. a psip population above a specified threshold, has been
Aswithfermionicsystems,frustratedHeisenbergmodelsdis- shown to dramatically reduce the number of psips re-
play a population plateau. quired to obtain excellent results in many systems.15 In
i-FCIQMC,theeffectiveHilbertspacegrowsandchanges
dynamically as the simulation proceeds, including only
models with many different choices of N and N can determinants that have a significant psip population or
thereforebe foundusingtheFCIQMCmet↑hodwith↓very lie close (in terms of applications of the Hamiltonian) to
small numbers of psips. determinantswithasignificantpsippopulation. Asare-
Thisbringsusontoanimportantpoint: thesignprob- sult,psips arepreventedfromspawninginregionsoflow
lem in FCIQMC is not constant for a given system but psipdensityandtheannihilationrateisgreatlyincreased
is dependent upon the choice of basis. Consider two for a given psip population. No plateau is observed in
Hamiltonian matrices, H and H , which describe the i-FCIQMC calculations14,15 because the growth of the
1 2
same system but have been constructed using different in-phase combination is suppressed by annihilation even
basis sets, one of which is obtained from the other by a when the total psip population is low. The i-FCIQMC
transformation matrix S, so that H = S 1H S. The approximationbecomesincreasinglygoodasthenumber
1 − 2
correspondingtransitionmatrices T andT arerelated of psips is increased because more of the Hilbert space
1 2
in the same way. However, when we partition the two becomes accessible and the psip dynamics more closely
transition matrices into T+ and T , the partitions are resembles that of the true Hamiltonian.
−
not necessarily related by the same unitary transforma- IncoupledclusterMonteCarlo39 (CCMC),theexcita-
tion. Inotherwords,althoughT1 =S−1T2S, andhence tion amplitudes of the coupled cluster wave function are
tTh+1at−TT±1−1 ==SS−−11T(T±2+2S−. TC−2on)Sse,qiuteinstnlyo,ttihnegepnlaetreaaluthheecigahset sptloincghaosftcicoanlfilygusaramtpiolnedaminpmlitaundneesriannFaCloIgQouMsCt.oAthltehosaumgh-
(and thus the ease with which FCIQMC can be used to theCCMCalgorithmismuchmorecomplicatedthanthe
find the groundstate) depends on the basis in which the FCIQMC algorithm, we believe that a similar analysis
Hamiltonian matrix is expressed. holds. The particles or “excips”40 that sample the dis-
An FCIQMC simulation of the 18-site half-filled 1D crete excitor space in CCMC also have positive or nega-
Hubbard model at U = t exhibits no plateau when the tive signs and annihilationis crucialfor the desired(and
Hamiltonian matrix is expressed in a basis of determi- physicallymeaningful)solutionto emerge. Itis observed
nants of local orbitals (as in Eq. 19), since this partic- thatCCMCcalculationshaveahigherplateauthancon-
ular problem is sign-problem free as explained above. figuration interaction quantum Monte Carlo (CIQMC)
Yet an FCIQMC simulation of exactly the same Hamil- calculations at the same truncation level.39 It is likely
tonian expressed in a basis of determinants of Bloch that this is because the excips have to explore a larger
functions (as in Eq. 20) exhibits a plateau at 6.9 106 effective Hilbert space than the psips in the equivalent
×
psips. As expected, the sign problem increases the dif- CIQMC simulation, so more excips are requiredto make
ficulty of the FCIQMC simulation in the Bloch repre- the annihilation rate high enough to suppress the in-
sentation. A short simulation using just 2.8 105 psips phase signal.
×
in the local orbital basis gave a ground-state energy of In summary, we have shown that the sign problem in
18.8423(3)t,whereas2.3 107psipswererequiredtoob- FCIQMC simulations is caused by the growth of an un-
− ×
tainaground-stateenergyof 18.84248(8)tintheBloch physical dominant solution of the coupled ODE’s that
−
orbitalbasis. Furthermore,ourimplementationcurrently describe the time evolution of the densities of positive
conserves crystal momentum when using Bloch orbitals andnegativepsips. Thisunphysicalsolutiongrowsfaster
but does not make use of symmetry when using local than the ground-state solution we seek and masks the
orbitals; the Hilbert space used in the calculation with ground-state signal. A similar problem arises in DMC
localorbitalsthereforecontains2.36 109 determinants, simulations of fermionic systems, where the dominant
×