The sign problem in density matrix quantum Monte Carlo
Hayley R. Petras, William Z. Van Benschoten, Sai Kumar Ramadugu, James J. Shepherd
TThe sign problem in density matrix quantum Monte Carlo
Hayley R. Petras ( i ) , ( ii ) , a) William Z. Van Benschoten ( i ) , ( ii ) , a) Sai Kumar Ramadugu ( i ) , ( ii ) , andJames J. Shepherd ( i ) , ( ii ) b) ( i ) Department of Chemistry, University of Iowa ( ii ) University of Iowa Informatics Initiative, University of Iowa (Dated: February 2, 2021)
Density matrix quantum Monte Carlo (DMQMC) is a recently-developed method for stochastically samplingthe N -particle thermal density matrix to obtain exact-on-average energies for model and ab initio systems. Wereport a systemic numerical study the sign problem in DMQMC based on simulations of atomic and molecularsystems. In DMQMC, the density matrix is written in an outer product basis of Slater determinants and has asize of space which is the square of the number of Slater determinants. In principle this means DMQMC needsto sample a space which scales in the system size, N , as O [(exp( N )) ]. We find that, while DMQMC suffersa sign problem that scales with the size of the outer product space, the interaction picture modification ofDMQMC (IP-DMQMC) has a sign problem which only scales with the size of original determinant space. Weachieve this by establishing that the critical walker population to converge a simulation is the same betweenIP-DMQMC and full configuration interaction quantum Monte Carlo (FCIQMC). This equivalence betweenIP-DMQMC and FCIQMC seems to extend to the initiator approximation, which is often required to studylarger basis sets and other systems. This suggests IP-DMQMC ameliorates the cost of moving between aSlater determinant space and an outer product basis. I. INTRODUCTION
In a recent study, we showed that the density matrixquantum Monte Carlo (DMQMC) method could be ap-plied to molecular systems, extending it beyond origi-nal applications to model systems in condensed matterphysics. The use of finite temperature electronic struc-ture methods are becoming increasingly important in ap-plications such as plasmonic catalysis, the study ofplanetary interiors, and solid-state materials , wherethe temperature dependence is key in obtaining phys-ical and chemical properties, such as phase diagramsand excitation energies. The inclusion of temperaturein quantum chemistry methods is difficult because at fi-nite temperatures, more than one state is often occupied,increasing the difficulty of solving the Schrodinger equa-tion. DMQMC joins a growing set of methods includ-ing other quantum Monte Carlo methods, many bodytheories and others in attempting to solve the finitetemperature problem that has attracted recent attentionamongst quantum chemists. Many of these methods,like DMQMC, continue to undergo development.
Widespread adoption of all methods in the FCIQMCfamily, including DMQMC, is hindered, in part, dueto the sign problem. In FCIQMC-based methods, co-efficients in the wavefunction (or density matrix, inDMQMC) are sampled by a distribution of walkers. Theoriginal FCIQMC paper found that simulations that ex-ceeded a critical walker population were able to success-fully resolve the signs of the wavefunction and generatean energy estimate that was exact-on-average; lower pop-ulations were instead biased. In FCIQMC, walkers ar- a) These authors contributed equally to this paper b) Electronic mail: [email protected] riving at the same site can be exactly annihilated dueto having the discrete basis set; this contrasts a contin-uous real-space basis where the same approach can bemuch more difficult. If walkers never meet because thereare too few, then annihilation cannot occur, leading to abias in the energies. A simulation with a growing walkerpopulation will have its growth briefly stall out, form-ing a plateau in the total walker population ( N w ) as afunction of the simulation iteration, known as the “anni-hilation plateau”, as the simulation establishes the signof critical elements of the wavefunction. When the pop-ulation has grown above the plateau, the sign problem isresolved and exact energies can be collected.The sign problem in FCIQMC was discussed in depthin the early developmental papers for the methodbefore subsequently being systematically studied bySpencer et al. , whose work we refer to throughout. Thiswork established the origin of the sign problem as anunphysical Hamiltonian (whose solution does not have asign problem) and that is unavoidably encountered in un-dersampled dynamics. There are also attempts to lever-age this understanding directly using a fixed-node or trialwavefunction approach. The development of the ini-tiator approach in FCIQMC removed the annihilationplateau at a cost of introducing a small error in the en-ergy (removed by increasing the number of walkers in thesimulation). The motivation for and derivation of thisapproximation was related to the alleviation of the signproblem and allowed for a much broader scope for appli-cations. The development of the initiator approximationin DMQMC achieved a similar outcome allowing for ourprevious work on the uniform electron gas and ab ini-tio molecular systems.
Subsequently, there were alsoa wide variety of FCIQMC or FCIQMC-like methods de-velopments which are beyond the scope of this work toreview in detail.
Large scale implementations of the a r X i v : . [ phy s i c s . c h e m - ph ] J a n FCIQMC method and related methods have also been de-veloped and these papers review current challenges anddevelopments for the interested reader.
Here, we conduct a systematic investigation of thesign problem in density matrix quantum Monte Carlo(DMQMC). We find that the annihilation plateau comesfrom the same unphysical Hamiltonian as in FCIQMC.We measure these critical walker populations for a testset from the FCIQMC literature and find that DMQMCplateau heights are proportional to the square of theFCIQMC plateau height. However, we also show thatmoving to the interaction picture (IP-DMQMC) causesthe DMQMC plateau heights to scale linearly with theFCIQMC plateau heights. Despite being able to controlthe sign problem, IP-DMQMC can have an issue with thetrace population as there is no global estimator for theenergy in DMQMC unlike in FCIQMC. To address thecollapse of the trace population that occurs even whenthe sign problem is overcome, we examine importancesampling and the initiator adaptation. We show that im-portance sampling has a more severe sign problem, whilethe initiator adaptation has similar performance as theinitiator adaptation in ground-state calculations usingFCIQMC. We believe this shows DMQMC is as effectiveat solving for finite-temperature energies as FCIQMC isat solving zero-temperature energies. We see this work ascomplementary to our previous and future studies whichdevelop and apply DMQMC as well as the related work ofRubenstein et al. discussing the sign problem for finite-temperature auxiliary field quantum Monte Carlo.
II. METHODS
In this section, we provide a summary of the meth-ods used here. We begin with the three methods pri-marily used in this work: DMQMC, interaction pictureDMQMC and FCIQMC. We then describe the initiatoradaptation and importance sampling. We note now thatHartree atomic units are used throughout.
A. Density matrix quantum Monte Carlo
We begin with original formulation of DMQMC. Starting with the unnormalized thermal density matrixˆ ρ = e − β ˆ H (1)where ˆ H is the Hamiltonian operator and β = ( k B T ) − ,we can show that the density matrix satisfies the sym-meterized Bloch equation d ˆ ρdβ = −
12 ( ˆ H ˆ ρ + ˆ ρ ˆ H ) . (2)by differentiating ˆ ρ ( β ) with respect to β . A Euler updatescheme, or finite difference approach, with a finite time step, ∆ β , can then be used to find the density matrix atany β , followingˆ ρ ( β +∆ β ) = ˆ ρ ( β ) − ∆ β H ˆ ρ ( β )+ ˆ ρ ( β ) ˆ H )+ O (∆ β ) . (3)We then rewrite Eq. (3) in a basis of outer products ofSlater determinants to obtain a matrix form that can besolved stochastically, by evolving a population of parti-cles through the inverse temperature regime. The resultis ρ ij ( β + ∆ β ) = ρ ij ( β ) + ∆ β (cid:88) k ( T ik ρ kj + ρ ik T kj ) (4)where T ij = − ( H ij − Sδ ij ) is the update matrix, and S isa variable shift for population control of the particles inthe simulation, explained later in this section.The matrix elements ρ ij = (cid:104) D i | ˆ ρ | D j (cid:105) are representedby particles in the simulation, where | D i (cid:105) are Slater de-terminants in the defined finite basis set. A populationof particles, N w is then used to sample elements of thedensity matrix by evolving with respect to β , accord-ing to Eq. (4). Integer weights were used in the originalFCIQMC algorithm (which DMQMC is based off of) andbecause we want to make comparison with previous re-sults, this is what we use here.The simulation starts at β = 0, where the density ma-trix is the identity matrix. The simulation is then propa-gated to the desired value of β . At each step, the popula-tion is updated following rules for spawning and death ofparticles, summarized below, while particles of oppositesigns on each matrix element are annihilated; these stepsare closely analogous to FCIQMC. There are three rules for evolving particles that candescribed as follows: • Spawning: occurs from one matrix element ( ρ ik ) toanother ( ρ ij ), along both the rows and columns. • Cloning and death: occur on single matrix elementsonly, and are designed to increase and decrease thepopulation respectively. • Annihilation: particles of opposite signs on singlematrix elements are removed from the simulation.Spawning will occur occur with the probability p s ( ik → ij ) = ∆ β | T kj | and the sign will correspond tosign( ρ ij ) = sign( ρ ik ) × sign( T kj ). The same equations willhold for spawning from ρ kj to ρ ij . The sign of the newlyspawned particle is important because as can be seen,the sign of the new particle depends on both the signof the matrix element where the particle spawned from and the update matrix, T , connecting the two elements.Because of this, the signs of newly spawned particles willnot always be sign-coherent, resulting in the manifesta-tion of the sign problem. In order to resolve the signs ofthe particles on the matrix elements, a system dependentnumber of particles is required. This will be explored fur-ther throughout this work.Cloning and death occur with a probability given by p d ( ij ) = ∆ β | T ii + T jj | . The population increases ifsign( T ii + T jj ) × sign( ρ ij ) >
0, and decreases otherwise.Annihilation also occurs on single matrix elements andis used to control the sign problem and particle growthwithin the simulation, and has been show to be key inovercoming the sign problem.
A population control must be used, so we introduce avariable shift parameter, that is controlled by S ( β + A ∆ β ) = S ( β ) − ξA ∆ β ln (cid:16) N w ( β + A ∆ β ) N w ( β ) (cid:17) (5)The shift update is dependent on N w ( β ), the total num-ber of walkers at β , A , the number of ∆ β steps betweenshift updates, and ξ , a shift damping parameter.The steps outlined above are repeated until the de-sired inverse temperature is reached. To obtain esti-mates of thermodynamic quantities, one averages overmany independent simulations, termed “ β loops”. Then,to find the energies the following expression is used, (cid:104) ˆ H (cid:105) = T r (ˆ ρ ˆ H ) /T r (ˆ ρ ), where the numerator and denom-inator of this equation are sampled separately over thecourse of the propagation through β , and averaged overthe desired amount of β loops. B. Interaction Picture DMQMC (IP-DMQMC)
The interaction picture variant of DMQMC (IP-DMQMC throughout) was developed to overcome twosampling issues present in the original DMQMC method:the initial density matrix rarely contains the importantdeterminants and the distribution of weight fluctuatesrapidly as a function of β . Replacing the density ma-trix with an auxiliary matrix, ˆ f , means that the simula-tion can be started at a non-interacting density matrix, e − β ˆ H , rather than the identity, providing a good firstapproximation to the fully interaction density matrix forweakly-correlated systems. The auxiliary matrix can bewritten: ˆ f ( τ ) = e − ( β − τ ) ˆ H e − τ ˆ H (6)where ˆ H = ˆ H + ˆ V and ˆ H is a mean-field Hamiltonian.In this work, we use the Hartree-Fock Hamiltonian forˆ H , though it is possible to use a more general mean-fieldHamiltonian. In practice, ˆ H only has diagonal matrixelements in a Slater determinant basis, and e − ( β − τ ) ˆ H only has diagonal matrix elements at any temperature.It is important to note that this matrix evolves from e − β ˆ H at τ = 0 to e − β ˆ H = ˆ ρ ( β ) at τ = β , which meansIP-DMQMC only samples the correct distribution at τ = β , so separate simulations are required for each β value.We can differentiate the matrix ˆ f with respect to τ tofind: d ˆ fdτ = ˆ H ˆ f − ˆ f ˆ H. (7) This equation can be simulated using the rules above,with one change: the cloning/death probability in thesecond rule changes to p d ( ij ) = ∆ τ | H ii − H jj | , becauseˆ H is diagonal in the chosen basis. The condition on in-creasing or decreasing the population is then based on thesign of ( H ii − H jj ) × sign( ρ ij ), with population increasingif this expression is greater than 0.IP-DMQMC is the same as DMQMC, in that manysimulations need to be averaged to obtain estimates forobservables. When introduced, it was said that one ma-jor benefit of this variant is that as long as H ii > H ii ,there is little to no death along the diagonal; this over-comes one problem with large systems in DMQMC,where the distribution along the diagonal approacheszero with β . When H is based on Hartree–Fock, H ii = H ii , the initial condition must also be changedand this is described in detail in the original paper. The grand canonical density matrix corresponding to ˆ H is used to obtain the desired distribution according to e − β ˆ H . C. Full configuration interaction quantum Monte Carlo
Next, we describe the FCIQMC method in brief asit will be used for comparison throughout this work.FCIQMC begins with the imaginary time Schrodingerequation d | Ψ (cid:105) dτ = − ˆ H | Ψ (cid:105) (8)where | Ψ (cid:105) is the ground state wavefunction, ˆ H is theHamiltonian operator and τ represents imaginary time.Here, the wavefunction is represented as a sum overSlater determinants, | D i (cid:105) , | Ψ (cid:105) = (cid:88) i c i | D i (cid:105) (9)where c i is the coefficient on the i th determinant and theHamiltonian is represented as H ij = (cid:104) D i | ˆ H | D j (cid:105) . (10)In the same vein as DMQMC, we can obtain a finitedifference equation c m +1 i − c mi = τ ( − H ii + S ) c mi − (cid:88) j (cid:54) = i τ H ij c mj (11)by substituting the sum over Slater determinants (fromEq. (9)) into Eq. (8), where c mi is the coefficient of the i th determinant at iteration m of the simulation. Notehere that the total population of particles, N w , is givenby N w = (cid:80) i | c i | . To obtain an estimate of the groundstate energy, S is varied to keep the particle populationconstant, and can be averaged to obtain the estimate.The rules for evolving particles are those on which theDMQMC algorithm was subsequently based. At eachstep of the simulation, the particles on each element willundergo spawning, death/cloning and annihilation, asthey do in DMQMC. Particles spawn from site i withweight c i to connected sites, j , where i (cid:54) = j , where theprobability is uniform in j . In the death/cloning step,particles on site i increase or decrease their populationaccording to | S − H ii | τ . Particles on site i with oppositesigns are removed from the simulation.The particle population is evolved through imaginarytime following the rules above, through a system de-pendent number of iterations. After the wavefunctionemerges, the correlation energy is found by averagingover the iterations in the simulation, in a similar fash-ion to how β loops are averaged in DMQMC. D. Initiator Adaptation
The sparsity of the thermal density matrix (or thewavefunction coefficient matrix in FCIQMC) can be uti-lized through use of the initiator approximation varia-tion of both methods, here represented as i-DMQMC and i-FCIQMC . The initiator approximation worksby setting a threshold “ n add ” value, where spawningto unoccupied matrix elements is limited to occur onlyfrom matrix elements with particle populations largerthan n add , called “initiator determinants” (or from co-incident spawns of particles of the same sign from twonon-initiator sites). This approximation limits the num-ber of density matrix elements (or vector elements inFCIQMC) that need to be sampled over the course of thesimulation. Increasing the total number of particles, N w ,can reduce the magnitude of the approximation. Both ofthe original algorithms are obtained as N w → ∞ . Theinitiator adaptation can be used with or without the in-teraction picture. E. Importance Sampling
Finally, we give a brief description of importance sam-pling in DMQMC (represented as IS DMQMC through-out); additional details can be found in Ref.54. Thegoal of importance sampling in DMQMC is to reduce theprobability of particles spawning far from the diagonal .Concentrating sampling on the diagonal matrix elementshelps with the convergence of stochastic error as theenergy expression is focused on these elements. Den-sity matrix weights stored on excitations which are morethan one Hamiltonian action away from the diagonal (i.e.more distant than two-particle excitations) are less likelyto contribute back to the energy directly. The approach isto give particles on higher excitation levels larger weights,in order to avoid changes in the expectation values of thedesired operator. More particles with lower weights onor near the diagonal will help to decrease stochastic er- ror. The number of pairs of opposite signs that mustbe flipped in order to reach | D i (cid:105) from | D j (cid:105) is defined asthe excitation level, as first described by Ref. 54. Im-portance sampling can be used with both DMQMC andIP-DMQMC. F. Kernel Density Estimation
The plateau height in this work is defined as the popu-lation that occurs with the highest frequency in the sim-ulation, and we call this population the critical popula-tion ( N c ). The Scott kernel density estimation (KDE)method is used in this work to assign critical popu-lations through a systematic and reproducible protocol.This is a continuous adaptation from prior work. TheKDE method works by calculating the probability thata certain walker population is present in the DMQMCsimulation through the use of a KDE kernel, K . If welet f ( x ) be a continuous function representing the popu-lation dynamics in one trajectory, we can use the kerneldensity estimatorˆ f h ( x ) = 1 nh K (cid:0) x − x i h (cid:1) (12)where h is a smoothing parameter and n is the numberof data points to find the KDE kernel. The KDE ker-nel itself gives a probability distribution of the numberof walkers as a function of the number of walkers. Themaximum value of the kernel will correspond to the crit-ical walker population.Simulations to measure the plateau height are per-formed with a single β loop, and the output files are an-alyzed using the Python scripts provided in the HANDEsoftware package, producing one analysis file per out-put file. The plateau assignments are performed on thedata sets with the total walker population ( N w ) on a log-arithmic axis. For the plateau heights in DMQMC, therewere some cases where the simulations entered variableshift before the annihilation plateau occurred, or the to-tal population collapsed to zero and did not recover. Ifeither of these situations occurred in the simulation, itwas not used when measuring plateau heights.The maximum value of the KDE kernel is assigned asthe critical population, and these are collected in a sepa-rate file. Graphs of the KDE kernel and the total walkerpopulation are produced, and checked visually to ensurethe critical population was assigned correctly. Once allplateaus have been validated by visual inspection, thecritical populations are averaged and the standard errorcalculated.We note here that the FCIQMC critical walker popu-lations used throughout this work are from Ref. 29, andwere not recalculated for this work. III. RESULTS AND DISCUSSION
Calculations were performed on a variety of linear hy-drogen chains, and other small atoms and molecules us-ing the HANDE-QMC package, version 1.4. All simu-lations in this work were performed using a timestep of0 .
001 and a shift damping value of 0 .
30. Integral dumpfiles were generated using MOLPRO, in the form ofan FCIDUMP. The single particle eigenvalues for thesystems are then calculated in using an in-house codefrom the orbitals in the FCIDUMP according to stan-dard equations. These single particle eigenvalues arethen added to the FCIDUMP before the core Hamilto-nian energy.The equilibrium H n chains used in this study had abond length of 0.945110567 Angstroms, and the stretchedH n chains had a bond length of 1.270025398 Angstroms.These correspond to 1.786 a.u. and 2.4 a.u. respectively,which come from a previous study using auxiliary fieldquantum Monte Carlo. The H O system used a O-Hbond length of 0.975512 Angstroms and an H-O-H angleof 110.565 degrees. The CH system used a C-H bondlength of 1.087728 Angstroms, and a H-C-H bond angleof 109.47122 degrees. The bond lengths for the diatomicsystems are as follows: HF, 0.91622 Angstroms; NaH,1.885977 Angstroms; C , 1.27273 Angstroms; N , 2.068a.u. and stretched N , 4.2 a.u. These come from a pre-vious study using FCIQMC. The critical walker populations, or “plateau heights”,were measured by the Scott KDE method usingNumPy in Python3. The DMQMC calculations tomeasure critical populations for the H systems were per-formed with initial populations of 5 × and target pop-ulations of 5 × , and for H , the simulations used initialpopulations of 5 × and target populations of 5 × .These simulations were propagated to β =25. The IP-DMQMC and FCIQMC simulations for measuring thecritical populations were performed with an initial popu-lation of 1 and a target population of 5 × . All calcu-lations used the integer walker algorithm in all methodsto maintain comparability with the plateaus reported inthe first FCIQMC paper. The following results are now arranged as follows: inSec. III A, we begin by confirming the presence of the an-nihilation plateau and compare the critical walker pop-ulation in DMQMC to FCIQMC for stretched H . InSec. III B, we then explore the connection to the un-physical Hamiltonian that was found to contaminate pre-annihilation plateau energy in FCIQMC. Next, we gen-eralize our finding from Sec. III A to a wide range ofatomic and molecular systems in Sec. III C, and also ex-plore the interaction picture variant of DMQMC. As IP-DMQMC requires a target β value, we explore this nextin Sec. III D. Finally, we compare the initiator adaptionsto IP-DMQMC and FCIQMC in Sec. III E where we alsodiscuss importance sampling. β (Ha − )10 N w ( β ) (a) β (Ha − ) − − − − E ( β )( H a ) DMQMCft-FCI (b)
Figure 1. The total walker population, N w ( β ) (a) and energy, E ( β ) (b) from a single β loop, propagated to β = 25 forstretched H /STO-3G. The simulation was started at β = 0,and used a shift of 0.343 to ensure the plateau was exited by β = 25. In (b), the exact diagonalization (ft-FCI) is shown asa red dotted line. A. An example of a DMQMC annihilation plateau
We first begin by describing and then reproducing theoriginal finding of the DMQMC annihilation plateau,where we offer an example of an ab initio system. Thefirst paper on DMQMC described the sign problem in thismethod as similar to that of FCIQMC, due to the closesimilarities between the population dynamics within themethod. Of particular interest is the annihilation step,which is identical between the two methods, and is foundto be key in overcoming the sign problem, as describedearlier. One difference between the two methods is thatthe rate of annihilation is likely less frequent in DMQMCthan in FCIQMC because there are more density matrixelements than there are terms in the wavefunction coef-ficient vector. Blunt et al. suggested that because ofthis slower rate, a higher number of walkers would beneeded in DMQMC to overcome the sign problem – ap-proximately the square of the size of the FCIQMC criticalwalker population – and this observation was based on aHubbard model calculation.The annihilation plateau for stretched H /STO-3G isshown in Fig. 1. This plateau occurs after the first expo-nential growth, when the population reaches a system-specific population of walkers, as seen in Fig. 1(a) (forone β loop) between β = 0 and β = 5. When the spe-cific population of particles is reached, the spawning andannihilation rates are approximately equal, resulting inno population growth, i.e., the plateau. After exiting theplateau around β = 15, we observe a second exponentialgrowth phase. The plateau can almost always be visu-ally identified by its distinctive appearance, although inpractice, we have also automated this initial measure-ment (see Sec. II).When inspecting Fig. 1(b), the instantaneous energyestimate begins the simulation in reasonable agreementwith finite temperature full configuration interaction (ft-FCI) energy , but quickly thereafter the energy fluctuatesconsiderably. After the simulation exits the plateau re-gion, we see a return to agreement between the DMQMCenergy and the ft-FCI energy.We can also compare the plateau heights forstretched H /STO-3G (200 determinants) in DMQMCand FCIQMC. Here, the plateau height is measured at2 . × particles for DMQMC. This system inFCIQMC has a smaller plateau height, at only 2 . × particles. This is consistent with the description ofBlunt et al. , where the authors commented that theDMQMC plateau height is approximately the square ofthe FCIQMC plateau height. The plateau occurs be-tween β = 5 and β = 15 in this simulation and this tem-perature range is something that we cannot easily con-trol as an independent variable. Thus, while the criticaltemperature is something we could measure we generallyneglect it for this study.In summary, the annihilation plateau in DMQMC for ab initio systems follows previous observations basedon model Hamiltonians and the FCIQMC annihilationplateau. B. Connection to unphysical Hamiltonian and annihilationrate
The sign problem arises in DMQMC because spawningevents are affected by the sign of the Hamiltonian, H ik ,connecting two density matrix elements. In general, thesign of the matrix element H ik ( i (cid:54) = k ) can be positiveor negative. One way to think about how this arises isthat the Slater–Condon rules are applied by bringing theoccupied orbitals in determinant i and k into maximumcoincidence by permuting the electron indices with eachpermutation causing a change in sign. It follows that ρ kj can also have any sign. A sufficient number of walkersmust be present to allow for the efficient cancellation ofsigned spawning events arriving at ρ kj to resolve the signof ρ kj .Spencer et al. proposed that the sign problem inFCIQMC was: (1) due to an unphysical Hamiltonian( ˆ˜ H ) whose off-diagonal matrix elements have been whollynegated while leaving the magnitude unchanged, i.e. ,˜ H ik = δ ik H ik − (1 − δ ik ) | H ik | , where δ ik is the Kronecker delta, and (2) tending to be as severe as the energy ofthe dominant eigenvalue of the unphysical Hamiltonian.The authors found that these can be summarized by thefollowing equation for the critical walker population, N c : N c ≈ V max κ (13)Here, κ is the annihilation rate constant and V max is theenergy of the highest energy eigenstate of V = − ˜ˆ H ac-counting for the shift correlation energy and the HF en-ergy (i.e. V max = V + S + E HF ). The approximationin the equation refers to this being valid in the limit ofa small population and to first order in V max . Below,we test the same observations for DMQMC using thestretched H /STO-3G system.It is first useful to identify the sign structure of den-sity matrix for both the physical and unphysical Hamil-tonians. These are shown in Fig. 2(a) and Fig. 2(b),respectively for β = 3. It can be seen in this figure thatthese matrices differ in both the signs of their elementsas well as the distribution of the occupied elements. Inthe physical Hamiltonian, there is a mixture of both posi-tively and negatively signed elements, distributed denselyacross the entire matrix. The combination of the heav-ily signed and densely packed elements explains why thisinverse temperature is difficult to sample. In contrast,we see in the unphysical Hamiltonian matrix that onlynegatively signed elements exist, and is more evenly dis-tributed compared to the physical Hamiltonian matrix.Now, in the DMQMC simulations of both Hamiltonians,different dynamics are seen. For the physical Hamilto-nian, DMQMC exhibits a characteristic plateau shape asthe total population growth rises exponentially, pauses,and then resumes (Fig. 2(c)). Only when the popula-tion growth resumes does the growth of walkers on thediagonal of the density matrix start in earnest. In thedynamics of the simulation, we see that the walkers onthe diagonal tend to spawn and then die, depleting thediagonal population. It is only when enough of a popula-tion exists on the off-diagonal part of the density matrixand the sign structure has been established that the di-agonal population can be sustained. By contrast, for theunphysical Hamiltonian, DMQMC exhibits largely unin-terrupted growth in both the total population and thepopulation of walkers on the diagonal. Therefore the an-nihilation plateau, signalling an unresolved sign problem,occurs when ˜ˆ H (cid:54) = ˆ H .To show that the sign problem is also related to thedominant eigenvector of the unphysical Hamiltonian, ft-FCI results are shown in Fig. 2(d). We can see here thatin general, the energies obtained from the two Hamil-tonians are different, where the unphysical Hamiltonianenergy is lower than that of the physical Hamiltonianenergy. The one exception we see is at β = 0, when thetwo solutions are degenerate owing to the trace being thesame between the physical and unphysical Hamiltonian.In the low β regime is exactly where the dynamics ap-pear to be the most similar in terms of population growth − − − − − − − − − − − − − − (a) − − − − − − − − − − − − − − (b) β (Ha − )10 N o . o f P a r t i c l e s ( w a l k e r s ) N w ( β ), physical N Tr ( β ), physical N w ( β ), unphysical N Tr ( β ), unphysical (c) β (Ha − ) − − − − − E ( β )( H a ) ft-FCI, physical ft-FCI, unphysical (d) Figure 2. For the stretched H /STO-3G system we show the deterministic density matrices for β = 3 expressed as heatmapsfor (a) the physical Hamiltonian and (b) the unphysical Hamiltonian. In (a) and (b), blue corresponds to negatively signedelements, and red corresponds to positively signed elements. The darker the color, the larger the weight of the element (basedon a log scale, e.g. , -7 represents elements of 10 − ). In (c), for the same system, the population on the diagonal ( N Tr ) and totalwalker population ( N w ) are shown for the physical (blue) and unphysical (green) Hamiltonians from a single β loop simulationin DMQMC. In (d), the exact temperature dependent diagonalizations (ft-FCI) of the physical (red) and unphysical (green)Hamiltonians are shown. (Fig. 2(c)) which is consistent with the idea that the dom-inant eigenvalue of the unphysical Hamiltonian causes achange in the population dynamics.To provide further data to make this point, we scaledthe off-diagonal matrix elements linearly by a factor of C , starting from the true Hamiltonian, resulting in˜ H ik = δ ik H ik + C (1 − δ ik ) H ik (14)where H ik and δ ik follow previous definitions, with anadditional factor of a positive constant C . The plateauheight for low C follows a linear trend (Fig. 3) which fits the form of Eq. (13) as V max is linear in C for thissystem (assuming a constant κ ). This observation is alsoconsistent with that of Spencer et al. showing that theplateau height varies linearly with U/t in the Hubbardmodel where U is the on-site interaction strength and t is the hopping integral. Thus the analog to U in ourre-scaled molecular Hamiltonian is C .For completeness, the last component of the plateauexpression given in Eq. (13) we want to test for DMQMCis the dependence on the shift parameter, S . We collecteddata for the equilibrium H system shown in Fig. 4. It C N c ( w a l k e r s ) × Figure 3. The plateau height in DMQMC ( N c ) for equilib-rium H as a function of the scaling factor C , described inthe main text (Eq. (14)). Here, the plateau heights shownare averages from four β loops, and the values of C used are1 , . , , ,
10. The y = mx + b fit is N c = (5 . × ) × C ,where the y-intercept is assumed to equal zero. V max (Ha)0123 N c × Equilibrium-H Figure 4. The critical population ( N c , walkers) for equilib-rium H as a function of the energy value V max , which includesthe shift, S , as V max = V + S + E HF . The critical populationshere were obtained from averaging over 25 β loops. can be seen from these data that at low S , the plateauheight is linear in S which is consistent with the form ofEq. (13).In this section, we found that the unphysicalHamiltonian that contaminates pre-plateau unconvergedFCIQMC calculations is the same as the Hamiltonianthat causes the sign problem in DMQMC. This is impor-tant because, as in FCIQMC, the contamination preventsus from obtaining accurate energies; overcoming the an-nihilation plateau is the only way to prevent this fromoccurring. Thus, the annihilation plateau is the (mem-ory and time) bottleneck in converging DMQMC calcu-lations and allows us to explore the cost scaling of themethod in relation to FCIQMC. C. How the DMQMC and IP-DMQMC plateau heightsscale in relation to the FCIQMC plateau height
In Sec. III A, we observed that the plateau height inDMQMC was approximately the square of the plateauheight in FCIQMC for the stretched H system. In thissection, we attempt to generalize this observation to awide range of atomic and molecular systems for bothDMQMC and IP-DMQMC. To achieve this, we study therange of closed-shell systems that were previously con-sidered by Booth et al. (various atoms and moleculescomprised of first-row atoms) and supplement these with1D hydrogen chains. The latter set are of interest be-cause they are approximate analogs of the Hubbard mod-els, which are also used for plateau studies . Thus,the total test set is comprised of: Ne (aug-cc-pVDZ),H O (cc-pCVDZ), HF (cc-pCVDZ), NaH (cc-pCVDZ),C (cc-pVDZ), CH (cc-pVDZ), N (cc-pVDZ), stretchedN (cc-pVDZ), as well as stretched and equilibrium H n (STO-3G) for even n between 4 and 16 inclusive. The Beatom is excluded from the test set as it has no measurableannihilation plateau in IP-DMQMC. This test set repre-sents a variety of chemical systems including hetero- andhomonuclear diatomics with single and multiple bonds.Our preliminary observation was that the DMQMC andIP-DMQMC plateau heights were a system-dependentfraction of the size of the space similar to FCIQMC. Thismade it difficult to establish a specific trend with systemsize.Figure 5 shows DMQMC plateau heights plottedagainst their relative FCIQMC plateau heights for equi-librium and stretched H , and for equilibrium andstretched H . All of the other systems in our test sethad critical populations in DMQMC that were > × particles (our choice of the cutoff in population in ourexperimental design). What we see in this data is thatfor these four systems, the DMQMC plateau height is ap-proximately the square of the FCIQMC plateau height.Figure 5 also shows plateaus heights from the inter-action picture variant of DMQMC (IP-DMQMC) whichwas introduced in Sec. II B. By means of a brief sum-mary, IP-DMQMC differs from DMQMC because: (1) Atarget beta must be chosen (here, β = 25, with varia-tions in β considered in the next section); (2) the den-sity matrix is replaced with an auxiliary matrix ˆ f ( τ ) = e − ( β − τ ) ˆ H e − τ ˆ H ; and (3) the propagator changes in or-der to simulate ˆ f . Our calculations in Fig. 5 show thatthe IP-DMQMC plateau height is approximately equalto the plateau height in FCIQMC for the systems stud-ied here. For example, for the stretched H system, theIP-DMQMC plateau height is 2 . × particles, andthe FCIQMC plateau height is also 2 . × particles.This finding is remarkable, as it shows that the criticalwalker population in IP-DMQMC is directly related tothe same in FCIQMC.To further emphasize this, Fig. 6 shows the criticalwalker population in IP-DMQMC plotted next to the FCIQMC N c (walkers)10 N c ( w a l k e r s ) Benchmark Systems (IP-DMQMC)Equilibrium H n (DMQMC)Stretched H n (DMQMC)Equilibrium H n (IP-DMQMC)Stretched H n (IP-DMQMC) Figure 5. The plateau heights (N C ) for DMQMC (red) andIP-DMQMC (green and blue) simulations are shown with re-spect to the plateau height in FCIQMC, with both axes ona logarithmic scale for the benchmark systems from Booth etal. (circle), equilibrium H n chains (even n between 6 and16 inclusive, square symbol) and stretched H n chains (even n between 6 and 16 inclusive, x symbol). IP-DMQMC simula-tions used a target β = 25. Straight lines are plotted for both y = x (solid) and y = x (dashed) to help guide the eye. Theplateau heights were measured using the KDE method, andwere averaged over 25 simulations. The exception to this wasthe water molecule, where data collection follow the captionin Fig. 7 for consistency with that data set (see text for expla-nation). The FCIQMC critical walker populations are frompublished data. Error bars are shown, and, in some cases,are smaller than the size of the marker. size of the Slater determinant space in FCIQMC. It canbe seen that almost all of these systems have plateausheights lower than the number of determinants and,therefore, lower than the square root of the number ofelements in the density matrix.So why does IP-DMQMC have such a different plateauthan DMQMC? To analyze this, we can use Eq. (13). Forthe simulations of stretched H , we can find V max = 1 . V max is thesame or approximately the same then κ can be calcu-lated for each method. We then find that the κ valuesfor FCIQMC, DMQMC, and IP-DMQMC are 7 . × − ,5 . × − , and 7 . × − respectively. Here, we can seethat the IP-DMQMC rate of annihilation is the same asFCIQMC and approximately the square root of that inDMQMC i.e. IP-DMQMC encourages a similar rate ofannihilation as in FCIQMC. This explains the appear-ance of Fig. 5. The DMQMC plateau height scales howwe would expect it to: the annihilation rate is the squareroot of FCIQMC because the size of space is squared, N e /a u g - cc - p V D Z H O / cc - p V D Z H F / cc - p C V D Z N a H / cc - p C V D Z C / cc - p V D Z C H / cc - p V D Z E q u ili b r i u m - N / cc - p V D Z S t r e t c h e d - N / cc - p V D Z E q u ili b r i u m - H / S T O - G S t r e t c h e d - H / S T O - G N o . o f S t a t e s N c (IP-DMQMC) N dets Figure 6. The critical population ( N c , blue) in IP-DMQMCfor the atoms and molecules in our test set, compared to theestimated size of space ( N dets , red). These N c were collectedin the same way as Fig. 5. and therefore the plateau height is the square of that inFCIQMC. In contrast, we find the IP-DMQMC annihi-lation is identical to FCIQMC.We attribute the enhanced annihilation in IP-DMQMCto the more compact form of ˆ f in comparison to ˆ ρ ; thelatter is spread fully throughout the space. This can beillustrated by comparing the equations that are simu-lated in each method, that were presented in Sec. II. InDMQMC, d ˆ ρdβ = − ( ˆ H ˆ ρ + ˆ ρ ˆ H ) is propagated through thesimulation, whereas in IP-DMQMC d ˆ fdτ = ˆ H ˆ f − ˆ f ˆ H ispropagated. As defined previously, ˆ H is diagonal in ourchosen basis set, increasing the compactness compared toDMQMC, where neither ˆ ρ or ˆ H are diagonal matrices.Furthermore, the IP-DMQMC propagator only spawnsdown one of either the rows or columns of the density ma-trix and it appears that this contributes to IP-DMQMCbeing similar to FCIQMC in its annihilation.The key finding of this section is that, while theDMQMC plateau height scales with the square of theFCIQMC plateau height, the IP-DMQMC plateau heightscales linearly. It is possible that the saving in the costof IP-DMQMC only applies when a large β correspond-ing to the ground state is used as the target β ; the nextsection tests whether this is the case for all β .0 β (Ha − )4 . × . × . × × . × N c ( w a l k e r s ) Figure 7. The critical population ( N c ) in IP-DMQMC asa function of the target β value, averaged over 10 β loops( N β ) for H O/cc-pCVDZ. As the β is decreased, simulationsprogress through the annihilation plateau to different extents,causing an increase in the stochastic error. D. The IP-DMQMC plateau at different β values IP-DMQMC differs from DMQMC in that a specific(target) β must be specified and then each β has a uniquesimulation. This means that there is the potential fordependence of the IP-DMQMC plateau height on thespecified β . In this section, we explore whether the IP-DMQMC plateau heights across intermediate β valuesare the same as the FCIQMC plateau height.We first test this dependence using the H O system,and we present critical populations as a function of tar-get β in Fig. 7. As each target β progresses throughthe plateau to a varying extent by the time the target isreached, to ensure a fair test we used a wall time limitof 4 hours instead. We found that for all β values simu-lated the plateau heights are between 4 × and 5 × ,showing evidence that the IP-DMQMC critical popula-tion is not strongly β dependent. The FCIQMC criticalpopulation for this system was found to be 4 . × ,and so these results confirm that the IP-DMQMC crit-ical population is approximately the same as FCIQMCacross all β values.When moving to a different system, namely stretchedH , we found that it was more challenging to measurea plateau height at intermediate β values directly, aschanges in input parameters would be required to makesure that a plateau even exists by the time the target β is reached. When comparing the two systems, thewalker growth is slower in the simulations of the stretchedH system compared to the growth in the simulations ofH O, and we expect this may be why the plateau heightsin H O can measured directly, whereas they cannot bemeasured directly in stretched H . By means of an alter-native, we instead study how the DMQMC energy con-verges to the exact result with walker number and howthis convergence changes a function of the target β . Wedo so as an alternative to changing input parameters inthe IP-DMQMC simulations. β (Ha − ) − E n e r g y D i ff e r e n ce ( H a ) N w = 10 , N β = 2 . × N w = 10 , N β = 2 . × N w = 10 , N β = 2 . × N w = 10 , N β = 2 . × N w = 10 , N β = 2 . × (a) N β − − − − − E n e r g y ( H a ) ft-FCIIP-DMQMC (b) Figure 8. For the stretched H system, the absolute en-ergy difference from exact diagonalization (ft-FCI) is shownin panel (a) for IP-DMQMC at four different target popula-tions: 10 (orange circle), 10 (blue star), 10 (teal diamond),10 (green x symbol), and 10 (purple cross). The closer tozero the energy difference is, the more converged we considerthe energy. Error bars are shown but may be smaller than thesize of the marker in some cases.In panel (b) we show the IP-DMQMC energy (blue) at β =7 as a function of the numberof β loops in the simulation ( N β ), showing the energy with atarget population of 10 particles does converge to the ft-FCIenergy (black dashed line). We have, so far, established that for the stretchedH system the critical populations in DMQMC andIP-DMQMC (target β = 25) are 2 . × and2 . × particles, respectively. Simulations were per-formed at varying populations between 10 and 10 walk-ers with variable shift used throughout the simulation.We note that the random initialization algorithm of IP-DMQMC means the population can vary slightly fromthe starting population. For this test, the number of β loops was reduced as the walker number was increased so1the error stated is then the stochastic error for a givencomputational cost. Looking first at the IP-DMQMCresults (Fig. 8), the β = 10 results are consistent withthe plateau height being ∼ , as all walker numbersbetween 10 and 10 converge. That 10 walkers con-verges is somewhat at odds with a plateau height of being2 . × . However, we know that being in fixed shiftmode adds a factor linear in S to the plateau so it seemslikely that variable shift convergence is slightly improvedover fixed shift plateaus. When looking at intermedi-ate β values the stochastic error ( σ ) rises considerably,complicating the picture of convergence. All β values aretechnically converged within stochastic error (i.e., ∼ σ ),but the errors are of the order 1 Ha at β = 7. To checkthat this does converge without systematic bias, we canenlarge the number of beta loops and look at where theenergy converges to in the N β → ∞ limit. This is shownin Fig. 8(b), demonstrating that the energy does con-verge. Thus, these data support the observation that theplateau height for IP-DMQMC is ∼ —the same asthat in FCIQMC.To provide further corroboration, we can examine thesame trends in DMQMC. This is shown in Fig. 9. Here,similar patterns are seen to IP-DMQMC. When 10 walk-ers are used, a trend in the energy can be seen which issimilar to when 10 walkers are used in IP-DMQMC,which supports the idea that the IP-DMQMC plateauheight is the square root the height of the DMQMCplateau. When 10 walkers are used for DMQMC, wecan see that the energies still agree within stochastic er-ror with the number of β loops used. However, whenwe use more β loops we can see that this apparent con-vergence is false and the DMQMC result converged withrespect to stochastic error falls significantly away fromthe exact result. These observations are consistent withDMQMC having a plateau height of ∼ walkers.This result is significant because it shows that in IP-DMQMC, the target population of 10 is sampling thedensity matrix associated with the physical Hamiltonian,described in Sec. III B. This contrasts our finding forDMQMC, where because the energy converges to some-thing different than the ft-FCI energy, we can concludethat the simulation is not sampling the density matrixassociated with the physical Hamiltonian. Together thisgives corroborating evidence that the critical populationin IP-DMQMC is lower than that in DMQMC, and is infact cheaper than DMQMC across all β values.From the evidence presented here, it seems that IP-DMQMC does not follow FCIQMC scaling just becausethe β values we were using in the previous section corre-sponded to ground states. While the number of β loopsrequired to converge an intermediate β simulation is, ad-mittedly, very large, for the stretched H system studiedhere, the simulation is indeed free of systematic bias whenthe population is above the critical population we mea-sured. For calculations that were performed at the low-temperature (high β ) limit, where the ground state is oc-cupied, the annihilation plateau follows the explanation β (Ha − ) − E n e r g y D i ff e r e n ce ( H a ) N w = 10 , N β = 2 . × N w = 10 , N β = 2 . × (a) N β − − − − E n e r g y ( H a ) ft-FCIDMQMC (b) Figure 9. For the stretched H system, the absolute en-ergy difference from exact diagonalization (ft-FCI) is shownin panel (a) for DMQMC at two different target populations:10 (empty, orange edge circle), and 10 (empty, teal edge di-amond). The closer to zero the energy difference is, the moreconverged we consider the energy. Error bars are shown butmay be smaller than the size of the marker in some cases.For N w = 10 , the energy differences at β = 7 and β = 8,are outside the range of the plot, and are -1.814 Ha and -9.808 Ha, respectively. In panel (b) we show the DMQMCenergy (blue) at β =7 as a function of the number of β loopsin the simulation ( N β ), showing that the energy converges toa different value from the ft-FCI energy (black dashed line). posited above related to the annihilation rate κ . Then,as β is lowered, the sign problem remains remarkablyconstant while the sampling problem (and its stochasticerror) becomes more severe (as there are more elementsto be sampled). Our next two sections explore possi-ble ways to improve and resolve the stochastic samplingproblem which emerges. E. Importance sampling and the initiator approach
Importance sampling was developed along withDMQMC to prevent the escape of walkers from the trace2 H ( S T O - G ) H ( S T O - G ) N e ( a u g - cc - p V D Z ) H F ( cc - p C V D Z ) N c N detIP-DMQMCIS IP-DMQMC DMQMCIS DMQMC Figure 10. The plateau height ( N c ) on a logarithmic scale forequilibrium H , equilibrium H , Ne and HF in IP-DMQMC(blue), IS IP-DMQMC (green), DMQMC (dark blue) and ISDMQMC (dark green). For comparison, the number of deter-minants ( N det ) for each system are shown (red). The criticalpopulations in the importance sampling simulations were av-eraged over 25 β loops. Importance sampling generally raisesthe plateau height with the exception of Ne atom where theplateau heights agree within error bars. of the simulation and improve statistical sampling. Wefind that, in general, that while the importance samplingapproach does keep walkers on the diagonal of the den-sity matrix, it also generally slightly raises the heightof the annihilation plateau (Fig. 10). So, while it haspromise in terms of converging the energy with reducednoise (due to having an increased trace population) wedo not investigate it any further here.The second method which maintains population on thediagonal is the initiator adaptation (Sec. II D), whichis popular in FCIQMC because it removes the require-ment that the simulation has to have a total walker num-ber greater than the critical walker population ( i.e. theplateau is removed), introducing only a modest error.Unfortunately, the removal of the plateau means we can-not compare how i-FCIQMC and i-IP-DMQMC scale us-ing this measure alone. We can, therefore, use a previousstudy of i-FCIQMC where a walker population thresh-old measure was defined as the total number of walkersin the simulation required to have 50,000 walkers on theHartree–Fock determinant. This threshold is analogousto measuring the plateau height because it was shown fora variety of atoms that the energy did not vary after thisthreshold was reached and the simulation was convergedwith respect to stochastic sampling, which is the sameidea as the canonical method needing to reach a criti-cal walker population to obtain a converged energy. We considered but did not attempt a “growth witness” mea-sure ( G ) introduced by other authors. So to adapt thismeasure for DMQMC, we consider a threshold of 50,000walkers on the trace of the density matrix which controlsfor systematic and stochastic errors simultaneously. Asin the previous section, we compare i-IP-DMQMC withi-FCIQMC.In Fig. 11, we show the total walker population at thesimulation iteration at which the population thresholdwas met on the diagonal of the matrix (or HF deter-minant for i-FCIQMC), and we will refer to this walkervalue as N thresh throughout this section. We find thatthe i-IP-DMQMC cost in terms of walker number is the same as i-FCIQMC for low temperatures ( i.e. , β =10),which makes sense because we are simulating the groundstate when our choice of target β is large.In the intermediate temperature range ( β =2 and β =5),we find that for the smaller hydrogen chains, the costis slightly higher in i-IP-DMQMC, but as the lengthof the chain is increased, i-IP-DMQMC returns to be-ing approximately the same cost as i-FCIQMC. At thetwo higher temperatures ( β =1 and β =0.1), we find thecost is lower than that of i-FCIQMC. We attribute thisto the initiator adaption itself. This variation tends tokeep walkers on the diagonal of the density matrix forDMQMC and at high temperatures more particles on thediagonal is closer to the physical solution, making it eas-ier for DMQMC to simulate. In general, these data showthat the initiator approximation in IP-DMQMC controlsthe walker population in a similar manner to the initia-tor approximation in FCIQMC. This gives us confidencethat the initiator approximation can be used in futureapplications. IV. CONCLUSIONS
DMQMC has been shown to be a promising methodfor finite temperature applications, and in this work wehave confirmed that DMQMC (especially in its interac-tion picture variant) shows the potential to be as ef-fective for finite temperature work as FCIQMC is forground state simulations. We confirmed that the criti-cal walker population in DMQMC scales as the squareof that in FCIQMC. Additionally, we found that thecritical walker population in IP-DMQMC is the sameas that of FCIQMC across all β values. The latter isa very exciting result, as it shows that we can obtaina temperature-dependent energy at roughly the samememory and walker cost as FCIQMC, allowing us totreat systems with IP-DMQMC that cannot be treatedby DMQMC. With respect to the critical walker pop-ulation, this implies that IP-DMQMC has more utilitycompared to DMQMC, because a smaller population ofparticles is required to obtain the density matrix associ-ated with the physical Hamiltonian. Finally, we showedthat the initiator adaption with IP-DMQMC performs ina similar or better way than i-FCIQMC, again allowing3 N thresh , i-FCIQMC10 N t h r e s h ,i - I P - D M Q M C β = 10 β = 5 β = 2 β = 1 β = 0.1 10 . × N . . × N . . × N . . × N . . × N − . Figure 11. The total walker population at the simulationiteration at which the population threshold was met on thediagonal of the matrix (HF determinant for i-FCIQMC) inthe i-IP-DMQMC simulation ( N thresh ) as a function of thesame in the i-FCIQMC simulation ( N thresh ), for the equilib-rium H n systems studied (n= 4, 6, 8, 10, 12). i-IP-DMQMCsimulations were run with increasing target populations start-ing from 1 × up to 5 × , with initial populations of 10particles, and with target β values of 0.1 (magenta x sym-bols), 1 (blue circles), 2 (green diamonds), 5 (gold squares)and 10 (purple + symbols). The y = 10 b × x m fits are shownin the legend, on the same line as the marker correspondingto the β value. us to expand upon the systems we can treat with thismethod.As such, we now know that IP-DMQMC will be moreuseful than DMQMC for systems with a severe sign prob-lem. One disadvantage of using IP-DMQMC, which wedid not explore here, is that IP-DMQMC requires sep-arate simulations to obtain energies for different inversetemperature values. Whether the computational over-head is then more expensive to obtain a full β spectrumin IP-DMQMC compared to DMQMC is still an openquestion. We note in passing that we did not explore theconnection between this observation and the Krylov pro-jected FCIQMC, as we felt that this was beyond thescope of this work.Overall, this strongly suggests that the IP-DMQMCalgorithm has the same potential as FCIQMC, and givesa focus for future development. A natural place for fu-ture work to begin is to explore the uses (and limitations)of the initiator approach in a systematic way as well asexamining ways to modify and lower the IP-DMQMCplateau height. For example, as it is known that basissets rotations do affect the plateau height, we are in-clined to explore basis sets that are optimized for a giventemperature. V. ACKNOWLEDGEMENTS
Research was primarily supported by the U.S. Depart-ment of Energy, Office of Science, Office of Basic EnergySciences Early Career Research Program (ECRP) underAward Number DE-SC0021317 (calculations and analy-ses by WZV and HRP). This work was also supported bythe University of Iowa through start-up funding (researchfacilitation by SKR, computer time). For the purposesof providing information about input options for the cal-culations used, files will be deposited with Iowa ResearchOnline (IRO) with a reference number [to be inserted atproduction].
REFERENCES H. R. Petras, S. K. Ramadugu, F. D. Malone, and J. J. Shepherd,J. Chem. Theory Comput. , 1029 (2020). S. Mukherjee, F. Libisch, N. Large, O. Neumann, L. V. Brown,J. Cheng, J. B. Lassiter, E. A. Carter, P. Nordlander, and N. J.Halas, Nano Lett. , 240 (2013). L. Zhou, C. Zhang, M. J. McClain, A. Manjavacas, C. M.Krauter, S. Tian, F. Berg, H. O. Everitt, E. A. Carter, P. Nord-lander, and N. J. Halas, Nano Lett. , 1478 (2016). G. Mazzola, R. Helled, and S. Sorella, Phys. Rev. Lett. (2018), 10.1103/PhysRevLett.120.025701. E. Gull, O. Parcollet, and A. J. Millis, Phys. Rev. Lett. ,216405 (2013). Y. Liu, M. Cho, and B. Rubenstein, J. Chem. Theory Comput. , 4722 (2018). Y. Liu, T. Shen, H. Zhang, and B. Rubenstein, J. Chem. TheoryComput. , 4298 (2020). D. M. Ceperley, Journal of Statistical Physics , 1237 (1991). D. M. Ceperley, Phys. Rev. Lett. , 331 (1992). T. Dornheim, S. Groth, J. Vorberger, and M. Bonitz, Phys. Rev.Lett. (2018), 10.1103/PhysRevLett.121.255001. J. P. F. LeBlanc, S. Li, X. Chen, R. Levy, A. E. Antipov, A. J.Millis, and E. Gull, Phys. Rev. B (2019), 10.1103/Phys-RevB.100.075123. G. Sanyal, S. H. Mandal, and D. Mukherjee, Chem. Phys. Lett. , 55 (1992). W. Li and P. Piecuch, The Journal of Physical Chemistry A ,6721 (2010). X. He, S. Ryu, and S. Hirata, J. Chem. Phys. , 024702(2014). A. A. Rusakov and D. Zgid, J. Chem. Phys. , 054106 (2016). A. E. Doran and S. Hirata, J. Chem. Theory Comput. , 6097(2019). D. Neuhauser, R. Baer, and D. Zgid, J. Chem. Theory Comput. , 5396 (2017). D. Hirschmeier, H. Hafermann, E. Gull, A. I. Lichtenstein,and A. E. Antipov, Phys. Rev. B (2015), 10.1103/Phys-RevB.92.144409. B. Hirshberg, M. Invernizzi, and M. Parrinello, J. Chem. Phys. , 171102 (2020). A. Yilmaz, K. Hunger, T. Dornheim, S. Groth, and M. Bonitz,J. Chem. Phys. , 124114 (2020). T. Dornheim, S. Groth, and M. Bonitz, Contributions to PlasmaPhysics , e201800157 (2019). G. Harsha, T. M. Henderson, and G. E. Scuseria, J. Chem.Theory Comput. (2019), 10.1021/acs.jctc.9b00744. G. Harsha, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. , 154109 (2019). P. Shushkov and T. F. Miller, J. Chem. Phys. , 134107 (2019). A. F. White and G. K.-L. Chan, J. Chem. Theory Comput. ,5690 (2018). A. F. White and G. Kin-Lic Chan, J. Chem. Phys. , 224104(2020). F. Hummel, J. Chem. Theory Comput. , 6505 (2018). A. Roggero, A. Mukherjee, and F. Pederiva, Phys. Rev. B ,115138 (2013). G. H. Booth, A. J. W. Thom, and A. Alavi, J. Chem. Phys. , 054106 (2009). G. H. Booth, D. Cleland, A. J. W. Thom, and A. Alavi, J. Chem.Phys. , 084104 (2011). J. S. Spencer, N. S. Blunt, and W. M. C. Foulkes, J. Chem.Phys. , 054110 (2012). M. Kolodrubetz and B. K. Clark, Phys. Rev. B , 075109 (2012). D. Cleland, G. H. Booth, and A. Alavi, J. Chem. Phys. ,041103 (2010). F. D. Malone, N. Blunt, E. W. Brown, D. Lee, J. Spencer,W. Foulkes, and J. J. Shepherd, Phys. Rev. Lett. , 115701(2016). J. E. Deustua, I. Magoulas, J. Shen, and P. Piecuch, J. Chem.Phys. , 151101 (2018). N. S. Blunt, J. Chem. Phys. , 174103 (2019). K. Ghanem, A. Y. Lozovoi, and A. Alavi, J. Chem. Phys. ,224108 (2019). K. Ghanem, K. Guther, and A. Alavi, J. Chem. Phys. ,224115 (2020). R. J. Anderson and G. H. Booth, J. Chem. Phys. , 184103(2020). E. Vitale, A. Alavi, and D. Kats, J. Chem. Theory Comput. ,5621 (2020). R. J. Anderson, T. Shiozaki, and G. H. Booth, J. Chem. Phys. , 054101 (2020). G. Li Manni, W. Dobrautz, and A. Alavi, J. Chem. TheoryComput. , 2202 (2020). H. R. Petras, D. S. Graham, S. K. Ramadugu, J. D. Goodpaster,and J. J. Shepherd, J. Chem. Theory Comput. , 5332 (2019). W. Dobrautz, S. D. Smart, and A. Alavi, J. Chem. Phys. ,094104 (2019). N. S. Blunt, A. J. W. Thom, and C. J. C. Scott, J. Chem. TheoryComput. , 3537 (2019). H. Luo and A. Alavi, J. Chem. Theory Comput. , 1403 (2018). N. S. Blunt, J. Chem. Phys. , 221101 (2018). G. Li Manni, S. D. Smart, and A. Alavi, J. Chem. Theory Com-put. , 1245 (2016). N. M. Tubman, J. Lee, T. Y. Takeshita, M. Head-Gordon, andK. B. Whaley, J. Chem. Phys. , 044112 (2016). N. S. Blunt, S. D. Smart, J. A. F. Kersten, J. S. Spencer, G. H.Booth, and A. Alavi, J. Chem. Phys. , 184107 (2015). J. S. Spencer, N. S. Blunt, S. Choi, J. Etrych, M.-A. Filip,W. M. C. Foulkes, R. S. T. Franklin, W. J. Handley, F. D. Mal-one, V. A. Neufeld, R. Di Remigio, T. W. Rogers, C. J. C. Scott,J. J. Shepherd, W. A. Vigor, J. Weston, R. Xu, and A. J. W. Thom, J. Chem. Theory Comput. , 1728 (2019). K. Guther, R. J. Anderson, N. S. Blunt, N. A. Bogdanov, D. Cle-land, N. Dattani, W. Dobrautz, K. Ghanem, P. Jeszenszki,N. Liebermann, G. L. Manni, A. Y. Lozovoi, H. Luo, D. Ma,F. Merz, C. Overy, M. Rampp, P. K. Samanta, L. R. Schwarz,J. J. Shepherd, S. D. Smart, E. Vitale, O. Weser, G. H. Booth,and A. Alavi, J. Chem. Phys. , 034107 (2020). T. Shen, Y. Liu, Y. Yu, and B. M. Rubenstein, J. Chem. Phys. (2020), 10.1063/5.0026606. N. S. Blunt, T. W. Rogers, J. S. Spencer, and W. M. C. Foulkes,Phys. Rev. B , 245124 (2014). F. D. Malone, N. S. Blunt, J. J. Shepherd, D. K. K. Lee, J. S.Spencer, and W. M. C. Foulkes, J. Chem. Phys. , 044116(2015). D. Cleland, G. H. Booth, and A. Alavi, J. Chem. Phys. ,024112 (2011). D. W. Scott,
Multivariate density estimation: theory, practice,and visualization , second edition ed. (Wiley, Hoboken, New Jer-sey, 2014). J. J. Shepherd, G. E. Scuseria, and J. S. Spencer, Phys. Rev. B , 155130 (2014). P. J. Knowles and N. C. Handy, Computer Physics Communica-tions , 75 (1989). A. Szabo and N. S. Ostlund,
Modern quantum chemistry: intro-duction to advanced electronic structure theory (Dover Publica-tions, Mineola, N.Y, 1996). C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers,P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg,N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk,M. Brett, A. Haldane, J. F. del R’ıo, M. Wiebe, P. Peterson,P. G’erard-Marchant, K. Sheppard, T. Reddy, W. Weckesser,H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature , 357(2020). IP-DMQMC is currently limited to treat only systems with M s = 0. M. Yang, E. Pahl, and J. Brand, J. Chem. Phys. , 174103(2020). N. Blunt, A. Alavi, and G. H. Booth, Phys. Rev. Lett.115