Ab initio computations of molecular systems by the auxiliary-field quantum Monte Carlo method

Ab initio computations of molecular systems by the auxiliary-field quantum Monte Carlo method

Mario Motta Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA    Shiwei Zhang shiwei@wm.edu Department of Physics, the College of William and Mary, Williamsburg, VA 23187, USA
July 3, 2019
Abstract

The auxiliary-field quantum Monte Carlo (AFQMC) method provides a computational framework for solving the time-independent 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 non-orthogonal Slater determinants, the method resembles an ensemble of density-functional theory (DFT) calculations in the presence of fluctuating external potentials. Its computational cost scales as a low-power of system size, similar to the corresponding independent-electron calculations. Highly efficient and intrinsically parallel, AFQMC is able to take full advantage of contemporary high-performance computing platforms and numerical libraries. In this review, we provide a self-contained 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 quantum-mechanical 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 many-electron wave function can be found by solving the time-independent Schrödinger equation for the Born-Oppenheimer 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 independent-electron approaches, based on the celebrated density-functional 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 DFT-based approaches has been exceptional. Their difficulties are also well-documented, especially in so-called strongly correlated systems, for example in many transition-metal 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 coupled-cluster 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 density-matrix renormalization group (DMRG) White (1992a); White and Martin (1999); Chan and Head-Gordon (2002); Olivares-Amaya et al. (2015); Chan et al. (2016) and full-configuration 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 auxiliary-field 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 non-perturbative approach that balances accuracy and computational scaling. Utilizing a field-theoretic framework and casting the projection of the many-body ground state or the finite-temperature partition function as a random walk in non-orthogonal 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 independent-particle methods such as DFT or Hartree-Fock (HF). AFQMC is naturally parallelizable and very well suited for large-scale 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 step-by-step algorithmic outline is then provided. Then in the next section, we describe several recent methodological advances with AFQMC, including frozen-core, downfolding and embedding approaches, back-propagation 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 inter-nuclear repulsion, are fermionic creation and annihilation operators associated to an orthonormal basis of one-electron 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 two-body 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 atom-centered 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 high-level quantum chemistry methods such as CCSD(T), this is advantageous for reaching the CBS. We typically use well-established bases (e.g. cc-pVxZ for light atoms Jr. (1989); Woon and Jr. (1993), or cc-pVnZ-PP 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 one-electron basis suitable for the problem. In addition to Gaussians, AFQMC has also been implemented using plane-wave 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 Kohn-Sham orbitals derived from a plane-wave DFT calculation as basis sets, which provides a kind of hybrid between the plane-wave and the MO approach.

ii.1 Properties of (non-orthogonal) 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 one-body Hamiltonians, including those of non-interacting systems and independent-electron Hamiltonians from DFT or mean-field 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 non-orthogonal 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 non-orthogonal SDs are summarized below. (Additional details can be found in the Appendix.)

Theorem 1.

The following properties hold for any two non-orthogonal SDs and , parametrized by the matrices and , respectively, and for any spin-independent one-particle operator

  1. Overlap:  .

  2. One-particle reduced density matrix:

    (5)
  3. Generalized Wick’s theorem Wick (1950); Balian and Brezin (1969):

    (6)
  4. Thouless theorem Thouless (1960, 1961): the state is a SD parametrized by the matrices

    (7)
  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 Gram-Schmidt procedure).

ii.2 Ground-state 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 imaginary-time Schrödinger equation

(8)

where is any initial wavefunction non-orthogonal to , and is the ground-state energy of . The formal solution of (8) is obtained applying the imaginary-time propagator to ,

(9)

The ground-state 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 imaginary-time propagation in Eq. (9) as it is to solve the time-independent Schrödinger equation. However, for an independent-electron method such as HF, where is replaced by

(10)

Eq. (7) guarantees that the imaginary-time 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 non-interacting system, when the 2-body part of the Hamiltonian vanishes.) Similarly, a DFT calculation can be cast in this way, with the projection in imaginary-time as an alternative for the non-linear minimization of the energy .

AFQMC maps the imaginary-time projection of a many-body 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 many-body correlation effects is recovered by random sampling around a typical “mean-field”. Alternatively, the mapping can be seen as a different and complementary approach to fixed-node diffusion Monte Carlo Reynolds et al. (1982b); Foulkes et al. (2001b), which simulates the imaginary-time 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, non-orthogonal 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.

Figure 1: (color online) Pictorial illustration of the AFQMC algorithm, and of its beyond-mean-field nature. Left: independent-electron methods (e.g., HF or DFT) provide an approximation of the ground state wave function through a deterministic evolution in a manifold (gray surface) of Slater determinants (black curves converging to ). Right: AFQMC represents the ground state as a stochastic linear combination of Slater determinants by mapping the electron-electron interaction onto a fluctuating external potential, and the imaginary-time evolution onto an open-ended random walk in (colored curves departing from ).

ii.2.1 The Hubbard-Stratonovich transformation

The first step to provide a stochastic representation of the imaginary-time projection is often a decomposition in short-time propagators,

(11)

Within AFQMC, one first expresses the Hamiltonian as

(12)

with , one-body operators, explicitly given in the next Section. The short-time propagator is then approximated by the use of a Trotter breakup Trotter (1959); Suzuki (1976)

(13)

and the Hubbard-Stratonovich transformation Hubbard (1959); Stratonovich (1958) (see Appendix),

(14)

Combining Eqs. (13) and (14) yields the following representation of the short-time propagator,

(15)

with and

(16)

The short-time propagator is thus written as an integral over normally distributed parameters , called auxiliary fields, of exponentials of one-body operators, . Eq. (15) maps the original interacting system onto an ensemble of non-interacting systems coupled with fluctuating external potentials. The integration over the auxiliary fields recovers the effect of the two-particle interaction. While for a completely general two-body Hamiltonian the number of auxiliary fields is Al-Saidi 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 real-valued orbitals are used as one-electron basis set, the interaction tensor can be reshaped into a real-valued, symmetric and positive-definite 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 next-iteration 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 non-negative. When , the Cauchy-Schwarz 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.

Figure 2: (color online) The ratio between Cholesky vectors and basis functions, for atoms He to Kr and dimers He to Kr at bondlength , using STO-6G and cc-pVxZ bases with x=2,3,4,5,6 and threshold . Asymptotically . Inset: dependence of the AFQMC energy on the threshold for Cl at experimental equilibrium bondlength, . yield energies within 1 m of converged value.

ii.2.3 Mean-field 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 ground-state average of each operator ,

(23)

where the modified one-body operators are , and the constant shift can be obtained, for example, using the trial wave function,

(24)

This is called mean-field background in the AFQMC literature. HS transformation applied to the modified one-body operators can lead to smaller fluctuations in free-projection calculations and smaller systematic errors when a constraint is applied Al-Saidi et al. (2006a); Purwanto et al. (2009b).

ii.3 The free-projection 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 ground-state energy:

(25)

where is a trial wavefunction, typically the SD from HF or DFT, or a linear combination of SDs from a complete-active space self-consistent 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 fixed-length (of -segments) path in the auxiliary-field space, for example, by the Metropolis Monte Carlo (MC) algorithm. This is the form in which the auxiliary-field approach was first proposed, both at finite-temperature 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, open-ended random walks in the space of SDs. For free-projection 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 open-ended 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 ground-state energy can be computed as weighted average of the local energy functional,

(27)

The free-projection AFQMC estimator of the ground-state energy is exact. The walkers will need to be periodically re-orthonormalized to maintain numerical precision. The only potential systematic errors come from the Trotter error of finite time-step, and population control bias. These can be easily removed by standard extrapolation procedures, as illustrated in Fig. 3. The Trotter error in free-projection is typically quadratic in . Population control bias is typically , and should be small in a free-projection 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 re-orthonormalization of the walkers, which are shared by free-projection and constrained calculations, will be discussed under Implementation issues.

Figure 3: (color online) Left: time step extrapolation of free-projection and phaseless AFQMC energies for HO (cc-pVDZ level, triangular geometry with , ), using a population of walkers. (inset: extrapolation of phaseless AFQMC energies vs. the inverse population size, illustrating population control bias at small . A time step of was used here.) Right: emergence of the phase problem in free-projection AFQMC calculations of ethane (STO-6G basis). (inset: comparison between free-projection and phaseless AFQMC for short projection time).

Iii Constrained AFQMC calculations

iii.1 The phase problem

The free-projection 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 long-range 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 free-projection 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 free-projection 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 free-projection 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 over-complete space of non-orthogonal 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 ground-state 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 half-filled repulsive Hubbard model Hirsch (1985), where particle-hole 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 ground-state 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 signal-to-noise 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 imaginary-time 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 hyper-surface . The resulting approximation defines the constrained-path AFQMC (cp-AFQMC) 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 Hubbard-Stratonovich transformations result in complex-valued 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 constrained-path approximation. The approach to real materials has been referred to in the literature as the phaseless or phase-free 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 auxiliary-field paths in the random walk, we will often refer to them as cp-AFQMC for simplicity. Below we will introduce the key ingredients that enable the formulation of the cp-AFQMC 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 variance-reduction 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 importance-sampling transformation for the most general case below. The case of real-valued , such as in the Hubbbard model, will follow as a special example. Recalling that the Hubbard-Stratonovich transformation is defined up to a complex-valued 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 importance-sampling transformations. Introducing a short-hand 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 auxiliary-fields 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 complex-valued operators, the random walk dynamics is invariant under the transformations , . As the projection time grows, the overlaps undergo a rotationally-invariant 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 back-propagation 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 real-valued, the projection reduces to the constrained-path AFQMC, because or and , respectively.

The resulting method, the phaseless AFQMC (ph-AFQMC), is the state-of-the-art technique for performing AFQMC calculations in real materials. A ph-AFQMC calculations differs from a fp-AFQMC 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 multi-determinant 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 trial-wavefunctions, 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 Sherman-Morrison-Woodbury 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 ground-state 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 on-the-fly 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 cp-AFQMC 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 round-off 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 lowest-energy orbital. Within AFQMC, the instability is eliminated by periodically orthogonalizing the orbitals with, e.g., a modified Gram-Schmidt (mGS) procedure, by sweeping through the orbitals and of with the following steps:

(44)

for , and similarly for the spin-down 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 one-particle instability against the Bose state to be carried out analytically in AFQMC.

iii.5 The AFQMC algorithm

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

  2. If the walker weight is non-zero, compute its overlap with the trial wavefunction and force bias.

  3. Sample an auxiliary-field configuration and propagate the walker in imaginary time.

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

  5. Update the walker’s weight and phase as .

  6. In a constrained-path calculation, ignore the phase and apply the projection.

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

  8. Periodically perform the branching and population control procedure.

  9. Periodically re-orthonormalize the orbitals of the walkers.

  10. 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 imaginary-time correlation Motta et al. (2014b); Vitali et al. (2016), the treatment of spin-orbit coupling Rosenberg et al. (2017), the use of self-consistent 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 ground-state properties: the backpropagation algorithm

The importance sampling transformation provides a stochastic representation of the ground-state 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 back-propagation (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: