Twenty Years of Auxiliary-Field Quantum Monte Carlo in Quantum Chemistry: An Overview and Assessment on Main Group Chemistry and Bond-Breaking

In this work, we present an overview of the phaseless auxiliary-field quantum Monte Carlo (ph- AFQMC) approach from a computational quantum chemistry perspective, and present a numerical assessment of its performance on main group chemistry and bond-breaking problems with a total of 1004 relative energies. While our benchmark study is somewhat limited, we make recommendations for the use of ph-AFQMC for general main-group chemistry applications. For systems where single determinant wave functions are qualitatively accurate, we expect the accuracy of ph-AFQMC in conjunction with a single determinant trial wave function to be between that of coupled-cluster with singles and doubles (CCSD) and CCSD with perturbative triples (CCSD(T)). For these applications, ph-AFQMC should be a method of choice when canonical CCSD(T) is too expensive to run. For systems where multi-reference (MR) wave functions are needed for qualitative accuracy, ph-AFQMC is far more accurate than MR perturbation theory methods and competitive with MR configuration interaction (MRCI) methods. Due to the computational efficiency of ph-AFQMC compared to MRCI, we recommended ph-AFQMC as a method of choice for handling dynamic correlation in MR problems. We conclude with a discussion of important directions for future development of the ph-AFQMC approach.


I. INTRODUCTION
As computational power grows and theoretical algorithms have improved, the importance of the role played by electronic structure theory in chemistry, physics and biology has increased dramatically. 1 For a number of important applications, the predictive power of quantum chemistry has reached the level of "chemical accuracy," namely relative errors of one kcal/mol. 2 However the most accurate methods are generally the most computationally expensive, and for systems with many electrons, often one must resort to approaches where convergence to predictive accuracy is far from guaranteed. When "goldstandard" wave function-based approaches are impractical, Kohn-Sham density functional theory (DFT) has been the method of choice, due to its combination of accuracy and mild scaling with system size. However, there are many variants of DFT which depend on the choice of functional and it can be difficult a priori to select the most accurate density functional for a given application. In addition, DFT generally suffers from self-interaction errors and has difficulty in the treatment of the effects of strong static correlation prevalent in the breaking of chemical bonds and in the electronic structure of transition metal-containing systems. 3 Indeed, even more expensive wave function approaches such as coupled-cluster theory 4 can have difficulty accurately describing strongly correlated cases. Clearly there is still room for the development of novel methods in quantum chemistry which have a high ratio of accuracy to expense and may be reliably applied to challenging chemical problems.
Since the viability of an electronic structure method depends critically on the trade-off between accuracy and scaling with respect to system size, the development of new techniques that push the envelope of accuracy and scalability is the primary goal of the field. Adoption of a method by the community cannot occur before careful and extensive benchmark studies are carried out to establish the domain of validity as well as the benefits and drawbacks of the method in question. It is our aim here to review the foundations of one particularly promising electronic structure framework, namely the auxiliaryfield quantum Monte Carlo (AFQMC) approach, and to present, for the first time, the performance of this approach on a large-scale chemical benchmark set.
Unlike DFT or conventional wave function methods, AFQMC employs a statistical as opposed to deterministic route to calculating electronic structure. This necessitates consideration of an entirely different approach to the extraction of ground state energies, and requires the consideration of the role played by statistical fluctuations and error in the solution of the Schrödinger equation. The use and development of AFQMC has notably accelerated in recent years due to new algorithmic advances 5-9 and the publication of promising results on non-trivial chemical problems. [10][11][12][13][14][15][16][17] Unlike standard DFT approaches, AFQMC is free of self-interaction error 16 (as it is a wavefunction method) and is often capable of describing strong correlation effects, 9,[18][19][20][21][22][23][24][25][26][27][28][29][30][31] although our focus is not limited to such problems in this work.
AFQMC has its roots in the determinant QMC (DQMC) method proposed by Blankenbecler, Scalapino, and Sugar. 32 DQMC is an exact, unbiased technique, initially developed to describe interacting bosons and fermions at finite temperatures in lattice field theory models of high energy physics. Since its invention, DQMC has been applied to numerous condensed matter problems, most notably the Hubbard model. 33 More recently this approach has been applied to more realistic chemical systems with long-ranged interactions. [34][35][36][37] As is the case for DQMC, the scalability of AFQMC is limited by the fermionic sign problem, rendering exact, brute force applications to large chemical systems com-putationally infeasible. 38 To be more precise, the sign problem manifests in the following way: Per statistical sample, the complexity of the calculation scales polynomially with system size. However, the sample complexity (i.e., the number of statistical samples required for a fixed statistical error) scales exponentially with system size, rendering practical calculations impossible.
To remedy this fundamental barrier in DQMC, Fahy and Hamann proposed the approximate positive projection method, 39 which imposes a constraint to control the sign problem, resembling the fixed-node approximations used in Green's function Monte Carlo (GFMC) 40,41 and diffusion Monte Carlo (DMC). 42 Adopting this approximation for open-ended random walks, Zhang, Carlson, and Gubernatis 43 proposed the constrained path approximation, and this method will be referred to as cp-AFQMC. cp-AFQMC has been successfully applied to non-trivial electronic lattice problems, [44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61] bosonic lattice problems, 62,63 and mixed fermions-bosons systems. 64,65 cp-AFQMC has been shown to be one of the more accurate and scalable many-body methods in numerous benchmark studies, 49,57,58 although it is important to bear in mind that like all constrained quantum Monte Carlo methods, cp-AFQMC is approximate and its domain of validity has not been fully mapped out.
The generalization of cp-AFQMC to chemical systems with long-range interactions yields the phaseless AFQMC (ph-AFQMC) which is the focus of this work. 66 Since its invention by Zhang and Krakauer in 2003, 66 ph-AFQMC has gained popularity in the quantum chemistry community due to its relative computational affordability and its accuracy, as demonstrated in recent calculations including those performed on small but challenging molecular systems, 9,[18][19][20][21][22][23][24][25][26][27][28][29][30][31] simple transition metal complexes, [10][11][12][13][14][15][16][17] and solids. 5,67-80 ph-AFQMC is highly flexible in that the accuracy of the approach may be systematically improved by the use of increasingly sophisticated trial wave functions. In addition, ph-AFQMC has become more accessible to users via the development of open-source codes. [81][82][83] Despite its growing usage in the chemistry community, it is not clear how well ph-AFQMC performs for a wide range of main group chemistry applications. While it has been speculated that ph-AFQMC with the simplest form of a single-reference trial wave function has an accuracy on par with coupled-cluster with singles, doubles, and perturbative triples (CCSD(T)), 18,31,73,84 there has been no extensive body of work that confirms this expectation. While ph-AFQMC has been compared to other commonly used quantum chemistry methods such as CCSD(T) on simple bond breaking problems such as those that occur in H 2 O, 18 BH, 21 N 2 , 21 F 2 , 23 C 2 , 22 Cr 2 , 11 and Mo 2 , 12 studies comparing ph-AFQMC with both single-reference CCSD(T) and other multireference methods have been relatively scarce. Given the accuracy of ph-AFQMC for dynamic correlation problems, 15,31,69,75 it is crucial to compare its accuracy with other multi-reference dynamic correlation methods, such as perturbation theory and configuration interaction methods.
The goal of this manuscript is both to present a pedagogical review of ph-AFQMC from the quantum chemistry perspective as well as to investigate the performance of ph-AFQMC on a large chemical benchmark series, namely a thermochemistry benchmark set (W4-11), 85 a non-covalent interaction benchmark set (A24), 86 and model molecular bond dissociation problems (H 4 and N 2 ). [87][88][89][90][91][92][93][94][95] From this, we hope to establish fair expectations for ph-AFQMC for main group chemistry and simple strong correlation problems associated with the breaking of chemical bonds, and to further encourage the use and development of ph-AFQMC in quantum chemistry in the future. We note that there exist other pedagogical reviews written from a somewhat different perspective [96][97][98] which the interested reader is encouraged to consult.
It should be noted that the applications considered here, while chemically important, are somewhat uncommon for ph-AFQMC. In particular, the benchmark sets (W4-11 85 and A24 86 ) considered here are predominantly single reference problems. These do not involve strong electron correlation for which ph-AFQMC has been frequently used. Indeed, these benchmark problems are already known to be well-characterized by coupled-cluster techniques. Furthermore, we mostly consider the least accurate and the most efficient variant of ph-AFQMC. Namely we employ a simple single-determinant trial wave function, although we do discuss selected examples where the use of multi-determinant trials can greatly improve accuracy. These choices are purposeful ones: we aim to assess the breadth and baseline accuracy of the approach in situations where other quantum chemistry methods such as CCSD(T) are often well-suited. The main points that we will deliver in this manuscript are as follows: 1. For main group chemistry applications, with O(N 3 ) -O(N 4 ) scaling per sample, ph-AFQMC with a single-reference trial wave function is more accurate than CCSD (O(N 6 )) and is competitive with, albeit somewhat less accurate, than CCSD(T) (O(N 7 )) and state-of-the-art density functionals.
2. For bond breaking, while maintaining O(N 3 ) -O(N 4 ) scaling per sample, the accuracy of ph-AFQMC with a multi-reference trial wave function surpasses that of low-order multi-reference perturbation theory (O(N 5 ) -O(N 6 )) and is comparable in accuracy with multi-reference configuration interaction methods (O(N 6 )).
3. Given the relatively short history of ph-AFQMC in quantum chemistry, many aspects of the algorithm deserve further investigation and there is much room for improvements to the current implementation of the approach. We will conclude this work with a discussion of these opportunities.
This manuscript is organized as follows: In Section II we present an overview of the AFQMC method and its phaseless variant (ph-AFQMC). A discussion of trial wave functions, the calculation of observables, sizeconsistency and computational cost is contained in this section. In Section III we present benchmark results with a comparison with other approaches. Section IV presents a discussion of lessons learned from our calculations and prospects for future development. We conclude in Section V. The Appendices contain details omitted from the main text.

II. OVERVIEW OF PHASELESS AUXILIARY-FIELD QUANTUM MONTE CARLO
A. Ground-state calculations We refer the interested readers to Appendix A0 A1 for a more detailed formal exposition of AFQMC. Here, we focus on a short summary of the methodology. We assume the spin-orbital notation throughout this paper.
AFQMC is what is called a "projector" QMC algorithm because it projects towards the ground state from an initial wave function which has a non-zero overlap with the true, exact ground state of the system. Formally, the ground-state wave function |Ψ 0 of a HamiltonianĤ is found by via imaginary time propagation whereĤ is an ab initio Hamiltonian, τ is the imaginary time, and |Φ 0 is an initial wave function that is not orthogonal to |Ψ 0 . The form of |Ψ 0 in Eq. (1) resembles the CCSD wave function with generalized singles and doubles, 99,100 but in AFQMC there is no cluster amplitude to be determined, and formally Eq. (1) is exact. Discretizing τ and applying the Trotter decomposition 101 followed by the Hubbard-Stratonovich (HS) 102 transformation to the propagator e −∆τĤ , one obtains a high-dimensional integral formula for the imaginary time evolution operator over a short time interval where ∆τ is the infinitesimal time step, p(x) is the standard normal distribution andB(x) is an effective onebody propagator that is obtained from one-body operators coupled to x, a vector of N γ classical auxiliary fields. Since only one-body operators now appear, this representation of the short-time propagator is often viewed as the exact mapping of an interacting many-body system to an ensemble of non-interacting systems, each coupled to a classical, fluctuating external potential. The many-body correlation beyond mean-field theory is recovered upon the integration over the auxiliary fields. The phaseless AFQMC (ph-AFQMC) 66 approach utilizes importance sampling based on a trial wave function |Ψ T during the random walk process. To carry this out, one writes the global wave function at time τ as a weighted statistical sum over N walkers "walkers", Given this importance sampling, the walker weight update rule follows where the overlap ratio S is given by Since this weight update rule cannot guarantee that the walker weights are real and positive, {w i }, Eq. (5) leads to the fermionic phase (sign) problem with a diverging variance for any observable when calculated by the protocol outlined below. 38 To control the phase problem, ph-AFQMC imposes a constraint to ensure the positivity of the weights throughout the imaginary time propagation. The constraint that defines the standard ph-AFQMC approach is given by a modified overlap ratio, S ph , where the phase θ i (τ ) is given by This modified overlap ratio is used in ph-AFQMC to update the weights. It should be emphasized that this constrained random walk introduces systematic biases which induce deviations in the values of observables (such at the ground state energy) compared to the exact values. It should be noted that this bias disappears in the limit of |Ψ T becoming an exact eigenstate ofĤ. This fact introduces a means to improve the accuracy of ph-AFQMC by increasing the sophistication of the trial function |Ψ T . We will return to this point later in this work. A different way to view the constraint is that one imposes a gauge boundary condition (i.e., a global phase) for the wave function sampled through the imaginary time evolution. Within this framework, the global energy estimate at a given imaginary time τ is given by a weighted statistical average, where the local estimate for observableÔ of the i-th walker is given by IfÔ =Ĥ, these estimates are global energy and local energy estimates, respectively. Given the positivity of {w i }, the variance of this energy estimate grows linearly with system size, ensuring the polynomial-scaling sample complexity of the overall algorithm. It is important to remark here that the resulting ph-AFQMC energy computed via Eq. (9) is not guaranteed to be variational. 103

B. Trial wave functions
The accuracy of ph-AFQMC heavily depends on the quality of the trial wave function employed in the calculation. Conceptually, the role played by the trial wave function is very different from that of the reference wave function used in conventional quantum chemistry methods such as CCSD(T). Despite this difference, it is convenient to think of ph-AFQMC as adding correlation energy on top of an a priori chosen trial wave function. Here, we summarize existing strategies for generating these ph-AFQMC trial states.
Single determinant (SD) trial wave functions are the most widely used and they offer the most affordable variant of the ph-AFQMC algorithm. There are multiple approaches to obtain single determinant trial wave functions. One way is to employ the lowest energy spinunrestricted or spin-generalized Hartree-Fock (HF) wave function. 16,30 This approach can be well-suited for describing bond dissociation as shown in Section III C. While this approach is completely parameter-free, it often runs into issues associated with artificial symmetry breaking. 16,104 In particular, HF wave functions can exhibit an unphysical breaking of symmetry that makes post-HF calculations behave erratically in energetics and properties. [105][106][107][108][109][110][111][112] While ph-AFQMC can correct artificial symmetry breaking to some extent through imaginary time evolution, there are cases where the lowest energy HF solution is clearly not the best SD trial wave function. 16 Alternatively, one can obtain a single Slater determinant from methods that include electron correlation effects. For this, DFT 10,66 and regularized orbital optimized Møller-Plesset perturbation theory (OOMP2) have been used. 16,113 More broadly, one can employ any approximate Brueckner orbital wave functions 114,115 for this purpose. We also note that the recently-defined selfconsistent trial wave function method 98,116 also fits into the category of an SD trial wave function approach.
More elaborate trial wave functions may be used within ph-AFQMC, which confer greater accuracy for a greater computational expense. Most commonly, a linear combination of SDs, often referred to as multi-Slater determinants or MSDs for short, are employed. 9,21,23,30 The MSD trial can be obtained from a truncated configuration interaction (CI) expansion, complete active space CI (CASCI), or from selected CI methods. 117 While these trial wave functions can yield excellent ph-AFQMC energies, they are brute-force in nature and the cost of obtaining these trial wave functions generally scales exponentially with system size, which can limit the applicability of their use within ph-AFQMC.
Lastly, there are non-linear trial wave functions which may be used.
These include Jastrow factors, 51,118 coupled-cluster wave functions, 4,37,119 perfect-pairing wave functions, 23,119,120 and transcorrelated wave functions. 121 The cost of obtaining these wave functions in all such cases is polynomial-scaling, which makes them appealing for use as trials. Furthermore, these wave functions approximate the exact ground state much more accurately than do SD trial functions as they all include electron correlation inherently missing from the SD description. Unfortunately, using these non-linear wave function as a trial wave functions without any approximations is currently extremely difficult on both classical and quantum computing platforms. 119 The efficient and accurate use of such trial states within ph-AFQMC will require further work to become viable.

C. Physical properties extracted from ph-AFQMC
In practical applications, one may want to evaluate observables that do not commute withĤ. Most commonly, these observables are associated directly with reduced density matrices (RDMs) or some elements of them. The computation of such observables introduces additional challenges within the ph-AFQMC framework, because the mixed estimator in Eq. (9) is no longer unbiased. For such quantities, we need to use a different estimator known as the "pure" estimator. This problem is analogous to the issues associated with coupled-cluster expectation values 122 which are typically now handled by the coupled-cluster Lagrangian formalism. 123 A computationally simple approach to this problem is to use the approximate variational estimator which is frequently used in diffusion Monte Carlo. 124,125 Recently, Eq. (11) was used within the ph-AFQMC framework to estimate the dipole moments of simple molecules. 9 Here, it was found that the quality of the trial wave function is critically important for accurately approximating the pure estimates. For this quantity, using SD trial wave functions yielded very inaccurate results for the systems considered in Ref. 9.
Another approach that is more commonly used in ph-AFQMC calculations is the back-propagation algorithm. 25,62,126 Here, one propagates the bra state ( Ψ T |) in the mixed estimator (Eq. (A10)) backward in time such that the resulting estimator is symmetric with respect to both the bra and ket, ultimately resembling the pure (or variational) estimator. Formally, this entails computing observables using the following estimate: where the back-propagation time, κ, while formally taken to ∞, is in practice fixed to a long enough finite time length. This approach is computationally efficient, but its accuracy for at least some ab initio systems has been shown to be rather poor. 25,80 There are additional algorithmic ways to reduce the back-propagation bias, such as partially restoring the phase and cosine factors along the back propagation portion of the path. 25,98 However, even with these considerations, the resulting 1-RDMs were found to be inaccurate, at least in some systems. 80 Devising improved algorithms for accurately computing pure estimates for observables without greatly increasing the computational cost of the calculations is a worthy goal for future research.

D. Size-consistency
Size-consistency is a property of a wave function for isolated systems A and B that guarantees the product separability of a supersystem wave function (|Ψ AB = |Ψ A |Ψ B ) as well as the additive separability of energy (E AB = E A +E B ). 127 This property has important implications concerning the applicability of a given method to large systems and is therefore considered to be an important formal property in method development. The sizeconsistency of ph-AFQMC was first examined in Ref. 73, but some subtle issues were not fully discussed in that work. Here, we will provide the first rigorous analysis of size-consistency within ph-AFQMC.
We will assume that the trial wave function is product separable, and that for isolated systems A and B, the total Hamiltonian separates intoĤ A andĤ B which commute with each other. Given these conditions, the propagator is product separable, (13) which leads to the product separability ofB,B AB = B ABB .
Provided that the walker wave function is product separable, it can be shown that the overlap ratio in Eq. (6) of the total system, S AB , is product separable into monomer overlap ratios, S A and S B The product separability of the walker wave function can be satisfied as long as we start from a product separable wave function since the propagator itself is product separable. With this overlap ratio, walker weights are also product separable as w AB = w A w B . Since the local energy in Eq. (10) is additively separable (E AB However, the above analysis does not apply to the modified overlap ratio for arbitrary imaginary time steps in Eq. (7) used in the constrained ph-AFQMC formalism. To see this, noting the cosine projection which defines the constraint for the ph-AFQMC framework, we have where the cosine factor is not product separable, leading to size-inconsistency in the overall approach if the time step is not taken to zero. In other words, the walker weights of the total system are not product separable (w AB = w A w B ). The magnitude of the sizeinconsistency error is proportional to the magnitude of the sine terms in Eq. (16), and can be practically small. Importantly, if the time step ∆t is small, then θ A and θ B are small. Therefore, we expect the size-inconsistency error to vanish in the limit of ∆t → 0. We present numerical results to support this in Appendix A0 A2. In this sense, the size-inconsistency error is a part of the time step error with an O(∆t) error scaling as opposed to O(∆t 2 ) of the Trotter error. Like other time step errors, the above analysis suggests that size-inconsistency error can be controlled within ph-AFQMC by taking the limit of ∆t → 0. While we examined the size-consistency of ph-AFQMC assuming localized orbitals and productseparability of the trial wavefunction, it is possible that these assumptions may not be necessary to examine the size-consistency more generally. In summary, we have shown that ph-AFQMC is strictly size-consistent in the limit of ∆t → 0. The sizeconsistency is critically important to reliably apply ph-AFQMC to large systems and solid state problems, similarly to more traditional quantum chemistry methods such as coupled-cluster theory.

E. Computational cost
The computational expense and scaling with system size of ph-AFQMC is an important factor when considering suitable applications for this approach. The integral transform often necessary for ph-AFQMC scales as O(N 4 ). There are three additional considerations that need to be considered within ph-AFQMC to determine the scaling behavior and how it depends on the choice of a trial wave function. For the following discussion, we will assume that the walker wave functions are SDs. The first consideration concerns the propagation step. The central quantity to be computed is the modified overlap  (10))) for singledeterminant (SD), multi Slater determinant (MSD), and nonlinear trial wave functions. N denotes system size and Nc is the number of determinants in an MSD trial wave function. The cubic-scaling for the local energy evaluation in the case of SD follows from arguments presented in Refs. 5-8. ratio in Eq. (7). Specifically, the overlap between the trial and walker wave functions, Ψ T |ψ i (τ ) , must be efficiently evaluated if the method is to be computationally viable. In addition to the overlap, one needs to evaluate the one-body Green's function for the force bias evaluation (see Appendix A0 A1 for more details). The final consideration relates to cost of the local energy evaluation in Eq. (10) assuming that the ground state energy is the quantity of interest. In Table I we summarize the cost of these three parts for different types of trial wave functions discussed in Section II B.
For the local energy evaluation there has been a variety of recent algorithmic improvements. For example, the best algorithms now available reduce the standard quartic-scaling algorithm of the SD local energy evaluation to a cubic-scaling algorithm. This can be achieved by using double factorization, 6 tensor hypercontraction, 5 stochastic resolution-of-the-identity, 7 or via the use of localized orbitals. 8 Accelerating the local energy evaluation with multi-SD trials has not been explored as extensively as that of the simpler SD case, but recent explorations with localized orbitals 8 and with generalized Wick's theorem 9,37 are encouraging. Future work should be aimed at lowering the complexity of the overlap and the force bias evaluation since these computations often form the bottleneck for medium-sized molecules.
If one is to compare the cost of standard deterministic quantum chemistry methods to the cost of ph-AFQMC, a subtlety arises due to the statistical nature of ph-AFQMC. For bulk systems when energy per particle or other size-intensive quantities are relevant, our cost analysis given in Table Table I is sufficient, assuming the number of samples required for desired precision does not grow with system size. For finite molecular systems, it is sensible to estimate the cost of ph-AFQMC for a fixed statistical error as the system size, N , increases. Assuming that the auto-correlation time in the Markov chain does not grow with N and that the standard error in the energy estimate grows linearly with N , one crudely requires O(N 2 ) statistical samples to maintain a fixed statistical error. Therefore, we find a O(N 2 ) multiplicative factor in addition to the cost in Table I. While this is a correct formal asymptotic scaling, it is possible that in practice one may not experience O(N 2 ) sample complexity if the trial wave function is accurate such that the statistical fluctuations are suppressed below a desired error threshold for the range of system sizes under consideration.

III. PERFORMANCE OF PH-AFQMC FOR THERMOCHEMISTRY AND THE TREATMENT OF NON-COVALENT INTERACTIONS AND BOND DISSOCIATION
We will now assess the accuracy of ph-AFQMC using well-known thermochemistry and non-covalent interaction benchmark sets, along with several simple bond dissociation examples. For simplicity, we will refer to ph-AFQMC as AFQMC in this section. Computational details are available in Appendix A0 A3. Note that nearly all of the AFQMC calculations presented here are new and have not been previously published.
We note that not every data point in the W4-11 set is a simple single-reference (SR) problem. The TAE140 subset contains 16 multi-reference (MR) data points determined by the %TAE e [(T)] diagnostic. 85 Similarly, 202 energies in the HAT707 set and 16 energies in the BDE99 set have also been deemed MR data points. 2 We will discuss the performance of ph-AFQMC on 745 SR data points seperately from the 234 MR data points, as is often done for density functionals and other quantum chemistry methods. 2

Single atom total energies
The simplest systems in the W4-11 set are single atoms (Be, B, C, N, O, F, Al, Si, P, S, Cl). These are systems are small enough that brute-force methods such as SHCI can be reliably performed to obtain exact total energies. We first look at the performance of CCSD(T) and AFQMC for these in the aug-cc-pVTZ basis set as shown in as a reference state for CCSD(T) and as a trial state for AFQMC. The conclusions are unchanged if one uses a spin-unrestricted Hartree-Fock (UHF) trial instead. In Fig. 1, CCSD(T) achieves chemical accuracy (1 kcal/mol) for all the atoms in W4-11. However, AFQMC exhibits large errors above 2 kcal/mol for numerous atoms. There is an interesting trend that the AFQMC error becomes an order of magnitude smaller for atoms whose ground state angular momentum is of S symmetry. For all the other cases (ground state angular momentum, P), AFQMC errors are greater than 2 kcal/mol. We conducted a preliminary investigation of this trend using various different SD and MSD trials with the hypothesis that these large errors can be attributed to simple symmetry constraints, as seen in the Hubbard model. 48,50 Unfortunately, while one can converge all of the atomic energies to chemical accuracy with a sophisticated multi-SD trial, we were unable to find a simple, compact trial wave function that removes this bias. For instance, we investigated C in aug-cc-pVTZ further with an exact trial wavefunction generated by SHCI. To achieve chemical accuracy, we needed more than 200 determinants in the trial wavefunction which has more than 0.99 overlap with the exact ground state.
Very recent work suggests that using SD trial wave functions with a partially relaxed constraint may reduce the bias significantly although the accuracy of such a constraint on large systems remains unclear. 129 Alternatively, one could even perform free-projection 37 or release constraint 48 calculations with SD trial wave functions on these small atoms, but this is not a viable option for larger systems. We leave further investigation of these simple systems for future studies, and focus on the performance of AFQMC for W4-11 where we use CCSD(T) atomic energies but AFQMC molecular energies. Lastly, we note that due to the frozen core approximation, Be contains only two electrons that need to be correlated. Therefore, CCSD is exact for this system, as shown in Fig. 1. However, this exactness does not apply to AFQMC. We have performed CCSD(T) and AFQMC using both SD RHF and UHF states as references or trial functions, respectively. We report the raw energies of each case in the Zenodo repository. 130 For AFQMC, it was found that the use of UHF trials leads to more accurate results, while for CCSD(T) using an RHF references was statistically superior. Hence, we will compare AFQMC with UHF trials and CCSD(T) with RHF reference for all singlereference calculations.
In Fig. 2, we observe systematic improvements for both the root-mean-square deviation (RMSD) and meansigned error (MSE) as we increase the sophistication of the correlation treatment, with second-order Møller-Plesset perturbation theory (MP2) less accurate than CCSD which is less accurate than CCSD(T). However, both MP2 (9.01 kcal/mol) and CCSD (5.40 kcal/mol) do not achieve chemical accuracy in terms of RMSD for any of the W4-11 subsets. CCSD(T) achieves chemical accuracy in all subsets of W4-11 (total RMSD = 0.63 kcal/mol). We compare our results against a combinatorially optimized density functional, namely ωB97M-V, 131 which was found to be the best functional out of 200 functionals examined in Ref. 2. While there are now even more accurate functionals (e.g. double hybrid functionals 132,133 ), we believe that ωB97M-V serves as an example of an accurate functional for the types of systems we investigate here. Data for other functionals are available in Ref. 2. ωB97M-V provides chemical accuracy only for the SN13 subset and within the complete W4-11 set its RMSD is 2.52 kcal/mol. AFQMC performs quite similarly to ωB97M-V in that it achieves chemical accuracy for the SN13 subset but not for any other subsets. Overall, AFQMC's RMSD is 3.26 kcal/mol which is slightly worse than ωB97M-V. However, we note that this is still an improvement over CCSD (5.40 kcal/mol) by a sizable margin. As the data points in TAE140 completely determine relative energies in all the other subsets, we see qualitatively similar statistical results between the TAE140 subset and the entire W4-11 data set. For 234 MR data points, we performed the same analysis as shown in Fig. 3. With the exception of the error of CCSD for the TAE140 MR set, the relative performance between different methods over MR data points is the same as that over each of the SR data sets and within W4-11 overall. Similarly to the SR case shown in Fig. 2, MP2 (10.51 kcal/mol) and CCSD (8.59 kcal/mol) are found to be significantly less accurate than CCSD(T) overall. CCSD(T) also does not perform as well for these data points, yielding an RMSD of 1.44 kcal/mol. ωB97M-V has an RMSD of 4.27 kcal/mol overall for the MR W4-11 data set. AFQMC is more accurate than all but CCSD(T) for this subset, with an RMSD of 2.88 kcal/mol. Although AFQMC does not perform as well as CCSD(T) for these data points, it is significantly more accurate than the other approaches considered here.
The good performance of CCSD(T) for these MR data points may be understood by the fact that many of these MR calculations are not necessarily strongly correlated ones. There are multiple ways to diagnose MR character, with different metrics providing distinct classifications. For example, if we use regularized orbitaloptimized MP2 (κ-OOMP2) 113 and inspect the underlying spin-symmetry breaking of the solution, only 7 out of the 16 MR data points in the TAE140 set exhibit spinsymmetry breaking. Based on this, one would conclude that only 7 data points in TAE140 should be considered to carry MR character. We will discuss bond breaking examples in the later sections which will unambiguously fall into the strongly correlated category for stretched geometries. To emphasize the flexibility of AFQMC, we further study C 2 in aug-cc-pVTZ which is one of the representative MR examples in the MR16 subset. We generate a large MSD trial using SHCI for an active space with 8electron and 90-orbital. We then systematically converge AFQMC energies using that trial to the exact answer obtained by SHCI. As shown in Fig. 4, MSD-AFQMC quickly improves its accuracy as one adds more determinants to the trial and becomes chemically accurate with 10 4 determinants or so. Compared to more conventional quantum chemistry methods such as coupled-cluster theory, the generalization of AFQMC to MSD trials is rather straightforward and it can often be used to improve the results significantly beyond AFQMC with an SD trial. One can also try to find other MSD trials generated from a smaller active space to reduce the number of determinants needed to reach chemical accuracy. 22 We will see the power of MSD-AFQMC more later in Section III C 2.

MP2
CCSD In summary, based on this thermochemistry benchmark study, we conclude that the performance of AFQMC with a simple SD trial is expected to provide accuracy that lies between CCSD and CCSD(T). It is also quite competitive with a state-of-the-art density functional, ωB97M-V, which is a priori expected to be accurate for the class of systems studied here. These conclusions are summarized in Fig. 5. We thus recommend AFQMC with a simple SD trial for calculations where DFT is expected to struggle (either due to self-interaction error or due to strong correlation) and CCSD(T) is too expensive.

B. Non-covalent interaction benchmark (A24)
A24 is a high-quality non-covalent interaction benchmark set composed of 24 small, non-covalently bound complexes. 86 All molecules in the set should be well described by RHF wave functions. For this set we employ RHF wave functions for both CC and AFQMC methods. We also used the counterpoise correction to eliminate the basis set superposition error. The energy scale in this set is quite small, with interaction energies ranging from 1.115 kcal/mol to -6.493 kcal/mol. In Fig. 6, we present the RMSD values of MP2, CCSD, CCSD(T), ωB97M-V, and AFQMC for this set. MP2 has the largest RMSD (1.12 kcal/mol), while CCSD and CCSD(T) are nearly exact (RMSD = 0.05 kcal/mol). ωB97M-V is also very accurate with an RMSD of 0.09 kcal/mol. AFQMC exhibits more sizable errors than do the CC methods, with an RMSD of 0.32(7) kcal/mol. This is about a factor of 4 larger error than ωB97M-V, although still within chemical accuracy.
For relatively simple non-covalently bound systems, combinatorially optimized density functionals like ωB97M-V work well, approaching the accuracy of CCSD(T). 131 Therefore, we recommend the use of AFQMC for non-covalent interaction energies only if the system under study exhibits additional challenges for DFT, such as those arise in cases with open-shell electrons or self-interaction error. At the same time, the accuracy for AFQMC in such cases in comparison to CCSD(T) is still under-explored. Given the small energy scale and the small number of data points in A24, a more thorough investigation of of non-covalent interactions with AFQMC is highly desirable. The recent report 134 by Tkatchenko and co-workers of significant deviations between diffusion Monte Carlo and CCSD(T) for large, non-covalently bound complexes thus provides a worthy target for a future investigation performed with AFQMC.

C. Breaking chemical bonds
Bond breaking in simple molecular systems has been a test bed for the electronic structure treatment of strong correlation. This is because one can simply tune the strength of the electron correlation (i.e., the degree to which an SR wave function fails qualitatively) by changing the bond distance. This fact was recognized early in the development of AFQMC, as evidenced by the work of Zhang and co-workers. 11,12,[21][22][23] In this section, we add more data to this body of work by considering H 4 in STO-3G and cc-pVQZ bases and N 2 (all-electron) in STO-3G and cc-pVTZ bases. STO-3G is a minimal basis set and is known to exaggerate strong correlation effects. 135 Therefore, bond breaking in the STO-3G basis set has been a particularly popular test for strong correlation methods. [87][88][89][90][91][92][93][94][95] Larger basis set examples will exhibit a mixture of strong and weak correlation effects, which is a common challenge encountered in realistic problems. where we vary the distance between two stretched H 2 units, with each unit with a fixed bond distance of 1.23 A. This value is selected because for R H−H = 1.23Å , two RHF determinants become quasi-degenerate, maximizing the strong correlation aspect of this model. Thus RHF is qualitatively incorrect and UHF cannot fully remedy this problem. A natural active space for this problem is (4e,4o) (i.e., 4-electrons in 4-orbitals), and this is what we use to generate complete active space self-consistent field (CASSCF) states for our subsequent dynamic correlation treatment.
In the case of N 2 , we vary the bond distance, R N−N , between the two nitrogen atoms. As the bond is stretched, N 2 exhibits strong correlation effects, and ultimately dissociates to two quartet nitrogen atoms. RHF cannot correctly dissociate N 2 whereas UHF dissociates this molecule correctly. For large bond distances, UHF is therefore expected to be qualitatively correct, while RHF qualitatively fails. Nonetheless, for intermediate bond distances, UHF is incapable of describing spin-recoupling and will thus lead to quantitatively inaccurate results in this regime. For CASSCF, we use a minimal active space of (6e,6o) (i.e., 6-electrons in 6-orbitals).
To add dynamic correlation on top of the CASSCF calculation, we consider some of the more widely used multi-reference methods in addition to AFQMC. These methods are summarized in Table II along with their acronyms. For CASPT2 we did not use any shifts. [149][150][151][152] NEVPT2 is performed via the partially contracted algorithm, 139 whereas the rest of the deterministic MR methods employed fully internally contracted versions. These methods represent commonly used dynamic correlation approaches for MR systems.

Assessment of single-reference methods
We first discuss the performance of single-reference methods for H 4 . In Fig. 7 (a), we present the error for each single-reference method for the STO-3G basis set. We examine both RHF and UHF states for CCSD(T) and AFQMC for use as reference or trial functions. These different calculations are referred to as RCCSD(T)/UCCSD(T) and RAFQMC/UAFQMC, respectively.
The sharp derivative discontinuity at R H−H = 1.23 A in RCCSD(T) and RAFQMC is due to the two RHF solutions crossing at R H−H = 1.23Å. As noted previously in literature, 94  Error in the potential energy curve of H4 as a function of RH−H for methods based on HF states (a) using a STO-3G basis set and (b) using a cc-pVQZ basis set. Note that there are more data points than markers shown. See Fig. A2 for qualitative shapes of these potential energy curves in absolute energy scale. the exact energy, but the magnitude of the error is about nine to ten times larger than that of RCCSD(T). UCCSD(T) produces variational energies at every distance considered here. UAFQMC performs far better than RAFQMC at every bond distance and also outperforms UCCSD(T) up to R H−H = 1.64Å. This is highly encouraging because UAFQMC is very accurate near the distances where strong correlation is most pronounced. Nonetheless, UAFQMC does not dissociate this molecule into two independent H 2 correctly. Both RCCSD(T) and UCCSD(T) correctly describe dissociation, but neither does RAFQMC nor UAFQMC. In Section III A 1, we discussed the exactness of CCSD(T) and the non-exactness of AFQMC for two-electron systems. The H 4 → 2 H 2 dissociation problem is thus another good illustration of this point.
We observe similar behavior of all four methods in the larger basis set, cc-pVQZ, as illustrated in Fig. 7 (b). Since strong correlation is less pronounced in a larger basis set, the magnitude of the error produced by all methods is significantly smaller compared to that exhibited in the STO-3G basis. Nonetheless, qualitative failures of RCCSD(T), such as non-variationality and a derivative discontinuity at the square geometry are still observed in this basis set. As is the case for the STO-3G basis, RAFQMC exhibits a significantly larger error than RCCSD(T) does in this basis set. Finally, we observe quite similar levels of quantitative accuracy when comparing UCCSD(T) to UAFQMC.
We continue the discussion of the performance of single-reference methods for N 2 as shown in Fig. 8   rable between these two bases sets. While there are multiple RHF solutions for N 2 , we focused on the one that does not break spatial symmetry. RCCSD(T) exhibits incorrect turnover behavior and non-variationality at far distances, whereas UCCSD(T) has a large energy error in the spin-recoupling regime. UCCSD(T) ultimately dissociates correctly into two N atoms, but its error up to a bond distance R N−N = 2Å is quite sizable (7.5 kcal/mol for STO-3G and 15 kcal/mol for cc-pVTZ, respectively). RAFQMC again performs significantly worse than RCCSD(T). UAFQMC is quantitatively comparable to UCCSD(T) for both bases sets, but the sign of its error is opposite to that of UCCSD(T In summary, for bond breaking examples and other typical strongly correlated systems, it appears that RAFQMC exhibits significant errors, much larger than what is seen with RCCSD(T). However, UAFQMC is quantitatively better than UCCSD(T), although it is still far from chemical accuracy. This can be seen from the non-parallelity errors presented in Table III. Nonparallelity error is a commonly used metric for assessing the performance of different methods in the computation of the relative energetics on a potential energy surface. Therefore among single-reference methods, UAFQMC is a promising candidate for the treatment of relatively simple strongly correlated systems.

Assessment of multi-reference methods
We repeat the same analysis for MR methods where simple CASSCF states are employed either as a reference state or a trial wave function, and quantify the errors associated with the treatment of dynamic correlation. In Fig. 9, it is clear that the worst performer in terms of absolute energies in both H 2 and N 2 is NEVPT2. CASPT2 and MRMP2 perform much better than NEVPT2, by approximately 3 kcal/mol for H 2 and by more than 10 kcal/mol for N 2 . Nonetheless, in terms of non-parallelity error (NPE), see Table IV, NEVPT2 is comparable to CASPT2 and MRMP2 for H 4 and is better than these two methods for N 2 by more than 6 kcal/mol. CASPT2 and MRMP2 methods perform quite similarly (identically for H 4 ), and their performance is far from the level of chemical accuracy. The performance of DSRG-MRPT2 falls between that of NEVPT2 and CASPT2/MRPT2. The poor performance of secondorder perturbation theories is particularly worrying, as these methods are the most commonly employed MR methods due to their simplicity and their relatively inexpensive nature, namely a scaling of O(N 5 ) outside the active space. We have also investigated the performance of MR-CISD using different ways to correct for its sizeinconsistency error (namely MRACPF, MRAQCC, MR-CISD+Q). These methods are computationally more expensive than MRPT methods, with a scaling of O(N 6 ) outside the active space. This scaling also applies to two DSRG methods, DSRG-MRPT3 and MR-LDSRG(2). MRACPF, MRAQCC, and MRCISD+Q are all nearly exact for H 4 . MRAQCC becomes more inaccurate for N 2 , exhibiting errors beyond chemical accuracy at small bond distances. Nonetheless, MRACPF and MRCISD+Q provide near-chemical acuracy for all bond distances for N 2 . For both molecules, MR-LDSRG(2) is not as accurate as the MRCI methods for both absolute and relative energies. Nonetheless, for H 4 , its NPE is only 0.73 kcal/mol. DSRG-MRPT3 performs very similarly to MR-LDSRG(2) for H 4 , but it differs significantly from MR-LDSRG(2) for N 2 For N 2 , DSRG-MRPT3 produces absolute energies that are quite far from chemical accuracy, but its NPE is quite small and comparable to that of MRCISD+Q and MRAQCC.
AFQMC with CASSCF trial wave functions works remarkably well for both H 2 and N 2 achieving chemical accuracy at all bond distances. This accuracy comes with a far lower scaling cost (see Table I) than both MRPT and MRCI-based methods, where we emphasize the steep scaling of these methods within the active space due to the requirement of higher-order density matrices. We also note that none of the other MR methods except NEVPT2 and the DSRG methods are size-consistent. Since CASSCF is size-consistent, AFQMC with CASSCF trial wave functions is also size-consistent in the limit of ∆t → 0 (see Section II D).
Due to its computational efficiency and its accuracy in reproducing both absolute energies and NPE values (see Table IV), AFQMC would appear to be the method of choice when strong correlation can be handled by the CASSCF trial wave function and the remaining dynamic correlation needs to be corrected at scale. With the improved algorithm for using SCI trial wave functions, 9 AFQMC will likely become a commonly used MR dynamic correlation method in the future.

IV. OUTLOOK AND FUTURE OPPORTUNITIES
There are many future opportunities for developments and improvements of different aspects of ph-AFQMC. Due to its relatively new vintage and the comparatively limited use of the method in quantum chemistry, the development of ph-AFQMC when compared to, for example, CC methods has not yet reached full maturity. We note several promising opportunities for future development in this section.

Trial wave functions:
As mentioned in Section II B, the form of trial wave functions is currently limited to SD trials or relatively compact MSD trials. Using more accurate trial wave functions without greatly increasing the overall com-putational cost is a forefront area for future ph-AFQMC development. Compared to CCSD(T), where going beyond the SD reference is formally challenging, the flexibility afforded by ph-AFQMC will lead to important developments along this line. Our experience suggests that the sensitivity to the trial wave function within ph-AFQMC is far greater than is the sensitivity of CC methods to the reference wave functiton. This both provides opportunities to greatly enhance accuracy but also challenges for rationally developing strategies for the use of trial functions.
2. Alternative constraints: Due to the lack of extensive benchmark studies, we do not know if the constraint provided by the cosine projection in Eq. (7) is optimal in terms of accuracy and statistical efficiency. As noted in Section II D, the currently used phaseless constraint introduces an O(∆t) time step error, limiting the size of time steps and necessitating time step extrapolations beyond the Trotter error. Furthermore it may well be that different constraints provide different levels of accuracy even when used with the same trial wave function. Much more work in this direction is warranted.
3. Observables other than ground state energies: Some observable properties extracted from ph-AFQMC have been shown to be less reliable than is desirable. 25,80 The accuracy of properties computed via back propagation provides one concrete example. Furthermore, excited state energies also fall into this category where additional complications arise due to the nature of projector Monte Carlo. At this stage ph-AFQMC is an accurate method for the ground state energy in chemical systems, however ancillary algorithms that enable the computation of other properties with the accuracy consistent with that of the ground state energy within ph-AFQMC is an important goal for future research. 4. Computational cost: Local energy evaluation algorithms for ph-AFQMC with SD trial wave functions are relatively mature and provide multiple ways to reach cubic scaling per sample. 5-8 However, algorithms to accelerate walker propagation are scarce. In addition, local energy evaluation algorithms for ph-AFQMC with more complex MSD trial wave functions have been relatively less explored. 8,9 Thus, reducing the computational cost of all components in ph-AFQMC other than the local energy evaluation forms another direction worthy of future effort.

V. CONCLUSIONS
The purpose of this work has been two-fold: First, we have presented a self-contained overview of the formalism behind the ph-AFQMC approach from a quantum chemistry perspective, and have delineated the considerations needed to understand the computational implementation and cost of the approach. Second, we have assessed the performance of ph-AFQMC for a well-known main group chemistry benchmark set (W4-11) and a non-covalent interaction benchmark set (A24), leading to a total of 1004 relative energies. In addition, we have studied the potential energy curves of commonly studied model problems, namely H 4 and N 2 , documenting the performance of ph-AFQMC with different trial wave functions against other standard wave function-based quantum chemistry methods. This constitutes the largest quantum chemical benchmark study to date using ph-AFQMC. While the knowledge gained from our work is far from complete, we will conclude our work with some cautious recommendations for the use of ph-AFQMC in a broad chemical context.
Single Slater determinant trial wave functions: ph-AFQMC employed with a with single determinant trial is an accurate method that can be a method of choice if CCSD(T) is too expensive and if the target system does not exhibit open-shell electrons with sizable antiferromagnetic coupling. While studies on chemical systems using DMC and GFMC have been limited, AFQMC was shown to be more accurate than these for the same single determinant trial for several specific systems. 15,153 Such trials will generally take the form of HF wave functions, or DFT 10,66 or regularized OOMP2 16,113 wave functions when HF exhibits unphysical properties such as overlocalization. The performance of ph-AFQMC for noncovalent interaction energies is still unclear due to the limited number of data points we have considered. We hope that the community will take up the task of a more detailed investigation in the near future.
Multi-Slater determinant trial wave functions: For systems exhibiting open-shell electrons with antiferromagnetic coupling either due to bond breaking or arising from localized d-or f-electrons, ph-AFQMC must be used with more elaborate multi-Slater trial wave functions such as CASSCF or SCI wave functions. Using spin-unrestricted single determinant trials can be an option if there are only two electrons contributing to the strong correlation behavior. 16,154 However, the phaseless error needs to be carefully quantified by more sophisticated trial wave functions in realistic strongly correlated if possible. We stress that ph-AFQMC was found to be far more accurate than other MRPT methods, such as CASPT2, NEVPT2, and MRMP2. Given its high efficiency compared to MRCI methods, we anticipate ph-AFQMC to become a clear method of choice for including dynamic correlation in active-space calculations. With improved algorithms 9 and implementations leveraging graphical processing units (GPUs), 28,76 we believe that ph-AFQMC with ∼ 10 5 to 10 6 determinants will become routine calculations. Nonetheless, enabling the efficient use of sophisticated trial wave functions is one of the most pressing topics for near-term development of ph-AFQMC. 37,51,118,119 We hope that the numerical data presented in this work, as well as the existence of currently available opensource packages 81-83 will encourage many more developers and users to explore uncharted territory and contribute to method development within the ph-AFQMC framework. While the data provided here constitutes the most extensive data set for ph-AFQMC to date, it is important to continue producing data and comparisons for a wide class of systems. Furthermore, compared to other deterministic quantum chemistry methods, ph-AFQMC is not yet at the stage of development where it can be considered a "black-box" approach. In addition, its final energy is also statistical in nature, introducing additional barrier for users to use the method. It is, therefore, important to have and further develop open-source frameworks 82,83 and to provide straightforward user interfaces. Such efforts will enable the expansion of our understanding of the relative benefits and weaknesses of ph-AFQMC, and we hope this goal is taken up by the electronic structure community with enthusiasm in the coming years.

VI. ACKNOWLEDGMENT
We thank Luke Bertels, Garnet Chan, Francesco Evangelista, Fionn Malone, Adam Rettig, Sandeep Sharma, and Shiwei Zhang for helpful discussions. DRR would like to thank Richard Friesner and Shiwei Zhang for collaboration on AFQMC-related topics over the last several years. We thank Richard Friesner for the suggestion to carefully consider basis set extrapolations in this work. We thank Google for their gift that was used to support JL and HQP in developing a code used in this work. We also thank Francesco Evangelista for providing the bond dissociation data of H 4 and N 2 for the DSRG methods. We acknowledge computing resources from Columbia University's Shared Research Computing Facility project, which is supported by NIH Research Facility Improvement Grant 1G20RR030893-01, and associated funds from the New York State Empire State Development, Division of Science Technology and Innovation (NYSTAR) Contract C090171, both awarded April 15, 2010.
where N is the total number of imaginary time steps and the small imaginary time step is defined as ∆τ = τ /N . The efficiency of AFQMC relies on the representation of the short-time propagator. One first recastsĤ in Eq. (2) intoĤ whereĤ 1 is the one-body part of the Hamiltonian, and the two-body part is written as the sum of N γ squared operators, for example, by using the modified Cholesky decomposition 155,156 or density fitting, 157,158 v with In practice, we modify this Hamiltonian further by performing a mean-field subtraction 18 or using a shifted contour 36 . This amounts to writinĝ To realize Eq. (1) in a computationally efficient manner, we write the global wave function as a weighted summation over a statistical sample of wave functions, {|ψ i }, commonly referred to as "walkers": where w i (τ ) is the weight of the i-th walker at imaginary time τ . We also expressB(x) in Eq. (3) aŝ An instance of the short-time propagator in Eq. (A8) is then applied to a set of random walkers {|ψ i } to obtain the ground state |Ψ 0 from the initial state. Namely, each walker is assigned to a vector of Gaussian random variables x i (τ ) at imaginary time τ and one applies Eq. (A8) to advance the wave function to the next time step This procedure effectively computes the highdimensional integral in Eq. (3) by means of Monte Carlo sampling. Single Slater determinants are the most common form of walker wave functions, and the Thouless theorem ensures that the walkers will stay in the single Slater determinant manifold upon the action of the propagator. 159,160 Global estimates for operators that commute withĤ are evaluated by the following mixed estimator where |Ψ T is a wave function that allows for a straightforward evaluation of the above expression. IfÔ =Ĥ, this provides an unbiased way to estimate the total energy. The algorithm discussed so far is defined as the freeprojection AFQMC (fp-AFQMC) which is exact in principle, however, suffers from the phase (sign) problem. 66 Due to this problem, the mixed estimate in Eq. (A10) will generically have a variance that scales exponentially with system size. 38 One way to control the severe statistical fluctuation arising from the phase problem is to use the phaseless approximation, 66 resulting in the phaseless AFQMC method (ph-AFQMC). ph-AFQMC utilizes importance sampling based on a trial wave function |Ψ T throughout the random walk process and thereby works with the statistical representation of global wave functions as in Eq. (4). Given this importance sampling, the walker weight update rule follows where the importance function I which is proportional to the overlap ratio in Eq. (6) is given by × e xi(τ )xi(τ )−xi(τ )xi(τ )/2 . (A12) x i (τ ) is usually referred to as the optimal force bias 97 which provides a shift to the underlying normal distribution. The optimal force bias is computed as where the one-body Green's function G qp,i (τ ) is given by Beyond Eq. (A11), ph-AFQMC imposes a constraint to ensure the positivity of the weights throughout the imaginary time propagation. Such a constraint is achieved by a modified importance function, I ph , where the phase θ i (τ ) is given in Eq. (8). This modified importance function is used in ph-AFQMC to update the weights. We verify the size-inconsistency error of ph-AFQMC with a finite time step via numerical results on N 2 at a bond distance of 1.3Å and with an RHF trial wave function in the cc-pVDZ basis. 161 We add N 2 molecules at infinite separation. One should expect the correlation energy per molecule to be constant if the method was size-consistent. We used 640 walkers and checked that this result is not due to the population control bias. As shown in Fig. A1, for a relatively large time step (∆t = 0.05 a.u.), the correlation energy per molecule decreases as we add more N 2 molecules to the system. Going from one N 2 molecule to 10 N 2 molecules, we observe the change in the correlation energy per molecule by 10 mE h for ∆ = 0.05 a.u. This confirms that ph-AFQMC is in general not size-consistent as illustrated by Eq. (16). We do not see such a large error in the case of a smaller time step (∆t = 0.005 a.u.). For larger systems, we anticipate that ∆t = 0.005 a.u. would not be small enough to completely eliminate this error, and thus one must be cautious with respect to time step and the possibility of size-inconsistency. We also emphasize that when quantifying the size-inconsistency error in this numerical experiment, the typical time step error caused by the Trotter error in Eq. (3) is the same for each fragment. Therefore, our result is not due to the time step error during the Trotterization, and is solely caused by the lack of size-consistency in ph-AFQMC for finite ∆t.

A3. Additional computational details concerning benchmark data
We used 640 walkers in all systems considered here and the corresponding population control bias was found to be negligible with the pair-branch algorithm. 162 A time step of 0.005 au was used and the time step error (for both Trotterization and size-consistency) was found to be also negligible at the energy scale we focus on in this work. The frozen core approximation was used throughout unless specified otherwise.
For the W4-11 set, 85 we used the aug-cc-pV5Z basis set 161 to converge the Hartree-Fock (HF) energies to the basis set limit and used aug-cc-pVTZ 161 and aug-cc-pVQZ 161 to extrapolate the correlation energies to the basis set limit following Helgaker's two-point extrapolation. 163 We also use Karton and Martin's inner-shell correlation energy contribution to correct for the missing core correlation energy due to the frozen core approximation. 85 Furthermore, to estimate the residual basis set error due to the differences in our scheme and Karton and Martin's scheme, we also compared our CCSD(T) relative energies and their reported CCSD(T) relative energies. The root-mean-square-deviation of the difference was 0.66 kcal/mol which is small enough for the purpose of this paper. We note that Karton and Martin used the composite W4 scheme to estimate reference energies which involved up to CCSDTQ5 or CCSDTQ6. For the A24 set, 86 we used aug-cc-pVTZ 161 with counterpoise corrections, although the reference data (i.e., CCSD with triples and perturbative quadruples) were obtained in the basis set limit. For H 4 and N 2 dissociation data, we used cc-pVQZ and cc-pVTZ, respectively, bases 161 without the frozen core approximation.
Molecular integrals for ph-AFQMC were generated by PySCF, 164 all deterministic single-reference quantum chemistry calculations were performed with Q-Chem, 165 and AFQMC calculations were performed with the ipie 82,83 and QMCPACK. 81 Cholesky factorization and integral transformations necessary for AFQMC calculations were performed by a script in ipie 82,83 with a Cholesky threshold of 10 −5 . Semistochastic heat-bath configuration interaction (SHCI) calculations were performed with Dice. 166 Multi-reference dynamic correlation methods other than driven similarity renormalization group (DSRG) methods were run on Orca. 167 DSRG methods were run through Psi4 with the flow parameter σ = 0.5 a.u.. 168 DSRG-MRPT2 and DSRG-MRPT3 energies were obtained via the partial relaxation scheme. 147 MR-LDSRG(2) energies were obtained by relaxing the reference state twice. 148 We used SHCI to generate near-exact benchmark energies for bond dissociation examples. SHCI energies are converged more accurately than 0.1 mE h , except for N 2 in the cc-pVTZ basis. For N 2 in the cc-pVTZ basis, we estimate the remaining bias in our SHCI energies to be less than 1 kcal/mol. The second-order perturbation theory (PT2) contribution was found to be at most 10 mE h for this system and we did not perform any extrapolation to reduce the potential bias further. The SHCI energies used in this work are within 1 mE h of the SHCI energies from our recent investigation of N 2 in the cc-pVTZ basis at a limited set of bond distances 119 where the PT2 contribution in SHCI was converged to better than 2 mE h . We believe that the residual bias in SHCI energies is small enough that our conclusions are not affected by this bias. ωB97M-V results are taken from ref. 2 All relevant raw energies are available in the Zenodo repository. 130

A4. Potential energy curves of H4 and N2 in STO-3G
We present the potential energy curve of H 4 and N 2 using an STO-3G basis set with single-reference methods in Fig. A2 and multi-reference methods in Fig. A3. This is to give the reader a visual sense for the potential energy curves for each of these problems.    A3. The potential energy curve of (a) H4 in cc-pVQZ and (b) N2 in cc-pVTZ as a function of the bond distance using multi-reference methods.. Note that there are more data points than markers shown.