⇦Chinese Physics Letters, 2017, Vol. 34, No. 8, Article code 080201⇨ Validation of the Ability of Full Configuration Interaction Quantum Monte Carlo for Studying the 2D Hubbard Model * Su-Jun Yun(贠素君)1,3**, Tie-Kuang Dong(董铁矿)2, Shi-Ning Zhu(祝世宁)3 Affiliations 1School of Electronic Engineering, Nanjing Xiaozhuang University, Nanjing 211171 2Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008 3School of Physics, National Laboratory of Solid State Microstructures, Nanjing University, Nanjing 210093 Received 17 February 2017 ^*Supported by the Natural Science Foundation for Colleges and Universities of Jiangsu Province under Grant No 16KJB140008, the National Natural Science Foundation of China under Grant Nos 11447204 and 11647164, the Natural Science Foundation of Jiangsu Province under Grant No BK20151079, and the Scientific Research Foundation of Nanjing Xiaozhuang University under Grant No 2015NXY34.
^**Corresponding author. Email: yunsujun@163.com Citation Text: Yun S J, Dong T K and Zhu S N 2017 Chin. Phys. Lett. 34 080201 Abstract To validate the ability of full configuration interaction quantum Monte Carlo (FCIQMC) for studying the 2D Hubbard model near half-filling regime, the ground state energies of a $4\times4$ square lattice system with various interaction strengths are calculated. It is found that the calculated results are in good agreement with those obtained by exact diagonalization (i.e., the exact values for a given basis set) when the population of psi particles (psips) is higher than the critical population required to correctly sample the ground state wave function. In addition, the variations of the average computational time per 20 Monte Carlo cycles with the coupling strength and the number of processors are also analyzed. The calculated results show that the computational efficiency of an FCIQMC calculation is mainly affected by the total population of psips and the communication between processors. These results can provide useful references for understanding the FCIQMC algorithm, studying the ground state properties of the 2D Hubbard model for the larger system size by the FCIQMC method and using a computational budget as effectively as possible. DOI:10.1088/0256-307X/34/8/080201 PACS:02.70.Ss, 71.10.Fd, 02.70.Uu © 2017 Chinese Physics Society Article Text Among various methods that have been developed to investigate the quantum many-body systems, the quantum Monte Carlo (QMC) methods have been attracting more and more attention in recent years. Though QMC calculations are very resource-consuming, the extremely high accuracy has inspired many scientists to search for effective algorithms. In the past decades, several QMC methods have been developed to treat various systems.[1-5] As for many-electron systems, the projector quantum Monte Carlo (PQMC) methods have achieved great successes, such as Green's function quantum Monte Carlo (GFQMC), diffusion quantum Monte Carlo (DQMC), and auxiliary-field quantum Monte Carlo (AFQMC). In particular, the DQMC with the fixed-node approximation has been implemented in the software CASINO. This software can be used to calculate the ground state properties of solid state materials.[6] PQMC methods will suffer from a sign problem for fermionic systems. The sign problem stems from the antisymmetry property of many-electron wave functions when exchanging two electrons. Therefore, the wave functions of many-electron systems have both positive and negative amplitudes. The minus sign of wave function reduces the rate of convergence exponentially in real Monte Carlo calculations. This problem is regarded as one of the most important unsolved problems in the field of computational physics and chemistry. To obtain physically meaningful results, some approximated methods have been introduced, such as the fixed-node diffusion Monte Carlo method[7] and the phaseless auxiliary-field quantum Monte Carlo.[8] In principle, the exact ground state wave function can be obtained by the full configuration interaction (FCI) method in quantum chemistry, which is also known as the exact diagonalization for physicists. The FCI method converts the Schrödinger equation into a matrix eigenvalue problem, which can be solved by exact diagonalization. However, the computational cost of FCI scales factorially with the size of configuration space, thus FCI can only be used to investigate small systems. In 2009, Booth et al. combined the spirit of PQMC and the FCI method to develop the full configuration interaction quantum Monte Carlo (FCIQMC).[9] FCIQMC inherits some features from both PQMC and FCI. Similar to PQMC, FCIQMC performs a long-time integration of the imaginary-time Schrödinger equation, and similar to the FCI method, FCIQMC works in a determinant space. The weight of a given determinant is represented by the number of a given type of virtual particle with a charge ($q=\pm1$) on it, which is called the psi particle (psip). Each psip lives on a fixed determinant in which it was born until it is removed from the simulation. The ground state wave function is projected out through the evolution of population dynamics of psips rather than random walk, though psip is also called the walker in some studies. The FCIQMC algorithm contains three stages for each simulation cycle: spawning, death, and annihilation. It has been shown that it is the annihilation of opposite charged psips that allows FCIQMC to converge to the true ground state wave function. Here annihilation means a pair of psips with opposite charges on the same determinant are removed together from the simulation. This discrete sampling of the wave function allows efficient annihilation to take place between positive and negative psips on the same determinant. At present, the FCIQMC method has been proved much more successful than the DQMC method for computing the energies of small systems and systems where the sign problem is not severe. However, FCIQMC is a 'young' method, some questions still remain concerning how to use it effectively and automatically, how large a system can be calculated for a given computational resource, and so on. The FCIQMC algorithm is not a 'black box' method and the manner in which sampling is performed can impact the required computational resources. To understand how the FCIQMC method can be used to improve the fermion sign problem efficiently, and how to use the computational resources efficiently, we will take the 2D repulsive Hubbard model as an example in our discussions. The reasons for selecting the Hubbard model as a research subject are two-fold. On the one hand, the Hubbard model has been widely used to investigate the magnetism, metal-insulator transitions,[10] and superconductivity.[11] On the other hand, exact diagonalization results for some lattices are available, which provides a perfect benchmark to validate the ability of FCIQMC. Since there are many excellent reviews about the Hubbard model, we will give a brief description of it for the convenience of discussions. The single-band Hubbard model is defined by the Hamiltonian $$ \hat{H}=-t \sum_{\langle i j \rangle \sigma}C^{+}_{i \sigma} C_{j \sigma}+U \sum_{i} n_{i\uparrow} n_{i\downarrow},~~ \tag {1} $$ where $C^{+}_{i \sigma}$ ($C_{i \sigma}$) creates (annihilates) an electron with spin $\sigma$ on site $i$, and $n_{i\uparrow}$ ($n_{i\downarrow}$) is the particle number operator with spin $\uparrow$($\downarrow$). Since the aim of this work is to validate the ability of FCIQMC in investigating the Hubbard model, we select a relatively small system. As you know, the computation costs of QMC scale badly with the size of system. In this work, we restrict our discussions to the repulsive Hubbard model ($U > 0$) defined on a two-dimensional $4 \times 4$ square lattice, consider interaction strengths ranging from $ U=1.0$ to $U=10.0$, and focus on near half-filling regime ($N_{\rm e}=14$ ($7\uparrow7\downarrow$)), where $t$ is used as the unit of all energies. We briefly review the FCIQMC method, which was discussed in more detail in Refs. [9,12]. The FCIQMC algorithm is a type of PQMC. For each simulation cycle the wave function is updated by applying the following projection operator $$ \hat{P}={\boldsymbol I}-(\hat{H}-S{\boldsymbol I} )\delta \tau,~~ \tag {2} $$ where $S$ is an energy offset which is adjusted to control the population of psips after a critical population of psips required to correctly sample the ground state wave function, $\hat{H}$ is the Hamiltonian of the Hubbard model as shown in Eq. (1), ${\boldsymbol I}$ is the unit operator, and $\delta \tau$ is the step of imaginary time. On every iteration the representation of the wave function, $\varphi(\tau)$, at imaginary time $\tau$ is updated to $\varphi (\tau+\delta\tau)$ as follows: $$ \varphi (\tau+\delta\tau)=\hat{P}\varphi(\tau)=[{\boldsymbol I}-(\hat{H}-S {\boldsymbol I})\delta\tau]\varphi(\tau).~~ \tag {3} $$ As long as the initial wave function $\varphi(0)$ has a non-zero component of the ground state wave function, the ground state will manifest itself once $\tau$ is large enough, assuming that $\delta \tau$ is sufficiently small and $S$ is carefully controlled. In FCIQMC, the configuration interaction wave function can be written as $$ \varphi=\sum_{i} c_i |D_i \rangle,~~ \tag {4} $$ where ${|D_i\rangle}$ is a set of the Slater determinants of size $N_{\det}$, which is formed from $N_{\rm elec}$ electrons and $N_{\rm s.o.}$ spin–orbitals. The wave function coefficient $c_i=\langle D_{i}|\varphi(\tau)\rangle$ is represented by a number of (signed) unit weights located on $|D_{i}\rangle$. We call a single unit a psi-particle or psip. Similar to FCI, FCIQMC satisfies the variational principle, the optimal coefficients $\{c_{0i}\}$ gives the lowest energy, i.e., the ground state energy eigenvalue. To make the subsequent notation as simple as possible, we define a 'transition matrix' $$ {\boldsymbol T}=-(\hat{H} - S {\boldsymbol I}).~~ \tag {5} $$ Then Eq. (3) can be rewritten as $$ c_{i}(\tau+\delta\tau)=\sum _{j}(\delta_{ij}+T_{ij}\delta \tau)c_{j}(\tau).~~ \tag {6} $$ At each time step $\delta \tau$, Eq. (6) is sampled in three stages: (1) The spawning step: each psip (with weight $w_i$) attempts to spawn a child psip on a randomly selected $|D_{j}\rangle$ with probability $p_{\rm s}(j|i)=\frac{|T_{ij}|\delta \tau}{p_{\rm gen}(j|i)}$, where $p_{\rm gen}(j|i)$ is the probability that $|D_{j}\rangle$ is selected, provided that the parent psip is on $|D_{i}\rangle$ and with sign ${\rm sgn}(-\langle D_{i}|\hat{H}|D_{j}\rangle \omega_i)$. (2) The death/clone step: each psip dies with probability $|p_{\rm d}(i)|$ if $p_{\rm d}(i) < 0$, where $p_{\rm d}(i)=-\delta \tau T_{ii}$. Instead, if $p_{\rm d}(i)> 0$ a copy of the psip is created with probability $p_{\rm d}(i)$. (3) The annihilation step: in the final stage, we loop over all psips including newly-spawned, cloned and surviving parent psips, and annihilate pairs of psips with opposite signs on the same determinant. The annihilation step preserves the expected distribution of psips and crucially prevents the growth of exponential noise and collapses onto the sign problem-free ground state. Typically, the initial value of energy shift $S$ is set to be the Hartree–Fock energy, which is larger (less negative) than the FCI energy and this causes the population $N_{\rm p}(\tau)$ (total number of psips) to grow exponentially. However, in this stage the population contains physical and unphysical components, and the unphysical one dominates the population.[12] Then $N_{\rm p}(\tau)$ will enter the plateau phase at a system-dependent population, during which the ground state sign structure emerges, the component of physical solution becomes more and more important. At some critical point the physical component exceeds the unphysical component. It is only after this point that $\varphi(\tau)$ becomes a stochastic representation of the wave function and $N_{\rm p}(\tau)$ grows exponentially again, but the growth rate is lower than that of the first stage. In actual simulations, it is unnecessary to keep the exponential growth of particle population. Instead, the shift is then adjusted to keep the population approximately constant according to the expression $$ S(\tau+A\delta\tau)=S(0) - \xi \log\frac{N_{\rm p}(\tau+A\delta\tau)}{N_{\rm s}},~~ \tag {7} $$ where $S(0)$ is the initial value of the shift (in this work, the Hartree–Fock energy), $N_{\rm s}$ is the population at the end of the equilibration phase, $\xi$ is the damping factor, which is usually fixed during a simulation. For the constant population mode, the mean value of $S$ provides a measurement of the ground state energy. In addition to the shift, the energy of the system can be estimated via the projected estimator $$ E_{\rm project}=\frac{\langle D_{0}|\hat{H}e^{-\hat{H}\tau}|D_{0}\rangle}{\langle D_{0}|e^{-\hat{H}\tau}|D_{0}\rangle}=H_{00}+\frac{\sum _{i\neq0} H_{0i} N_i}{N_0},~~ \tag {8} $$ where $N_{i}(\tau)$ is the (signed) number of psips on the $i$th determinant for the state $|\varphi(\tau)\rangle$, $N_0$ represents population of psips on the reference determinant $|D_{0}\rangle$ (here taken to be the Hartree–Fock determinant), $H_{00}=\langle D_{0}|\hat{H}|D_{0}\rangle$ is the energy of the reference determinant $|D_{0}\rangle$, and $H_{0i}=\langle D_{0}|\hat{H}|D_{i}\rangle$.

Fig. 1. Population dynamics [(a)–(c)] and energy estimators (d) for the 2D 16-site square lattice Hubbard model at $U=4.0$, and $N_{\rm e}=14$ ($7\uparrow7\downarrow$). The ground state wave function has momentum $k=(0,0)$ and the Hilbert space formed from all determinants of this momentum contains about $0.8\times10^7$ states. (a) Evolution of the total number of H psips in the initial stage of simulation. (b) Evolution of the total number of $H$ psips in simulation. (c) Variations of $-\sum H_{0j}N_{j}$ and the population of psips on the reference determinant $N_0$ with the simulation iterations. (d) The shift $S$ and the correlation energy $E_{\rm corre}=N_0^{-1} \sum _{j\neq0} H_{0j} N_j$ in simulation.

Now, we take $U=4.0$ and $N_{\rm e}=14$ ($7\uparrow7\downarrow$) for the $4\times4$ square lattice as an example to discuss the FCIQMC algorithm in detail, for which the reference energy (the Hartree–Fock energy) $H_{00}=-11.75$ in Eq. (8). The initial value of the shift is set as a constant value ($S=0$) and a small starting population ($N_{\rm p}(0)=100$), and hence the population grows exponentially. The rapid exponential growth process is very short, and then a plateau phase in the population appears. The height of plateau is about $10^7$ as shown in Fig. 1(a). After the plateau population of psips is reached, we need to adjust the shift $S$ according to Eq. (7). When the current total population of Hamiltonian particles $N_{\rm p}(\tau)$ varies with a very small margin (see Fig. 1(b)), the population of Hamiltonian particles $N_0$ and $-\sum H_{ij} N_j$ becomes nearly constant from about the 10000th iteration (see Fig. 1(c)). Once the population is stable, both the shift $S$ and the instantaneous projected energy vary around a fixed value, namely the ground state correlation energy (see Fig. 1(d)).

**Table 1.** The simulated results for different interaction strengths.
$U$	$\delta \tau $	Plateau	H psips	$\sum H_{0j}N_j$	$N_0$	$S$	$E_{\rm corre}$
1.0	0.002	$1\times10^7$	12048220(60)	$-109640(30)$	$240790(60)$	$-0.45529(1)$	$-0.45533(2)$
2.0	0.001	$1\times10^7$	17012900(900)	$-124100(100)$	$95710(70)$	$-1.2959(1)$	$-1.29673(5)$
3.0	0.001	$1\times10^7$	27545000(2000)	$-172210(80)$	$69250(30)$	$-2.4873(1)$	$-2.4870(5)$
4.0	0.001	$1\times10^7$	49784000(4000)	$-280400(100)$	$70180(40)$	$-3.9947(2)$	$-3.9954(6)$
6.0	0.001	$2\times10^7$	471100000(10000)	$-2211400(500)$	$283650(60)$	$-7.79630(6)$	$-7.7963(3)$
8.0	0.0005	$3\times10^7$	363290000(20000)	$-1362700(400)$	$110150(30)$	$-12.3689(2)$	$-12.3710(3)$
10.0	0.0005	$3.5\times10^7$	1201600000(50000)	$-3925000(2000)$	$225200(100)$	$-17.4322(2)$	$-17.4308(1)$

In Table 1, we show the calculational results for different interaction strengths ($U=1$, 2, 3, 4, 6, 8 and 10 in units of $t$). In this table, the plateau denotes the plateau height, H psips denotes the total population of Hamiltonian particles, $N_0$ is the population of Hamiltonian particles on the reference determinant, $S$ is the mean value of energy shift, and $E_{\rm corre}$ is the correlation energy. From Table 1, one can find some characteristics. On the one hand, with the increase of the interaction strength $U$, the plateau height becomes higher, and the total population of psips becomes larger. When $U>6.0$, the population of H psips is too large to calculate if $\delta \tau$ is still 0.001, thus $\delta \tau=0.0005$ is taken in the case of $U=8.0$ and $U=10.0$, which however brings a larger difference between $S$ and $E_{\rm corre}$ as listed in Table 1. On the other hand, when $U$ becomes larger, the correlation energy $E_{\rm corre}$ becomes smaller. Here it should be noted that the mean and standard errors of these quantities shown in Table 1 are derived by the reblocking analysis. Because the state of a simulation at one iteration depends upon the state at the previous iteration, the data points are not independent of each other. The reblocking analysis can remove the correlation between near iterations. Furthermore, when we estimate the projected energy, the correlation between the numerator and denominator must be considered by taking the covariance into account. The results of ground state projected energies $E_{\rm project}=H_{00}+E_{\rm corr}$ for different interaction strengths are listed in Table 2. In this table $E_{\rm FCIQMC}$ and $E_{\rm exact}$ are the expected values of the ground state energy by FCIQMC and by the Lanczos diagonalizations,[13] respectively. We find that the shift $S$ in FCIQMC can give the ground state energy of the Hubbard model with exact diagonalizations quality, but the differences between $E_{\rm FCIQMC}$ and $E_{\rm exact}$ are slightly larger when $U\geq 8.0$. This is because $\delta \tau $ is set much smaller to decrease the population of psips in FCIQMC simulation and the computational cost (CPU time and memory requirements), which leads to the population below the plateau slightly and has a divergent error in the estimate of the energy just as the analytical results for $N_{\rm e}$ in Ref. [14].

**Table 2.** Computed ground state energy from FCIQMC, compared with the Lanczos diagonalizations results.[13]
system	$U$	$S$	$E_{\rm FCIQMC}$	$E_{\rm exact}$
$4 \times 4$ ($7\uparrow 7\downarrow$)	1.0	$-21.39279(1)$	$-21.39283(2)$	$-$21.39283
$4 \times 4$ ($7\uparrow 7\downarrow$)	2.0	$-19.1709(1)$	$-19.17173(5)$	$-$19.17135
$4 \times 4$ ($7\uparrow 7\downarrow$)	3.0	$-17.2998(1)$	$-17.2995(5)$	$-$17.30003
$4 \times 4$ ($7\uparrow 7\downarrow$)	4.0	$-15.7447(2)$	$-15.7454(6)$	$-$15.74459
$4 \times 4$ ($7\uparrow 7\downarrow$)	6.0	$-13.42130(6)$	$-13.4213(3)$	$-$13.42123
$4 \times 4$ ($7\uparrow 7\downarrow$)	8.0	$-11.8689(2)$	$-11.8710(3)$	$-$11.86883
$4 \times 4$ ($7\uparrow 7\downarrow$)	10.0	$-10.8072(2)$	$-10.8058(1)$	$-$10.80701

For any Monte Carlo simulation the computational efficiency must be considered. In the following we will discuss the average time consumption per simulation cycle of FCIQMC for different interaction strengths and different numbers of CPU processors. As listed in Table 1, the plateau height increases with the interaction strength $U$. Generally speaking, the plateau height can roughly reflect how hard the calculation is for both memory and computer time, since each Monte Carlo iteration loops over all the particles. Figure 2(a) shows the computer time cost for different interaction strengths under the same condition that the number of computer processors is set to be 36 and the time step $\delta \tau=0.001$ for $U=4.0$. It is clear that the computer time increases with the interaction strength. In particular, when $U > 5.0$, the computer time increases very quickly. To decrease the computer time, we need much more computer resources (i.e., more processors), but with the increase of the processor number, the communication cost between different processors will increase (see Fig. 2(b)). The reason is that in the parallel FCIQMC software implementation (one part of HANDE) the Hilbert space is partitioned over the invoked processors, which can distribute the required memory across the processors effectively. When a psip on a parent determinant located in one part of the Hilbert space spawn a child on a determinant located in another part of the space, communication between processors will occur. With the increase of the communication rate, the time cost will be very large. In other words, it is the inter-processor communication that limits the parallel scale of FCIQMC.

Fig. 2. The average time consumption (in units of second) per simulation cycle of FCIQMC for different interaction strengths (a) and different numbers of computer processors with $U=4.0$ (b).

In summary, we have mainly calculated the ground state energies of the 2D Hubbard model near half-filling regime on a $4\times4$ square lattice system, and have studied the computational efficiency of the FCIQMC method. It is found that the results of energy shift $S$ with $1.0 \leq U \leq10.0$ and the results of the projected energy with $U\leq 6.0$ are in good agreement with those by the exact diagonalization. When $U$ is increased, the number of psips will increase in FCIQMC simulation, which implies that the computational cost (CPU time and memory requirements) will increase. To decrease computational cost, the number of psips can be decreased by setting smaller time step, which however brings the larger variance of the projected energy estimator for $U=8.0$ and $U=10.0$. Moreover, it is found that much more computer resources (i.e., more processors) are helpful to decrease the computer time, but with the increase of the processor number the communication cost between different processors will increase. When one uses FCIQMC to calculate the ground state properties of the 2D Hubbard model for the larger system size which is beyond the ability of the exact diagonalization method, one can use higher computing resource. How to reduce the required number of psips to achieve convergence and how to decrease the statistical uncertainty at the same time need to be studied. As far as we know, FCIQMC with the initiator approximation[15] can dramatically reduce the required number of psips for molecular systems, and the semi-stochastic method[16] can reduce the stochastic error for molecular and lattice systems. However, the former introduces a systematic error and the latter usually slows the simulation down. Whether to treat larger systems by combining the initiator approximation and semi-stochastic method will be studied in the next works. We thank Professor Youjin Deng and Professor Qianghua Wang for useful discussions about the Hubbard model, and also thank Dr. James Spencer for helpful discussions about the software HANDE. References Insights into the structure of many-electron wave functions of Mott-insulating antiferromagnets: The three-band Hubbard model in full configuration interaction quantum Monte CarloSolutions of the Two-Dimensional Hubbard Model: Benchmarks and Results from a Wide Range of Numerical AlgorithmsSpin-Ice State of the Quantum Heisenberg Antiferromagnet on the Pyrochlore LatticeQuantum Monte Carlo Calculations in Solids with Downfolded HamiltoniansBenchmark study of the two-dimensional Hubbard model with auxiliary-field quantum Monte Carlo methodSpontaneous magnetic flux in disordered mesoscopic rings with interacting electrons: Monte Carlo simulationsSpin density wave and ferromagnetism in a quasi-one-dimensional organic polymerDoping-driven orbital-selective Mott transition in multi-band Hubbard models with crystal field splittingiQIST: An open source continuous-time quantum Monte Carlo impurity solver toolkitCriticality and superfluidity in liquid

^{4} He

under nonequilibrium conditionsA random?walk simulation of the Schr?dinger equation: H ⁺ ₃Quantum chemistry by random walk. H ² P , H ⁺ ₃ D _{3 h} ¹ A ′ ₁ , H ₂ ³ Σ ⁺ _u , H ₄ ¹ Σ ⁺ _g , Be ¹ SQuantum chemistry by random walk: H4 squareQuantum Monte Carlo Method using Phase-Free Random Walks with Slater DeterminantsFermion Monte Carlo without fixed nodes: A game of life, death, and annihilation in Slater determinant spaceMetal-insulator transitionsA common thread: The pairing interaction for unconventional superconductorsThe sign problem and population dynamics in the full configuration interaction quantum Monte Carlo methodHole-hole effective interaction in the two-dimensional Hubbard modelUnderstanding and improving the efficiency of full configuration interaction quantum Monte CarloCommunications: Survival of the fittest: Accelerating convergence in full configuration-interaction quantum Monte CarloSemistochastic Projector Monte?Carlo MethodSemi-stochastic full configuration interaction quantum Monte Carlo: Developments and application

[1]	Schwarz L R, Booth G H and Alavi A 2015 Phys. Rev. B 91 045139
[2]	LeBlanc J P F, Antipov A E, Becca F et al 2015 Phys. Rev. X 5 041041
	Huang Y, Chen K, Deng Y J, Prokofev N and Svistunov B 2016 Phys. Rev. Lett. 116 177203
[3]	Ma F J, Purwanto W, Zhang S W et al 2015 Phys. Rev. Lett. 114 226401
	Qin M P, Shi H and Zhang S W 2016 Phys. Rev. B 94 085103
[4]	Liang S D, Wang Q H and Wang Z D 1997 Z. Phys. B: Condens. Matter 102 277
	Liang S D, Wang Q H and Wang Z D 1997 Z. Phys. B: Condens. Matter 104 27
[5]	Wang Y L, Huang L, Du L and Dai X 2016 Chin. Phys. B 25 037103
	Huang L, Wang Y L, Meng Z Y, Du L, Werner P and Dai X 2015 Comput. Phys. Commun. 195 140
[6]	Foulkes W M C, Mitas L, Needs R J and Rajagopal G 2001 Rev. Mod. Phys. 73 1
[7]	Anderson J B 1975 J. Chem. Phys. 63 1499
	Anderson J B 1976 J. Chem. Phys. 65 4121
	Anderson J B 1979 Int. J. Quantum Chem. 15 109
[8]	Zhang S W and Krakauer H 2003 Phys. Rev. Lett. 90 136401
[9]	Booth G H, Thom A J W and Alavi A 2009 J. Chem. Phys. 131 054106
[10]	Imada M, Fujimori A and Tokura Y 1998 Rev. Mod. Phys. 70 1039
[11]	Scalapino D J 2012 Rev. Mod. Phys. 84 1383
[12]	Spencer J S, Blunt N S and Foulkes W M C 2012 J. Chem. Phys. 136 054110
[13]	Fano G and Ortolani F 1990 Phys. Rev. B 42 6877
[14]	Vigor W A, Spencer J S, Bearpark M J and Thom A J W 2016 J. Chem. Phys. 144 094110
[15]	Cleland D, Booth G H and Alavi A 2010 J. Chem. Phys. 132 041103
[16]	Petruzielo F R, Holmes A A, Changlani H J, Nightingale M P and Foulkes W M C 2012 Phys. Rev. Lett. 109 230201
	Blunt N S, Smart S D, Kersten J A F, Spencer J S, Booth G H and Alavi A 2015 J. Chem. Phys. 142 184107