Ab initio computations of molecular systems by the auxiliaryfield quantum Monte Carlo method
Abstract
The auxiliaryfield quantum Monte Carlo (AFQMC) method provides a computational framework for solving the timeindependent Schrödinger equation in atoms, molecules, solids, and a variety of model systems. AFQMC has recently witnessed remarkable growth, especially as a tool for electronic structure computations in real materials. The method has demonstrated excellent accuracy across a variety of correlated electron systems. Taking the form of stochastic evolution in a manifold of nonorthogonal Slater determinants, the method resembles an ensemble of densityfunctional theory (DFT) calculations in the presence of fluctuating external potentials. Its computational cost scales as a lowpower of system size, similar to the corresponding independentelectron calculations. Highly efficient and intrinsically parallel, AFQMC is able to take full advantage of contemporary highperformance computing platforms and numerical libraries. In this review, we provide a selfcontained introduction to the exact and constrained variants of AFQMC, with emphasis on its applications to the electronic structure in molecular systems. Representative results are presented, and theoretical foundations and implementation details of the method are discussed.
I Introduction
A central challenge in the fields of condensed matter physics, quantum chemistry, and materials science is to determine the quantummechanical behavior of many interacting electrons and nuclei. Often relativistic effects and the coupling between the dynamics of electrons and nuclei can be neglected, or treated separately. Within this approximation, the manyelectron wave function can be found by solving the timeindependent Schrödinger equation for the BornOppenheimer Hamiltonian Born and Oppenheimer (1927); Szabo and Ostlund (1989)
(1) 
where is the position of electron , and the position of nucleus , with charge . The numbers of electrons and nuclei are and , respectively. The nuclear positions are held fixed. We use atomic units throughout. The main obstacle to the investigation of chemical systems is that, in general, the computational cost of finding the exact ground state of Eq. (1) grows combinatorially with the size of the studied system Troyer and Wiese (2005); Schuch and Verstraete (2009). This limitation has so far precluded exact studies for all but small molecular systems and motivated the development of approximate methods.
By far the most widely used approximate methods are independentelectron approaches, based on the celebrated densityfunctional theory (DFT) Martin (2004); Kohn (1999). These are standard tools for electronic structure calculations in diverse areas across multiple disciplines, with sophisticated computer software packages available. The success of DFTbased approaches has been exceptional. Their difficulties are also welldocumented, especially in socalled strongly correlated systems, for example in many transitionmetal oxides.
For molecular systems, a hierarchy of quantum chemistry (QC) methods have been developed, which allow systematic improvement in accuracy, at increasing computational cost. For example, one of the most accurate QC methods is coupledcluster CCSD(T) Paldus and Li (1999); Bartlett and Musial (2007); Shavitt and Bartlett (2009), scaling as the seventh power of the system size. There are several very promising new approaches, including densitymatrix renormalization group (DMRG) White (1992a); White and Martin (1999); Chan and HeadGordon (2002); OlivaresAmaya et al. (2015); Chan et al. (2016) and fullconfiguration interaction quantum Monte Carlo (FCIQMC) Booth et al. (2009, 2013a), which have dramatically expanded the size of the molecular systems that can be treated with very high accuracy. However the computational cost of these methods, in the most general case, still scales exponentially with system size. A recent benchmark study has been performed on the linear hydrogen chain Motta et al. (2017), including many QC methods, diffusion Monte Carlo Reynolds et al. (1982a); Foulkes et al. (2001a), and auxiliaryfield quantum Monte Carlo (AFQMC).
As a general electronic structure approach, the AFQMC method Zhang and Krakauer (2003) has several unique and attractive features. We hope that this will become evident following the description of the method below. As a brief overview, AFQMC combines stochastic sampling with the machinery of DFT. It is a nonperturbative approach that balances accuracy and computational scaling. Utilizing a fieldtheoretic framework and casting the projection of the manybody ground state or the finitetemperature partition function as a random walk in nonorthogonal Slater determinant space, AFQMC takes the form of a linear combination of DFT calculations in fluctuating auxiliary fields. It can incorporate many techniques used by independentparticle methods such as DFT or HartreeFock (HF). AFQMC is naturally parallelizable and very well suited for largescale computing platforms. Though more demanding than DFT or HF, AFQMC has computational cost growing as the third (in some cases fourth) power of the system size, enabling applications to large systems.
The review is structured as follows. In the next section we discuss the AFQMC method in detail: the key technical backgrounds are described in several subsections; this is followed by discussions of the sign and phase problem, and of the phaseless AFQMC approach which removes the phase problem in electronic systems; a stepbystep algorithmic outline is then provided. Then in the next section, we describe several recent methodological advances with AFQMC, including frozencore, downfolding and embedding approaches, backpropagation to compute observables and correlation functions in electronic systems, and correlated sampling to improve efficiency in computing energy differences. Finally we briefly discuss the outlook of AFQMC as an electronic structure method, the advantages it brings and some directions for future development. Additional mathematical background and algorithmic details are included in the Appendix.
Ii The AFQMC method
Using the formalism of second quantization, the Hamiltonian (1) can be written as
(2) 
where is the internuclear repulsion, are fermionic creation and annihilation operators associated to an orthonormal basis of oneelectron molecular orbitals (MOs). These basis functions are typically real (although complex basis functions can be handled straightforwardly in AFQMC), and
(3) 
are the matrix elements of the one and twobody parts of (2) in that basis. is often denoted as in chemistry, and in physics.
Most applications of AFQMC in molecular systems so far have relied on basis sets of MOs obtained by orthonormalizing a set of atomcentered atomic orbitals (AOs) expressed by standard Gaussian basis functions Szabo and Ostlund (1989). The overlaps and matrix elements , in the AO basis are computed and output by standard quantum chemistry softwares Schmidt et al. (1993); Valiev et al. (2010); Sun et al. (). The matrix elements of Eq. (3) result from a straightforward change of basis.
The use of a finite basis set unavoidably introduces an approximation, which can be removed by extrapolating results to the complete basis set (CBS) limit. This is common to all quantum chemistry methods. For Gaussian bases, it is necessary to perform calculations with increasingly large basis sets designed to allow systematic convergence to the CBS limit. AFQMC scales as to , as further discussed below. Compared to highlevel quantum chemistry methods such as CCSD(T), this is advantageous for reaching the CBS. We typically use wellestablished bases (e.g. ccpVxZ for light atoms Jr. (1989); Woon and Jr. (1993), or ccpVnZPP Peterson (2003); Peterson et al. (2003)) and empirical extrapolation formulas Feller (1992); Helgaker et al. (1997) for CBS extrapolations.
AFQMC allows for the choice of any oneelectron basis suitable for the problem. In addition to Gaussians, AFQMC has also been implemented using planewave and pseudopotentials Suewattana et al. (2007); Zhang et al. (2005); Kwee et al. (2008); Purwanto et al. (2009a); Ma et al. (2015, 2017), more suited for calculations in solids. As we discuss below, one could also use KohnSham orbitals derived from a planewave DFT calculation as basis sets, which provides a kind of hybrid between the planewave and the MO approach.
ii.1 Properties of (nonorthogonal) Slater determinants
We first review some properties of Slater determinants (SDs) important to the formulation and implementation of AFQMC. A single SD is written as
(4) 
for and particles with spin up and down respectively, occupying the orbitals and . SDs are exact eigenfunctions of onebody Hamiltonians, including those of noninteracting systems and independentelectron Hamiltonians from DFT or meanfield HF.
There exists a key distinction between the SDs in AFQMC and in standard quantum chemistry methods. In quantum chemistry, SDs involved are orthogonal, as they are built from a common set of reference orbitals (e.g. occupied and virtual HF orbitals) through different excitations of the electrons. In AFQMC, SDs are nonorthogonal to each other, as the occupied orbitals of each SD are continuously changed (via rotations of the occupied orbitals) during the random walk. Some properties of nonorthogonal SDs are summarized below. (Additional details can be found in the Appendix.)
Theorem 1.
The following properties hold for any two nonorthogonal SDs and , parametrized by the matrices and , respectively, and for any spinindependent oneparticle operator

Overlap: .

Oneparticle reduced density matrix:
(5) 
Properties 2 and 3 remain invariant under orthonormalization of an SD, , where the orbitals in are obtained by orthonormalizing those in (e.g., by a modified GramSchmidt procedure).
ii.2 Groundstate projection
The phaseless AFQMC method is a projective quantum Monte Carlo (QMC) method. Projective methods are based on the observation that the ground state of a Hamiltonian is the asymptotic solution of the imaginarytime Schrödinger equation
(8) 
where is any initial wavefunction nonorthogonal to , and is the groundstate energy of . The formal solution of (8) is obtained applying the imaginarytime propagator to ,
(9) 
The groundstate energy is unknown but it can be estimated adaptively as the calculation progresses, for example with the growth estimator discussed in the Appendix.
In general it is often as difficult to realize the imaginarytime propagation in Eq. (9) as it is to solve the timeindependent Schrödinger equation. However, for an independentelectron method such as HF, where is replaced by
(10) 
Eq. (7) guarantees that the imaginarytime projection turns into the motion of a single SD along the manifold of Slater determinants, as illustrated in Fig. 1. (Of course a special case would be the noninteracting system, when the 2body part of the Hamiltonian vanishes.) Similarly, a DFT calculation can be cast in this way, with the projection in imaginarytime as an alternative for the nonlinear minimization of the energy .
AFQMC maps the imaginarytime projection of a manybody Hamiltonian onto a stochastic process in the manifold of Slater determinants, as sketched in Fig. 1. The procedure to realize this mapping is outlined below. There are different ways to view the mapping. It can be seen as a generalization of the aforementioned single SD projection in DFT or HF, where manybody correlation effects is recovered by random sampling around a typical “meanfield”. Alternatively, the mapping can be seen as a different and complementary approach to fixednode diffusion Monte Carlo Reynolds et al. (1982b); Foulkes et al. (2001b), which simulates the imaginarytime projection as a random walk in electron coordinate space, or Green’s function Monte Carlo Trivedi and Ceperley (1989); Calandra Buonaura and Sorella (1998) and FCIQMC Booth et al. (2009, 2013a), which carry out a similar projection [using ] in fermion configuration space (i.e., orthogonal SDs). As we show below, the use of overcomplete, nonorthogonal SDs in AFQMC significantly reduces the severity of the QMC fermion sign problem Feynman et al. (2010); Loh et al. (1990); Schmidt and Kalos (1984) in many situations.
ii.2.1 The HubbardStratonovich transformation
The first step to provide a stochastic representation of the imaginarytime projection is often a decomposition in shorttime propagators,
(11) 
Within AFQMC, one first expresses the Hamiltonian as
(12) 
with , onebody operators, explicitly given in the next Section. The shorttime propagator is then approximated by the use of a Trotter breakup Trotter (1959); Suzuki (1976)
(13) 
and the HubbardStratonovich transformation Hubbard (1959); Stratonovich (1958) (see Appendix),
(14) 
Combining Eqs. (13) and (14) yields the following representation of the shorttime propagator,
(15) 
with and
(16) 
The shorttime propagator is thus written as an integral over normally distributed parameters , called auxiliary fields, of exponentials of onebody operators, . Eq. (15) maps the original interacting system onto an ensemble of noninteracting systems coupled with fluctuating external potentials. The integration over the auxiliary fields recovers the effect of the twoparticle interaction. While for a completely general twobody Hamiltonian the number of auxiliary fields is AlSaidi et al. (2006a); Zhang (2013); Motta et al. (2014a), for real materials can be brought down to , an example of which is given next.
ii.2.2 Modified Cholesky decomposition
To perform the HS transformation described above, we can rewrite the Hamiltonian as
(17) 
Representing (exactly or approximately) as a sum of terms in which the indices , are uncoupled immediately leads to the form (12), with defined in (17) and
(18) 
There are many possible ways of manipulating , and taking advantage of this flexibility is essential to improving the efficiency and accuracy of AFQMC. When realvalued orbitals are used as oneelectron basis set, the interaction tensor can be reshaped into a realvalued, symmetric and positivedefinite matrix , which can then be approximated using a modified Cholesky decomposition (mCD) algorithm Beebe and Linderberg (1977); Koch et al. (2003); Aquilante et al. (2010). The goal of the mCD is to recursively reach a value when the maximum error in representing is less than a predetermined threshold value, . Given Cholesky vectors , can be expressed as
(19) 
where is the residual error matrix at the th iteration. Let denote the index of the largest diagonal element of :
(20) 
A damped prescreening operation is performed, i.e., elements satisfying
(21) 
are put equal to zero. As noted in Refs. Koch et al. (2003); Aquilante et al. (2010), the damping serves as a safeguard against rounding errors that may render the decomposition numerically unstable. The nextiteration Cholesky vector is then obtained as
(22) 
and the residual error matrix is updated accordingly. Note that the positivity of and the recursion relation (22) imply that is always nonnegative. When , the CauchySchwarz inequality implies for all , and the iteration terminates.
Typical values of range between and . Performing a mCD decomposition of requires operations and memory Purwanto et al. (2011), and leads to Cholesky vectors, as exemplified in Fig. 2. Direct diagonalization of , used in early implementations of AFQMC, leads instead to vectors, and requires operations and storage, which would be prohibitive for large molecules.
ii.2.3 Meanfield subtraction
Different forms of the HS transformation exist which can affect the performance of AFQMC. It is computationally advantageous to rewrite the Hamiltonian (12) prior to the HS decomposition, by subtracting an estimate of the groundstate average of each operator ,
(23) 
where the modified onebody operators are , and the constant shift can be obtained, for example, using the trial wave function,
(24) 
This is called meanfield background in the AFQMC literature. HS transformation applied to the modified onebody operators can lead to smaller fluctuations in freeprojection calculations and smaller systematic errors when a constraint is applied AlSaidi et al. (2006a); Purwanto et al. (2009b).
ii.3 The freeprojection AFQMC
With the formalism established so far, it is straightforward to realize the projection process by random walks in the space of and evaluate the groundstate energy:
(25) 
where is a trial wavefunction, typically the SD from HF or DFT, or a linear combination of SDs from a completeactive space selfconsistent field (CASSCF) calculation, for example. For convenience, is often chosen to be the same as , although this is not required. (We can take advantage of the flexibility to choose them differently in order to help impose certain symmetry properties in the calculations, see below and Refs. Purwanto et al. (2008); Shi and Zhang (2013)). We have written and to reduce clutter.
One way to evaluate the integrals in Eq. (25) is by sampling a fixedlength (of segments) path in the auxiliaryfield space, for example, by the Metropolis Monte Carlo (MC) algorithm. This is the form in which the auxiliaryfield approach was first proposed, both at finitetemperature Blankenbecler et al. (1981) and at Sugiyama and Koonin (1986), in model systems. Here we will describe a different approach by casting (25) into branching, openended random walks in the space of SDs. For freeprojection calculations per se, there is no particular advantage in using this approach over the more standard Metropolis procedure; in fact, it is more convenient to compute observables or correlations with the latter. However, in order to impose a constraint to control the sign problem (and phase problem) as we will discuss next, there is a fundamental difference in how the two algorithms scale, and it was necessary to recast the process as openended random walks Zhang et al. (1995).
In the random walk realization of Eq. (25), the wavefunction is represented, in a MC sense, by an ensemble of SDs, with weights and phases . The structure is called a random walker. The calculation begins by initializing a population of walkers, at . The population size is either held fixed at a preset value or allowed to fluctuate around it throughout the calculation, applying population control as needed.
As the random walk proceeds, for all walkers and times , we sample an auxiliary field and update the walker state, weight and phase as
(26) 
By Thouless’ theorem, is a SD if is a SD, so that all walkers maintain simple structure and efficient parametrization during the random walk. The groundstate energy can be computed as weighted average of the local energy functional,
(27) 
The freeprojection AFQMC estimator of the groundstate energy is exact. The walkers will need to be periodically reorthonormalized to maintain numerical precision. The only potential systematic errors come from the Trotter error of finite timestep, and population control bias. These can be easily removed by standard extrapolation procedures, as illustrated in Fig. 3. The Trotter error in freeprojection is typically quadratic in . Population control bias is typically , and should be small in a freeprojection calculation, since the population size is necessarily large; however, the calculation is prone to large fluctuations in the weights because of intrinsic instability from the sign/phase problem.
A detailed flowchart of the algorithm will be presented after constrained AFQMC calculations are introduced. Additional implementation details, including population control and reorthonormalization of the walkers, which are shared by freeprojection and constrained calculations, will be discussed under Implementation issues.
Iii Constrained AFQMC calculations
iii.1 The phase problem
The freeprojection AFQMC suffers from an asymptotic instability in , which is the manifestation of the fermion sign problem Feynman et al. (2010); Loh et al. (1990) in many model systems such as the doped repulsive Hubbard model. For more general interactions, such as the longrange Coulomb repulsion, a more severe instability appears Zhang and Krakauer (2003), which has been referred to as the phase problem in AFQMC literature. More specifically, in freeprojection AFQMC, is estimated by a ratio whose denominator is the expectation , with , with respect to the random variable over the manifold of SDs. (The numerator is given by .) Unless prevented by special symmetry, the expectation of the phase approaches zero exponentially with . This causes the variance in the MC estimator of the energy to grow exponentially, . The rate at which the variance increases depends on the nature of the interaction. For each system, the rate grows with and , making the instability of freeprojection AFQMC more severe for large basis sets and systems. As the Central Limit Theorem dictates, computing the average to a desired accuracy requires the number of samples to scale as . The exponential increase of noise with imaginary time translates into a formal exponential scaling of freeprojection AFQMC. The exponential growth in MC noise from the phase problem is illustrated in Fig. 3.
To gain a sense of how the phase problem arises, let us recall that we are seeking a representation of the solution in Eq. (9) in the form
(28) 
in the overcomplete space of nonorthogonal SDs, where the coefficients are represented stochastically by the probability distribution of MC samples of Fahy and Hamann (1990, 1991). The manifold of SDs contains both a generic determinant and its opposite, , as distinct points. Indeed it contains for any . This symmetry property has its root in the linearity of the Schrödinger equation: if is a solution to Eq. (8), so is ; similarly, if is a solution in Eq. (28) defined on a domain in the SD manifold comprising , so is , defined on a domain comprising .
In the framework of a numerical simulation, the symmetry between equivalent domains in implies that, for any statistical ensemble of walkers giving an MC representation of the groundstate wavefunction, there exist other ensembles providing equivalent representations. If is real, there are two equivalent ensembles, related by a sign change. If is complex, there are infinitely many ensembles given by different gauge choices.
Apart from special cases, such as the halffilled repulsive Hubbard model Hirsch (1985), where particlehole symmetry yields a that confines the random walk in one domain with a single choice of sign (either or ), the random walks will in general reach other domains. The MC dynamics given by does not have any means to detect the crossing of “boundaries” between different domains: the walkers evolve continuously in the random walk Zhang (2013); because of the high dimensionality of , we cannot tabulate a history of where the random walk has visited to estimate a sign or phase (in the spirit of cancellation algorithms Zhang and Kalos (1991); Anderson et al. (1991); Booth et al. (2009)), without incurring exponential scaling Zhang (2013). The exploration of being ergodic, in the limit the MC representation of the groundstate consists of an equal mixture of walkers from degenerate domains, regardless of whether the initial distribution was concentrated in one or not. This symmetry results in an exponentially fast decay in the signaltonoise ratio, i.e. in the average sign or phase, given by the expectation value .
In this framework, control of the sign problem can be achieved by modifying the dynamics so as to break the symmetry in the imaginarytime evolution or, equivalently, to prevent Slater determinants from crossing into other domains.
For real , the necessary symmetry breaking is achieved by deleting walkers whose overlap with the trial wavefunction changes sign, i.e., when crosses the hypersurface . The resulting approximation defines the constrainedpath AFQMC (cpAFQMC) method, which has delivered excellent results for lattice models of correlated electrons LeBlanc et al. (2015).
The phase problem affecting AFQMC calculations for real materials emerges as the HubbardStratonovich transformations result in complexvalued operators. (It is possible to maintain real values but these were found to have much larger fluctuations Zhang and Krakauer (2003)). The constraints used to remove that instability required an additional step beyond the constrainedpath approximation. The approach to real materials has been referred to in the literature as the phaseless or phasefree approximation. Once the latter is formulated, the former can be obtained as a special case. Here in discussing both these situations, where a constraint is applied on the sampling of the auxiliaryfield paths in the random walk, we will often refer to them as cpAFQMC for simplicity. Below we will introduce the key ingredients that enable the formulation of the cpAFQMC algorithm.
iii.2 Importance sampling transformation
Importance sampling is a variance reduction technique widely used in classical and quantum Monte Carlo Hammersley and Morton (1954); Rosenbluth and Rosenbluth (1955). Within AFQMC, importance sampling was introduced in the case of real propagators as a variancereduction measure Zhang et al. (1997); Zhang (1999). For complex , the importance sampling transformation is applied not just in the usual sense of variance reduction, but to define a gauge condition Zhang and Krakauer (2003).
We will outline the importancesampling transformation for the most general case below. The case of realvalued , such as in the Hubbbard model, will follow as a special example. Recalling that the HubbardStratonovich transformation is defined up to a complexvalued shift Zhang and Krakauer (2003); Purwanto and Zhang (2004), we can rewrite
(29) 
where now . Formally, we can apply a similarity transformation in each step along the random walk path using the overlap ratio , to obtain an equivalent form:
(30) 
where the importance function is defined as
(31) 
Note that the overlap ratio applied in our similarity transformation is in general complex, which is different from standard importancesampling transformations. Introducing a shorthand for the local “mixed” estimator: , and restricting our (arbitrary) choices of to be analytic functions of , we arrive at the following
(32) 
by expanding in a Taylor series. It is evident that the choice
(33) 
cancels fluctuations in the importance function to first order in , leading to a more stable random walk. The dynamic shift, which we will refer to as a force bias, modifies the MC sampling of the auxiliaryfields by shifting the center of the Gaussian (or Bernoullian in the case of discrete Ising fields Shi and Zhang (2013)) according to the overlap function . For complex , in general the force bias has an imaginary component. This is connected to the complex character of the overlap function in the similarity transformation discussed above, and turns out to play a key role for obtaining an optimal constraint to control the phase problem.
Inserting the optimal force bias into and expanding to first order in yields the expression
(34) 
Since the last term of (34) has zero average over auxiliary fields, leading to the approximate expression
(35) 
for the importance function. This is often referred to as the local energy formalism while Eq. (31) is referred to as the hybrid formalism.
By the introduction of the force bias, the importance sampling transformation guides the sampling of walkers away from regions of where , as diverges. In the case when is real, the procedure, taken at the limit of , imposes an infinite barrier at , and prevents the random walk from entering the other domain, thereby removing the sign problem. In the case of a phase problem, however, the transformation and the force bias alone do not lead to complete control of the phase problem. Achieving this goal requires combining importance sampling with the phaseless approximation described next.
iii.3 The phaseless approximation
From the expressions of the local energy and importance functions, (27) and (35), it is easy to see the role of in the phase problem: both quantities involve fractions with denominator , which diverges at the origin of the complex plane causing large fluctuations in AFQMC estimators. As mentioned, in the presence of complexvalued operators, the random walk dynamics is invariant under the transformations , . As the projection time grows, the overlaps undergo a rotationallyinvariant random walk in the complex plane, resulting in a finite concentration of walkers at the origin Zhang (2013).
The phase problem affecting AFQMC for real materials is unique, in that not only does the weight of the walker acquire a phase, but the orbitals of the SD become complex. Straightforward generalizations of the constrained path approximation, e.g., imposing the condition without the use of a complex force bias, lead typically to poor results Zhang and Krakauer (2003). The complex force bias of Eq. (33) keeps the overall phase stable to leading order in , thereby allowing a particular overall gauge to be chosen to break the rotational symmetry.
In the framework of electronic structure AFQMC calculations, control of the phase problem is achieved by the real local energy and the phaseless approximations. The real local energy approximation replaces the term in the importance function with its real part, thereby leading to real and positive weights. The impact of this approximation was found to be typically mild, and the imaginary part of the local energy can often be carried for extended projection time, with little effect on the computed energy. (See below on the effect for backpropagation and observables Motta and Zhang (2017a), however.) The phaseless approximation aims to break the rotationally invariance of the random walk undergone by the overlaps by projecting them onto an adaptively evolving straight line. The projection is accomplished by changing the importance function as follows:
(36) 
where is the phase difference between the overlaps, with the new walker state given by . The factor prevents stochastic trajectories from undergoing abrupt phase changes, and ensures that the probability distribution of vanishes at the origin.
With the optimal complex force bias, the ’’ form can be replaced by other choices which yield similar results Zhang et al. (2005). Notice that, when the operators are realvalued, the projection reduces to the constrainedpath AFQMC, because or and , respectively.
The resulting method, the phaseless AFQMC (phAFQMC), is the stateoftheart technique for performing AFQMC calculations in real materials. A phAFQMC calculations differs from a fpAFQMC in the update of walkers and weights, which is carried out as follows:
(37) 
with the force bias and importance functions given by (33), (36), and .
iii.4 Implementation issues
iii.4.1 The trial wavefunction
The choice of is important in AFQMC. First, the closer is to the ground state, the smaller the fluctuations in the local energy (if , ). In constrained calculations, also determines the efficiency of importance sampling and the accuracy of the phaseless constraint. Often using a single Slater determinant, such as the HF or DFT wavefunction, delivers results of excellent accuracy. Use of multideterminant trial wavefunctions,
(38) 
can improve accuracy, reduce local energy fluctuations and the imaginary time needed to reach equilibration of walker ensembles.
Theorem 1 can be straightforwardly generalized to multideterminant trialwavefunctions, by treating each configurations separately in the linear combination. Overlaps and Green’s functions can be computed as and
(39) 
For many practical purposes, it is useful to write
(40) 
where and are the matrices parametrizing . For multideterminants from many standard QC approaches, such as CASSCF, the sum over can be greatly accelerated using the ShermanMorrisonWoodbury formula Sherman and Morrison (1949); Press (2007).
iii.4.2 Propagation of walkers
Updating a walker requires computing the force bias . The cost of the operation can be brought down to by taking advantage of (40),
(41) 
where the tensor can be precomputed at the beginning of the simulation, and the trace in (41) costs operations. (We apply a cap to the force bias values if necessary Purwanto et al. (2009a) as discussed in the Appendix.)
Once the force bias is determined, the matrix is computed, at the cost of operations, and its exponential is applied to the matrices parametrizing . Computing as truncated Taylor series requires operations.
iii.4.3 Local energy evaluation
The local energy is needed for computing the groundstate energy, and it also controls the weights in constrained calculations. The most demanding part of its calculation comes from the interaction term, which via Wick’s theorem can be written as
(42) 
The scaling of (42) can be brought down to . Indeed, a simple calculation shows that the expression inside the square bracket in Eq. (42), the contribution to the local energy from the Cholesky vector and the configuration in can be written as
(43) 
The cost of the calculation results from computing onthefly the tensors , requiring operations, and the traces in (43), requiring operations. Additional reductions can be achieved for CASSCF type of . (A cap can be placed on the local energy to regulate spurious fluctuations Purwanto et al. (2009a) as discussed in the Appendix.)
iii.4.4 Population control
As the random walk proceeds, some walkers accumulate very large weights while some obtain very small weights. These different weights cause a loss of sampling efficiency because the algorithm spends a disproportionate amount of time keeping track of walkers that contribute little to the energy estimate. To eliminate such inefficiency, a branching scheme is introduced to redistribute weights without changing their statistical distribution. In such a scheme, walkers with large weights are replicated and walkers with small weights are eliminated with appropriate probability. However, because branching might cause the total population to fluctuate in an unbounded way (e.g. to grow to infinity or to perish altogether), we perform population control periodically to eliminate this instability. The population control can incur a bias when the total weight of the walkers is modified, as illustrated in Fig. 3. There exist approaches to reduce this bias by carrying a history of overall weight correction factors Calandra et al. (1998); Zhang et al. (1997); Umrigar et al. (1993). In typical cpAFQMC calculations, population sizes of order hundreds or thousands are used, and the population control bias is often negligible.
iii.4.5 Numerical stabilization
Repeated multiplications of the matrices onto the walker SD leads to accumulation of roundoff errors and a loss of numerical precision. This introduces spurious components in the walkers’ orbitals which eventually overwhelm the exact components. The instability may be interpreted as a tendency of walkers to collapse to a Bose ground state, in which all particles occupy the lowestenergy orbital. Within AFQMC, the instability is eliminated by periodically orthogonalizing the orbitals with, e.g., a modified GramSchmidt (mGS) procedure, by sweeping through the orbitals and of with the following steps:
(44) 
for , and similarly for the spindown orbitals. The mGS does not affect force biases or estimators of physical properties, thanks to Theorem 1.5. We comment that, in coordinate space methods such as diffusion MC or Green’s function Monte Carlo (GFMC), a similar tendency is present which is in fact the primary source of the sign problem in that framework. The procedure to stabilize the SDis in AFQMC is analogous to the cancellation of walkers in GFMC Zhang and Kalos (1991); Anderson et al. (1991), which requires a large density of walkers and does not scale well with system size. The structure of the SD allows the cancellation of the oneparticle instability against the Bose state to be carried out analytically in AFQMC.
iii.5 The AFQMC algorithm

Specify the initial state of each walker , and assign it a weight and phase .

If the walker weight is nonzero, compute its overlap with the trial wavefunction and force bias.

Sample an auxiliaryfield configuration and propagate the walker in imaginary time.

Compute the overlap between the updated walker and the trial wavefunction.

Update the walker’s weight and phase as .

In a constrainedpath calculation, ignore the phase and apply the projection.

Repeat steps 2 to 6 for all walkers. This forms one step of the random walk.

Periodically perform the branching and population control procedure.

Periodically reorthonormalize the orbitals of the walkers.

Repeat this process until an adequate number of measurements has been collected.
Iv Recent methodological advances
The development and application of AFQMC in molecular systems is an active area of research. Methods are maturing so that applications should start to increase rapidly. Recent methodological advances include the computation of observables, correlation functions, and geometry optimization; strategies for computing energy differences via correlated sampling; incorporation of functionalities of other quantum chemistry and electronic structure methods to extend the reach of AFQMC and improve computational efficiency. Concurrently, many other algorithmic advances have been made, including the computation of excitations and imaginarytime correlation Motta et al. (2014b); Vitali et al. (2016), the treatment of spinorbit coupling Rosenberg et al. (2017), the use of selfconsistent constraints Qin et al. (2016), etc, which to date have mostly been applied in lattice models of correlated electrons and which will not be covered here due to space limitations.
iv.1 Computing groundstate properties: the backpropagation algorithm
The importance sampling transformation provides a stochastic representation of the groundstate wavefunction, and gives the possibility of computing the mixed estimator of an observable as a weighted average
(45) 
Unless , the mixed estimator of is biased by the trial wavefunction used for importance sampling. In order to remove this bias, a backpropagation (BP) technique was proposed Zhang et al. (1997); Purwanto and Zhang (2004); Motta and Zhang (2017a) in the framework of AFQMC. The starting point of the BP algorithm is the observation that
(46) 
for large and . Inserting (15) into (46) yields Purwanto and Zhang (2004); Motta and Zhang (2017a)
(47) 
where we have introduced a shorthand for the collection of auxiliary fields along a path segment: