# Benchmark Phaseless Auxiliary-Field Quantum Monte Carlo Method for Small
Molecules

Zoran Sukurma  University of Vienna, Faculty of Physics and Center for
Computational Materials Science, Kolingasse 14-16, A-1090 Vienna, Austria
University of Vienna, Faculty of Physics & Vienna Doctoral School in Physics,
Boltzmanngasse 5, A-1090 Vienna, Austria  Martin Schlipf  VASP Software GmbH,
Sensengasse 8, 1090 Vienna, Austria  Moritz Humer  University of Vienna,
Faculty of Physics and Center for Computational Materials Science, Kolingasse
14-16, A-1090 Vienna, Austria  University of Vienna, Faculty of Physics &
Vienna Doctoral School in Physics, Boltzmanngasse 5, A-1090 Vienna, Austria
Amir Taheridehkordi  University of Vienna, Faculty of Physics and Center for
Computational Materials Science, Kolingasse 14-16, A-1090 Vienna, Austria
Georg Kresse  University of Vienna, Faculty of Physics and Center for
Computational Materials Science, Kolingasse 14-16, A-1090 Vienna, Austria
VASP Software GmbH, Sensengasse 8, 1090 Vienna, Austria

###### Abstract

We report a scalable Fortran implementation of the phaseless auxiliary-field
quantum Monte Carlo (ph-AFQMC) and demonstrate its excellent performance and
beneficial scaling with respect to system size. Furthermore, we investigate
modifications of the phaseless approximation that can help to reduce the
overcorrelation problems common to the ph-AFQMC. We apply the method to the 26
molecules in the HEAT set, the benzene molecule, and water clusters. We
observe a mean absolute deviation of the total energy of
$1.15\text{\,}\mathrm{k}\mathrm{cal}\mathrm{/}\mathrm{mol}$ for the molecules
in the HEAT set; close to chemical accuracy. For the benzene molecule, the
modified algorithm despite using a single-Slater-determinant trial
wavefunction yields the same accuracy as the original phaseless scheme with
400 Slater determinants. Despite these improvements, we find systematic errors
for the CN, CO2, and O2 molecules that need to be addressed with more accurate
trial wavefunctions. For water clusters, we find that the ph-AFQMC yields
excellent binding energies that differ from CCSD(T) by typically less than
$0.5\text{\,}\mathrm{k}\mathrm{cal}\mathrm{/}\mathrm{mol}$.

††preprint: AIP/123-QED

##  I Introduction

Density functional theory (DFT) approximately solves the many-body Schrödinger
equation by mapping it to a single-particle problem and has become routine in
quantum chemistry and condensed matter physics.Kohn and Sham (1965); Kresse
and Furthmüller (1996); Kresse and Joubert (1999) Despite its simplicity and
the fact that it can describe systems with hundreds of atoms, DFT is not
always accurate enough to serve as a general-purpose solution. The stretched
H2 molecule is a well-known example of the limitations of DFT, as its results
are qualitatively wrong. Kim  _et al._ (2015) The key problem is the treatment
of the electron-electron correlation effects via approximate exchange-
correlation (xc) energy functionals.Cohen, Mori-Sánchez, and Yang (2008)
Perdew established the analogy of _climbing Jacob’s ladder_ for the effort of
developing more advanced xc functionals to reach chemical accuracy.Perdew and
Schmidt (2001) Higher rungs offer higher accuracy at the price of higher
computational costs and lower transferability. Novel xc potentials based on
machine learning Nagai, Akashi, and Sugino (2020); Pederson, Kalita, and Burke
(2022) attempt to address the accuracy issue but still struggle with lower
transferability. It remains to be seen whether machine learning functionals
will replace traditional ones in the coming years. Machine-learned xc
potentials as well as machine-learned force fields Behler and Parrinello
(2007); Chen  _et al._ require highly accurate reference data. In addition,
more accurate methods are also important to calibrate new methods or to access
strongly correlated materials where DFT generally performs poorly.

For this reason, a whole spectrum of correlation-consistent methods has been
developed over the last 50 years.Szabo and Ostlund (1996); Helgaker,
Jørgensen, and Olsen (2000) Full configuration interaction (FCI) provides an
exact solution within a given basis set Handy (1980) but suffers from
exponential scaling. Therefore, it is only usable for modest system sizes (up
to $10^{10}$ Slater determinants). A clever partitioning of the Fock space led
to a class of methods called selected CI Holmes, Tubman, and Umrigar (2016);
Sharma  _et al._ (2017); Dash  _et al._ (2018) and its stochastic counterpart
full configuration interaction quantum Monte Carlo (FCIQMC). Booth, Thom, and
Alavi (2009); Cleland, Booth, and Alavi (2010); Ghanem, Guther, and Alavi
(2020) These methods enabled the treatment of the Fock space with up to
$10^{30}$ Slater determinants while maintaining the accuracy of the
traditional FCI. However, they still scale weakly exponentially making them
difficult to apply to more realistic systems.

Besides CI methods, the most prominent methods are based on the coupled-
cluster expansion. The most popular among them, the "gold-standard" coupled-
cluster singles, doubles, and perturbative triples (CCSD(T)),Purvis and
Bartlett (1982); Bartlett and Musiał (2007) is extensively used to treat
molecules and is known to give excellent results for systems with small static
correlation effects. High-accuracy benchmark data for small molecules are
usually calculated using the CC expansion. For instance, a series of papers
focused on high-accuracy extrapolated ab-initio thermochemistry (HEAT) include
electron-electron correlation up to full quintuples in the coupled cluster
expansion.Tajti  _et al._ (2004); Bomble  _et al._ (2005); Harding  _et al._
(2008) For the molecules in the HEAT set, these results are considered to be
as good as the FCI results. More importantly, the HEAT studies also show that
CCSD(T) alone is not accurate enough to achieve chemical accuracy
($<1\text{kcal}/\text{mol}$) but approaches it in many relevant cases.
Furthermore, CCSD(T) scales adversely with the 7th power of the system size
and performs poorly in the presence of strong static correlation effects, e.g.
bond dissociation. These drawbacks emphasize the need for alternative methods.

The density-matrix renormalization group (DMRG) is a variational method widely
used to treat systems with strong correlation effects.White and Martin (1999);
Schollwöck (2005); Chan and Sharma (2011) It uses matrix-product states to
encode the locality in one of the spatial dimensions and to reduce the
exponential scaling of the CI expansion. For low-dimensional quantum systems
(quantum dots, molecules extended in one dimension), DMRG is the method of
choice due to its efficiency and accuracy. Although DMRG is applicable to
arbitrary systems, its limitation of about 30 active orbitals prevents it from
being a general-purpose tool for more generic and compact molecules or solids.

Quantum Monte Carlo methods, such as variational Monte Carlo McMillan (1965)
(VMC) and diffusion Monte Carlo Anderson (1976); Ceperley, Chester, and Kalos
(1977); Foulkes  _et al._ (2001) (DMC), exhibit low polynomial scaling with
system size and are nowadays routinely applied to systems with hundreds of
electrons. However, Nemec and co-workersNemec, Towler, and Needs (2010) showed
that achieving chemical accuracy is difficult with DMC. They reported a mean
absolute deviation of atomization energies of $3.2\text{kcal}/\text{mol}$ for
55 molecules in the G2 set.Curtiss  _et al._ (1991) Additionally, VMC and DMC
require accurate models for multi-electron wavefunctions including linear
combinations of Slater determinants, Jastrow factors, backflow wavefunctions,
and many others. Finally, recent work that combines Deep Learning and VMC has
emerged as a promising way to tackle the quantum many-body problem.Pfau  _et
al._ (2020); Spencer  _et al._ (2020); Gerard  _et al._ (2022)

In this work, we use the auxiliary-field quantum Monte Carlo (AFQMC) method.
AFQMC is a well-established projector Monte Carlo method successfully used in
various applications. Al-Saidi, Zhang, and Krakauer (2007); Purwanto, Zhang,
and Krakauer (2009); Purwanto  _et al._ (2011); Virgus  _et al._ (2012);
Purwanto, Zhang, and Krakauer (2015); Motta  _et al._ (2017); Lee, Malone, and
Reichman (2020); Lee, Malone, and Morales (2020); Rudshteyn  _et al._ (2022)
The original AFQMC was formulated as a path-integral method using the
Metropolis algorithm.Hirsch (1985); Sugiyama and Koonin (1986); White  _et
al._ (1989); Sorella  _et al._ (1989) It behaved well for systems without
severe fermionic sign problem. Since almost all general fermionic systems
suffer from the sign problem, the reformulation into an open-ended random walk
led to the constrained path (cp)-AFQMC, which was remarkably successful for
model Hamiltonians.Zhang, Carlson, and Gubernatis (1995, 1997) In order to
handle general many-body Hamiltonians, the phaseless (ph)-AFQMC was proposed
by Zhang and co-workers.Zhang and Krakauer (2003) The ph-AFQMC was
successfully applied to isolated molecules using Gaussian-type orbitals (GTOs)
and to solids using a plane-wave basis (PWs).Suewattana  _et al._ (2007) With
the increasing popularity of the method, a wide range of ph-AFQMC applications
emerged in recent years.Al-Saidi, Zhang, and Krakauer (2007); Purwanto, Zhang,
and Krakauer (2009); Purwanto  _et al._ (2011); Virgus  _et al._ (2012);
Purwanto, Zhang, and Krakauer (2015); Motta  _et al._ (2017); Lee, Malone, and
Reichman (2020); Lee, Malone, and Morales (2020); Rudshteyn  _et al._ (2022)
One reason for its popularity is that the ph-AFQMC scales quartic with the
system size. This scaling stems from the local energy evaluation and can be
further reduced using tensor contraction,Motta  _et al._ (2019); Lee and
Reichman (2020) density-fitting techniques Malone, Zhang, and Morales (2019);
Weber  _et al._ (2022) or plane-wave basis.Suewattana  _et al._ (2007) In
contrast to DMC, AFQMC yields energies fully consistent with other traditional
quantum chemistry approaches, allowing correlation energies to be directly
compared.

Controlling the fermionic phase problem is crucial for the accuracy of the
AFQMC method. In the ph-AFQMC, the _phaseless approximation_ imposes a
constraint on the walker weights,Zhang and Krakauer (2003); Zhang (2013);
Motta and Zhang (2018) analogous to the fixed-node approximation in
DMC.Anderson (1976) The phaseless approximation can introduce somewhat
uncontrolled errors that often lead to overestimated correlation energies. In
the ph-AFQMC, improving the trial wavefunctions systematically reduces the
phaseless approximation error. As the trial wavefunction approaches the exact
ground-state wavefunction, the ph-AFQMC energy approaches the exact ground-
state energy and the phaseless error disappears. Usually, the trial
wavefunction consists of the single Slater determinant formed from Hartree-
Fock (HF) or Kohn-Sham (KS) orbitals. They are particularly popular because
they are easily accessible from all quantum chemistry / solid-state physics
codes and avoid additional computational costs introduced by more complex
trial wavefunctions. Recent research focused on multi-determinant trial
wavefunctionsShi and Zhang (2013); Chang, Rubenstein, and Morales (2016);
Mahajan and Sharma (2021); Mahajan, Lee, and Sharma (2022); Huggins  _et al._
(2022). For example, Mahajan _et. al._ used $10^{4}$ Slater determinants while
increasing the total cost of the ph-AFQMC computation by only a factor of 3 in
comparison to the single Slater determinant case.

In this study, we adopt a simpler approach. We will carefully scrutinize
whether modifications to the weight update reduce the need to go beyond
single-determinant trial wavefunctions. With these modifications, the ph-AFQMC
can provide close to chemical accuracy for the HEAT set. Bomble _et. al._
Bomble  _et al._ (2005) provided highly accurate CCSDTQP molecular energies at
the double-zeta basis set (cc-pVDZ). In this work, we benchmark the accuracy
of the ph-AFQMC energies against this reference data set of small molecules.
Borda _et. al._ recently performed a similar study on the G1 test set.Landinez
Borda, Gomez, and Morales (2019)

The rest of the paper is structured as follows: In Sec. II, the ph-AFQMC
method is briefly reviewed along with the proposed modifications, followed by
details on the implementation in Sec. III. In Sec. IV, we present and discuss
different applications of ph-AFQMC. Finally, Sec. V concludes our findings and
possibilities for future developments.

##  II AFQMC Formalism

In this section, we briefly introduce the AFQMC formalism. For a more detailed
overview of the theory, we refer the interested reader to one of the
reviews.Zhang and Krakauer (2003); Zhang (2013); Motta and Zhang (2018)
Consider the full many-body Hamiltonian written in any orthonormal one-
particle basis given by

$\displaystyle\begin{split}\hat{H}&=\hat{H}_{1}\ +\
\frac{1}{2}\sum_{g}\hat{L}_{g}^{2}\\\
&=\sum_{pq}h_{pq}\;\hat{a}_{p}^{\vphantom{\dagger}\dagger}\hat{a}_{q}\ +\
\frac{1}{2}\sum_{g}\
\sum_{pq}L_{g,pq}\;\hat{a}_{p}^{\vphantom{\dagger}\dagger}\hat{a}_{q}\
\sum_{rs}L_{g,rs}\;\hat{a}_{r}^{\vphantom{\dagger}\dagger}\hat{a}_{s},\end{split}$
(1)

where $\hat{a}_{p}^{\dagger}$ and $\hat{a}_{q}$ are fermionic creation and
annihilation operators, respectively. The single-particle Hamiltonian matrix
elements $h_{pq}$ include all one-body terms. The two-body Hamiltonian is
written as a sum of squares of one-body operators $\hat{L}_{g}$. In the
Gaussian basis, we can obtain these operators by iterative Cholesky
decomposition Beebe and Linderberg (1977); Koch, Sánchez de Merás, and
Pedersen (2003) of electron repulsion integrals (ERIs). The indices $p$, $q$,
$r$, and $s$ go over $N$ basis functions and the index $g$ goes over $N_{g}$
Cholesky vectors. Typically, $N_{g}\approx 10N$. We neglect the spin indices
for simplicity.

The exact many-body ground-state wavefunction $\ket{\Phi_{0}}$ is extracted
from the long-time imaginary propagation of an initial state
$\ket{\Psi_{\mathrm{I}}}$

$\ket{\Phi_{0}}=\lim_{n\to\infty}\left[e^{-\tau\left(\hat{H}-E_{0}\right)}\right]^{n}\ket{\Psi_{\mathrm{I}}},$
(2)

where $\tau$ is the imaginary time step and $E_{0}$ is the estimated ground-
state energy. This equation is exact in the limit of infinitesimally small
time steps. In this work, the Hartree-Fock orbitals form the initial
wavefunction.

To treat the electron-electron interaction efficiently, the Hubbard-
Stratanovich (HS) transformation Hubbard (1959); Stratonovich (1957) is used

$e^{-\frac{\tau}{2}\sum_{g}\hat{L}_{g}^{2}}=\int\text{d}x^{N_{g}}\;p(x_{g})\;e^{i\sqrt{\tau}\sum_{g}x_{g}\hat{L}_{g}}+\mathcal{O}(\tau^{2})~{},$
(3)

where $p(x_{g})$ is the standard normal distribution and $x_{g}$ are random
numbers drawn from this distribution. In other words, the Hubbard-Stratanovich
transformation maps the actual system of interacting particles onto a system
of non-interacting particles coupled to random fields. Since ERIs are positive
definite, random fields are purely imaginary, i.e., one can define a complex-
valued effective Hamiltonian $\bar{h}$ given by

$\bar{h}_{pq}^{w}=h_{pq}-\frac{i}{\sqrt{\tau}}\sum_{g}x_{g}^{w}L_{g,pq},$ (4)

where the superscript $w$ emphasizes that unique random fields are drawn for
each walker. The ensemble of $N_{w}$ walkers represents the exact many-body
ground state, where each walker is represented by the following elements:
real-valued weight $W$, phase $\theta$, and single Slater determinant
$\ket{\Psi}$. Therefore, the exact many-body ground-state wavefunction can be
approximated as

$\ket{\Phi_{0}}=\Bigl{(}\sum_{kw}W_{k}^{w}e^{i\theta_{k}^{w}}\Bigr{)}^{-1}\sum_{kw}W_{k}^{w}e^{i\theta_{k}^{w}}\frac{\ket{\Psi_{k}^{w}}}{\braket{\Psi_{\mathrm{T}}}{\Psi_{k}^{w}}},$
(5)

where the index $k$ stands for time steps and the index $w$ for walkers.
$\ket{\Psi_{\mathrm{T}}}$ denotes the trial wavefunction. The walkers are
initialized as follows:

$\displaystyle W_{0}^{w}$ $\displaystyle=1;$ $\displaystyle\theta_{0}^{w}$
$\displaystyle=0;$ $\displaystyle\ket{\Psi_{0}^{w}}$
$\displaystyle=\ket{\Psi_{\mathrm{I}}}.$ (6)

The equations of motion for the walkers are

$\displaystyle\ket{\Psi_{k+1}^{w}}=e^{-\tau\bar{h}^{w}}\ket{\Psi_{k}^{w}},$
(7) $\displaystyle
W_{k+1}^{w}e^{i\theta_{k+1}^{w}}=W_{k}^{w}e^{i\theta_{k}^{w}}\;\frac{\braket{\Psi_{\mathrm{T}}}{\Psi_{k+1}^{w}}}{\braket{\Psi_{\mathrm{T}}}{\Psi_{k}^{w}}}.$
(8)

During the AFQMC propagation, this approximate ground-state wavefunction
grants access to the observables of the system. The simplest and probably the
most important observable—the total energy—is defined as the mixed expectation
value of the Hamiltonian

$E_{0}=\frac{\braket{\Psi_{\mathrm{T}}}{\hat{H}}{\Phi_{0}}}{\braket{\Psi_{\mathrm{T}}}{\Phi_{0}}}\approx\frac{\sum_{kw}W_{k}^{w}e^{i\theta_{k}^{w}}E_{\text{loc}}(\Psi_{k}^{w})}{\sum_{kw}W_{k}^{w}e^{i\theta_{k}^{w}}},$
(9)

where the local energy estimator $E_{\text{loc}}(\Psi^{w})$ is computed using
the generalized Wick’s theorem Balian and Brezin (1969)

$E_{\text{loc}}(\Psi^{w})=\frac{\braket{\Psi_{\mathrm{T}}}{\hat{H}}{\Psi^{w}}}{\braket{\Psi_{\mathrm{T}}}{\Psi^{w}}}\\\
=\sum_{pq}h_{pq}G_{pq}^{w}+\frac{1}{2}\sum_{g}\sum_{pqrs}L_{g,pq}L_{g,rs}(G_{pq}^{w}G_{rs}^{w}-G_{ps}^{w}G_{rq}^{w}).$
(10)

The one-body reduced density matrix $G_{pq}^{w}$ is defined as

$G_{pq}^{w}\equiv
G_{pq}(\Psi^{w})=\frac{\braket{\Psi_{\mathrm{T}}}{\hat{a}_{p}^{\vphantom{\dagger}\dagger}\hat{a}_{q}}{\Psi^{w}}}{\braket{\Psi_{\mathrm{T}}}{\Psi^{w}}}=\left[\Psi^{w}(\Psi_{\mathrm{T}}^{\dagger}\Psi^{w})^{-1}\Psi_{\mathrm{T}}^{\dagger}\right]_{qp}.$
(11)

Here, $\Psi_{\mathrm{T}}$ and $\Psi^{w}$ represent the $N_{\mathrm{e}}$
occupied orbitals in the $N$ basis functions. Equations (7, 8, 9, 10) define
the free-propagation (fp)-AFQMC.

###  II.1 Mean-Field Subtraction

A decrease in the magnitude of the diagonal components of the Cholesky vectors
$L_{g,pq}$ leads to smaller fluctuations in AFQMC and thus better statistics.
To this end, one introduces a shift

$\bar{L}_{g}=\braket{\hat{L}_{g}}_{\mathrm{T}}=\frac{\braket{\Psi_{\mathrm{T}}}{\hat{L}_{g}}{\Psi_{\mathrm{T}}}}{\braket{\Psi_{\mathrm{T}}}{\Psi_{\mathrm{T}}}}$
(12)

to the two-body Hamiltonian

$\hat{H}_{2}=\frac{1}{2}\sum_{g}\left(\hat{L}_{g}-\bar{L}_{g}\right)^{2}+\sum_{g}\bar{L}_{g}\hat{L}_{g}-\frac{1}{2}\sum_{g}\bar{L}_{g}^{2}~{}.$
(13)

Here, we identify the second term as the classical Hartree potential
corresponding to the trial wavefunction and add it to the one-body
Hamiltonian. The last term is the Hartree energy corresponding to the trial
wavefunction. The mean-field subtraction incurs no computational cost for the
AFQMC propagation.

###  II.2 Importance Sampling

Introducing an importance function in AFQMC shifts the Gaussian probability
density in Eq. (3). This shift is referred to as the "force bias". Choosing
the force bias as

$f_{g}^{w}=-i\sqrt{\tau}\
\frac{\braket{\Psi_{\mathrm{T}}}{\hat{L}_{g}-\bar{L}_{g}}{\Psi^{w}}}{\braket{\Psi_{\mathrm{T}}}{\Psi^{w}}}.$
(14)

minimizes the statistical fluctuations in the random fields.Zhang (2013) After
applying the mean-field subtraction and importance sampling, the effective
Hamiltonian (Eq. (4)) changes to

$\bar{h}_{pq}^{w}=h_{pq}-\frac{i}{\sqrt{\tau}}\sum_{g}(x_{g}^{w}-f_{g}^{w})(L_{g,pq}-\bar{L}_{g}\delta_{pq}).$
(15)

The force bias is a walker-dependent quantity and must be calculated in each
time step for each walker. The force bias is crucial to successful AFQMC
simulations, however, computing it entails considerable computational effort
for the AFQMC procedure and may even be the computational bottleneck.Motta and
Zhang (2018)

###  II.3 Phaseless Approximation

The fp-AFQMC is formally exact but suffers from the fermionic phase problem.
The problem is observed in the expression for the energy evaluation (Eq. (9)).
After the equilibration time, during which excited states contained in
$\ket{\Psi_{\mathrm{I}}}$ decay exponentially, walkers will populate the
entire complex plane uniformly. Therefore, the denominator in Eq. (9) vanishes
and the whole expression becomes ill-defined. To circumvent the fermionic
phase problem, the phaseless approximation is introduced.

Since the effective Hamiltonian given in Eq. (15) is complex-valued, the
orbitals become complex-valued as well. In DMC, $\ket{\Phi}$ and $-\ket{\Phi}$
are both valid solutions of the Schrödinger equation. Similarly in AFQMC, if
$\ket{\Phi}$ is the valid solution, then $e^{i\theta}\ket{\Phi}$ is a valid
solution for any $\theta\in[0,2\pi)$. To prevent the exponential growth of the
noise in AFQMC simulation, Zhang _et. al._ Zhang and Krakauer (2003)
introduced a phaseless approximation

$\displaystyle\begin{aligned}
W_{k+1}^{w}=W_{k}^{w}\>\left|\frac{\braket{\Psi_{\mathrm{T}}}{\Psi_{k+1}^{w}}}{\braket{\Psi_{\mathrm{T}}}{\Psi_{k}^{w}}}I^{w}\right|\
\text{max}\left(0,\>\text{cos}(\Delta\theta)\right),\end{aligned}$ (16)
$\displaystyle\theta_{k}^{w}=0,$ (17)

where the importance sampling factor $I^{w}$ and the phase change
$\Delta\theta$ are defined as

$\displaystyle
I^{w}=\exp\Bigl{[}\sum_{g}x_{g}^{w}f_{g}^{w}-\frac{1}{2}f_{g}^{w}f_{g}^{w}\Bigr{]},$
(18) $\displaystyle\Delta\theta=\mathfrak{Im}\
\text{ln}\frac{\braket{\Psi_{\mathrm{T}}}{\Psi_{k+1}^{w}}}{\braket{\Psi_{\mathrm{T}}}{\Psi_{k}^{w}}}\approx\mathfrak{Im}\sum_{g}x_{g}^{w}f_{g}^{w}.$
(19)

The update equation (16) can be also rewritten as

$W_{k+1}^{w}=W_{k}^{w}\>\text{max}\left(0,\
\mathfrak{Re}\frac{\braket{\Psi_{\mathrm{T}}}{\Psi_{k+1}^{w}}}{\braket{\Psi_{\mathrm{T}}}{\Psi_{k}^{w}}}\right)\>|I^{w}|.$
(20)

Zhang _et. al._ Zhang and Krakauer (2003); Zhang (2013) state that different
approximations yield the same expectation values and slightly different
standard deviations.

In this work, we investigate this statement more quantitatively. To this end,
we modify the expression for the total energy evaluation (Eq. (9))

$E_{0}=\frac{\sum_{kw}W_{k}^{w}\cos(\theta_{k}^{w})E_{\text{loc}}(\Psi_{k}^{w})}{\sum_{kw}W_{k}^{w}\cos(\theta_{k}^{w})}$
(21)

and the update procedure for the walker weights

$\displaystyle
W_{k+1}^{w}e^{i\theta_{k+1}^{w}}=W_{k}^{w}e^{i\theta_{k}^{w}}\>\frac{\braket{\Psi_{\mathrm{T}}}{\Psi_{k+1}^{w}}}{\braket{\Psi_{\mathrm{T}}}{\Psi_{k}^{w}}}\mathfrak{Re}I^{w}~{}.$
(22)

We refer to this new approach as ph∗-AFQMC throughout this paper. In contrast
to the original approach, there are two key differences: First, we retain the
complex nature of the walker weights and use the real part for the population
control and evaluation of physical observables, such as the energy. Walkers
are killed explicitly if $|\theta_{k}^{w}|\geq\frac{\pi}{2}$. Secondly, we
include the real part of the importance weight instead of the absolute value.
We suggest that only the real part is meaningful and using the absolute value
increases the weights with the imaginary part. This may contribute to the
systematic over-correlation of ph-AFQMC. We will show in Sec. IV that
ph∗-AFQMC systematically yields correlation energies smaller in absolute
value.

##  III Implementation Details

In this section, we present our implementation of AFQMC in QMCFort.111QMCFort
is available from the authors upon a reasonable request. QMCFort is written in
Fortran and utilizes (threaded) BLAS/LAPACK for fast linear algebra and MPI
and OpenMP for parallelization. In addition to AFQMC, QMCFort offers
implementations of restricted, restricted open-shell, and unrestricted
Hartree-Fock (RHF, ROHF and UHF) and MP2. The electron repulsion integrals
(ERIs) are calculated using the McMuchie-Davidson schemeHelgaker, Jørgensen,
and Olsen (2000) and decomposed into Cholesky vectors. Therefore, QMCFort can
act as a standalone tool to setup and run AFQMC simulations. However, only the
AFQMC part is highly optimized for large-scale calculations. For this reason,
we also interface QMCFort with PySCFSun  _et al._ (2018) for isolated systems
and with VASPKresse and Furthmüller (1996) for extended systems. The
implementation of the computationally intensive parts of the AFQMC procedure
and the interface between QMCFort and VASP are detailed in the remainder of
this section.

###  III.1 Effective Hamiltonian

At each time step and for each walker, we build the effective Hamiltonian
according to Eq. (15). The computationally demanding part is the convolution
of the shifted random fields with Cholesky vectors
$\sum_{g}(x_{g}^{w}-f_{g}^{w})\cdot L_{g,pq}$. Treating the orbital indices
$p$ and $q$ as a combined index, we map this contraction to matrix-matrix
multiplication with a cubic scaling ($\sim N^{2}N_{g}N_{w}$) in system size.
Typically, we treat in the range of 10 to 100 walkers per MPI rank and compute
this contraction for all of them simultaneously. Although the exchange energy
evaluation (Sec.III.4) scales worse, the creation of the effective Hamiltonian
is often the most expensive operation in AFQMC because it is required at each
time step.

###  III.2 AFQMC Propagation

The effective Hamiltonian propagates the orbitals according to Eq. (7). The
matrix exponentials in the Trotter-Suzuki propagator Trotter (1959); Suziki
(1976) are approximated by the sixth-order Taylor expansion. This requires
several applications of the effective Hamiltonian to the orbitals, i.e.
matrix-matrix multiplications of the form

$\Psi^{w}_{pi}\to\bar{\Psi}^{w}_{pi}=\sum_{q}\bar{h}_{pq}^{w}\Psi^{w}_{qi}~{}.$
(23)

These scale as $N^{2}N_{e}$.

###  III.3 Force Bias

The force bias evaluation (see Eq. (14)) is equivalent to the evaluation of
the Hartree potential and energy. It scales cubically with respect to the
system size, i.e. $NN_{g}N_{e}N_{w}$. Similar to the effective Hamiltonian
(Sec. III.1), the force bias is computed for all walkers simultaneously. To
make it more transparent, we will first define the contracted Cholesky vectors

$\mathcal{L}_{g,pi}=\sum_{q}L_{g,pq}\Psi_{\text{T},qi},$ (24)

and the force bias is then

$f_{g}^{w}=\sum_{pi}\mathcal{L}_{g,pi}\Psi_{pi}^{w}~{}.$ (25)

Here, we can treat $\\{p,i\\}$ as a single index and evaluate the contraction
with a matrix-matrix multiplication. We use the force bias to compute the
Hartree energy for the walker $\Psi^{w}$

$E_{\text{H}}(\Psi^{w})=\frac{1}{2}\sum_{g}f_{g}^{w}f_{g}^{w}.$ (26)

###  III.4 Exchange Energy

Of all individual contributions to AFQMC, the computation of the exchange
energy incurs the largest computational cost. Although the naive
implementation of Eq. (10) scales as $N^{3}N_{g}N_{w}$, we simplify it by
computing intermediate arrays $\alpha_{g,ij}^{w}$

$\alpha_{g,ij}^{w}=\sum_{p}\mathcal{L}_{g,pi}\Psi_{pj}^{w}~{}.$ (27)

The exchange energy of the walker $\Psi^{w}$ is then given by

$E_{\text{x}}(\Psi^{w})=-\frac{1}{2}\sum_{g}\sum_{ij}\alpha_{g,ij}^{w}\alpha_{g,ji}^{w}~{}.$
(28)

The quartic-scaling ($NN_{g}N_{e}^{2}N_{w}$) construction of the intermediate
arrays $\alpha_{g,ij}^{w}$ is also mapped to a matrix-matrix multiplication.
Multiple walkers are treated simultaneously in a block of approximately 5-10
walkers for optimal performance. One can evaluate the exchange energy only say
every 10th time step. This has no significant impact on the statistics because
one needs to sample independent local energies so the alternative is block
averaging. Computing the exchange energy from the $\alpha$ tensors scales as
$N_{g}N_{e}^{2}$ but the computational effort is negligible in comparison to
other tasks discussed in this section.

###  III.5 Other tasks

Besides the operations mentioned above, there are a few additional parts of
the AFQMC procedure that should be mentioned. Since the AFQMC propagator is
non-unitary, we periodically reorthogonalize the walkers by a QR decomposition
of the walker matrices $\Psi^{w}$. The walker population is controlled to
avoid too small or too large weights. We opted for the comb methodCalandra
Buonaura and Sorella (1998), because it keeps the total number of walkers
constant.

###  III.6 Performance

Table 1:  Number of operations per walker and time step and asymptotic scaling of the computationally intensive parts of the AFQMC procedure. $N$ represents the number of basis functions, $N_{g}$ is the number of Cholesky vectors, and $N_{e}$ is the number of occupied states in the system. Task | # of operations | asymptotic scaling
---|---|---
effective Hamiltonian | $N^{2}N_{g}$ | $\mathcal{O}(N^{3})$
AFQMC propagation | $N^{2}N_{e}$ | $\mathcal{O}(N^{3})$
force bias | $NN_{g}N_{e}$ | $\mathcal{O}(N^{3})$
exchange energy | $NN_{g}N_{e}^{2}+N_{g}N_{e}^{2}$ | $\mathcal{O}(N^{4})$

In Table 1, we summarize the number of operations corresponding to the
computationally demanding tasks. Assuming that the number of walkers is
independent of the system size, we list the resulting asymptotic large system-
size scaling for these tasks. We gauge the asymptotic behavior for water
clusters of various sizes containing up to 5 water molecules. The calculations
were performed on a dual-socket Intel Skylake Platinum 8174 with 24 cores each
and a nominal base frequency of $3.1\text{\,}\mathrm{GHz}$. We average the
timings over 200 AFQMC steps, with each step processing 4 800 walkers on 48
MPI processes.

Figure 1: Elapsed time per AFQMC step of the computationally intensive
operations measured on water clusters of different sizes. AFQMC calculations
were performed on a dual-socket Intel Skylake Platinum 8174 with 4 800
walkers, and 48 MPI processes. Timings are averaged over 200 AFQMC steps.
Figure 2: Performance per CPU core of the computationally intensive operations
in the AFQMC procedure measured on water clusters of different sizes. Tuples
on the $x$-axis denote the number of electrons and orbitals. Executing the
AFQMC code on one CPU node (dual-socket Intel Skylake Platinum 8174) delivers
1–1.5 TFLOPS (30 GFLOPS per core, green line) on average.

Fig. 1 shows the elapsed times per AFQMC step of the computationally demanding
tasks. To estimate the effective scaling exponent for each task, we fit the
power function $f(x)=Cx^{\beta}$ through the data points. As predicted, the
construction of the effective Hamiltonian, AFQMC propagation, and force bias
calculation scale cubically. The measured exponent $\beta=3.1$ for the
exchange energy does not exactly match the expected exponent $\beta=4$. With
increasing system size, we observe an increasing per-core performance of the
exchange evaluation (Fig. 2). This increase effectively reduces the scaling
exponent. Still, the exchange evaluation is the least performant of the
relevant computational tasks. As depicted in Fig. 2, most other tasks run
efficiently with the performance of 20–35 GFLOPS/core. A notable exception is
the AFQMC propagation that reaches a performance above 50 GFLOPS/core for
large systems, very close to peak performance.

###  III.7 Interface with VASP

VASP provides a set of one-electron mean-field orbitals
$\\{\phi_{p}(\textbf{r})\\}$. From these, we compute the single-particle
Hamiltonian matrix elements $h_{pq}$ and Cholesky vectors $L_{g,pq}$ (see
Appendix for more details). We export the matrix elements and Cholesky vectors
to a file and read them with QMCFort. Then, the AFQMC calculation for extended
systems is equivalent to calculations on isolated molecules using the Gaussian
basis set.

Figure 3:  Convergence behavior of the CCSD(T) and ph-AFQMC correlation
energies as a function of the number of canonical orbitals for diamond. Linear
extrapolation to the complete basis set is performed using the first three
points ($1024^{-1}$, $512^{-1}$, and $320^{-1}$). The calculations were
performed using a 2x2x2 diamond supercell with an energy cutoff of
$700\text{eV}$.

As an example, we compute the 2x2x2 supercell of diamond. The convergence of
the AFQMC and CCSD(T) correlation energies with respect to the number of
canonical orbitals is studied. CCSD(T) energies are computed using Cc4s
code.Hummel, Tsatsoulis, and Grüneis (2017) A lattice constant of $a=3.567$ Å
is used with an energy cutoff of $E_{\text{cut}}=700\ \text{eV}$. In addition,
we use an energy cutoff $E_{\text{cut}}^{\chi}=500\ \text{eV}$ for the
truncation of the Coulomb kernel. VASP relies on PAW potentialsBlöchl (1994);
Kresse and Joubert (1999) and the potential referred to as C_GW with valence
$2s2p$ was used. The core radius of this potential is $r_{\text{core}}=1.5\
\text{a.u.}$; the local potential is equivalent to the d-pseudopotential and 4
projectors are used (two $s$-projectors with $r_{\text{cut}}=1.2\
\text{a.u.}$, and two $p$-projectors with $r_{\text{cut}}=1.5\ \text{a.u.}$).

Figure 3 shows the convergence of the correlation energy as a function of the
number of canonical orbitals for ph-AFQMC and CCSD(T). The number of orbitals
is varied from the small basis set with 64 electrons in 192 orbitals (64e,
192o) to a relatively large basis of 64 electrons in 1024 orbitals (64e,
1024o). A $1/N$ behavior is fitted through the data points. The ph-AFQMC and
CCSD(T) lines are almost indistinguishable suggesting that the methods yield
consistent energies for periodic systems. Further details on AFQMC
calculations for extended systems will be reported in a forthcoming
publication. Taheridehkordi (2023)

To compare the computational cost of the CCSD(T) and ph-AFQMC, we conducted
calculations using 256 and 1024 bands. For the 256-band calculation, the
CCSD(T) calculation took 20 minutes on a single node (2x AMD Epyc (Milan)),
while the ph-AFQMC calculation took 11 hours. However, for the 1024-band
calculation, CCSD(T) required 24 hours on 4 nodes (2x AMD Epyc (Milan)), while
ph-AFQMC calculation took 48 hours to complete. The timings we obtained for
the CCSD(T) and ph-AFQMC calculations demonstrate the favorable scaling of the
AFQMC method, which becomes increasingly advantageous as the number of
orbitals and the complexity of the system increase.

##  IV Results and Discussion

Table 2:  Total molecular energies (in Hartree units) for the HEAT setTajti  _et al._ (2004) molecules at different levels of theory: ph-AFQMCZhang and Krakauer (2003), ph∗-AFQMC, CCSD(T), CCSDTQPBomble  _et al._ (2005) and HCIHolmes, Tubman, and Umrigar (2016); Sharma  _et al._ (2017). All energies are calculated using the cc-pVDZ basis set. The last three rows of the table depict the root-mean-square deviation (RMSD), mean signed deviation (MSD), and mean absolute deviation (MAD). Molecule | ph-AFQMC | ph∗-AFQMC | CCSD(T) | CCSDTQP | HCI
---|---|---|---|---|---
H2 | $-1.16363(2)$ | $-1.16338(3)$ | $-1.163426$ | $-1.163426$ | $-1.163426$
CH | $-38.37764(6)$ | $-38.37705(7)$ | $-38.379655$ | $-38.380241$ |
CH2 | $-39.04181(5)$ | $-39.04150(6)$ | $-39.041196$ | $-39.04165$ |
NH | $-55.09087(5)$ | $-55.09077(7)$ | $-55.097448$ | $-55.0917$ |
CH3 | $-39.71651(6)$ | $-39.71592(7)$ | $-39.715543$ | $-39.71607$ |
NH2 | $-55.73265(6)$ | $-55.73222(8)$ | $-55.732506$ | $-55.733076$ |
OH | $-75.55854(6)$ | $-75.55837(8)$ | $-75.559233$ | $-75.559689$ |
HF | $-100.22933(7)$ | $-100.2290(1)$ | $-100.228131$ | $-100.228622$ |
H2O | $-76.24249(8)$ | $-76.24167(9)$ | $-76.241018$ | $-76.241649$ | $-76.241637$
NH3 | $-56.40334(7)$ | $-56.4024(1)$ | $-56.401913$ | $-56.402517$ |
C2H | $-76.3999(2)$ | $-76.3986(2)$ | $-76.398558$ | $-76.401217$ |
CN | $-92.4997(2)$ | $-92.4973(2)$ | $-92.488695$ | $-92.49276$ | $-92.492788$
C2H2 | $-77.1118(2)$ | $-77.1092(2)$ | $-77.109249$ | $-77.110678$ |
CO | $-113.0594(2)$ | $-113.0567(2)$ | $-113.054431$ | $-113.055892$ |
HCN | $-93.1912(2)$ | $-93.1881(2)$ | $-93.188321$ | $-93.18991$ |
N2 | $-109.2781(2)$ | $-109.2750(2)$ | $-109.275298$ | $-109.277012$ | $-109.277005$
HCO | $-113.5778(2)$ | $-113.5757(2)$ | $-113.575706$ | $-113.577384$ |
CF | $-137.4764(1)$ | $-137.4754(2)$ | $-137.474848$ | $-137.476019$ |
NO | $-129.5970(2)$ | $-129.5945(2)$ | $-129.597778$ | $-129.599737$ |
HNO | $-130.1736(2)$ | $-130.1697(2)$ | $-130.170989$ | $-130.172906$ |
O2 | $-149.9792(2)$ | $-149.9787(2)$ | $-149.985684$ | $-149.987773$ |
HO2 | $-150.5606(2)$ | $-150.5603(2)$ | $-150.558481$ | $-150.56038$ |
OF | $-174.5004(2)$ | $-174.4998(2)$ | $-174.497924$ | $-174.50009$ |
H2O2 | $-151.1963(2)$ | $-151.1940(2)$ | $-151.19363$ | $-151.195266$ |
F2 | $-199.0968(2)$ | $-199.0965(2)$ | $-199.097448$ | $-199.099328$ |
CO2 | $-188.1561(2)$ | $-188.1525(2)$ | $-188.147429$ |  | $-188.149551$
RMSD (in $\mathrm{m\text{ $E_{\textup{ h }}$ }}$) | $2.89$ | $2.73$ | $1.65$ |  |
MSD (in $\mathrm{m\text{ $E_{\textup{ h }}$ }}$) | $-0.38$ | $1.17$ | $1.39$ |  |
MAD (in $\mathrm{m\text{ $E_{\textup{ h }}$ }}$) | $1.85$ | $1.84$ | $1.39$ |  |
Figure 4: CCSD(T), ph-AFQMC and ph∗-AFQMC total-energy differences $\Delta E$
relative to the reference CCSDTQP energies for the HEAT set. Calculations use
the cc-pVDZ basis set and a frozen-core approximation including only the last
occupied shell. Trial wavefunctions are single RHF/UHF Slater determinants.
The solid black line represents the chemical accuracy range ($\Delta
E<$1\text{\,}\mathrm{k}\mathrm{cal}\mathrm{/}\mathrm{mol}$$)

In the following four subsections, we report (i) the total molecular energies
for 26 molecules from the HEAT set, (ii) the AFQMC energy of the benzene
molecule using the cc-pVDZ basis set, and (iii) the AFQMC binding energies of
the water clusters.

###  IV.1 Heat Molecules

We calculated total energies for 26 molecules from the HEAT set using the
frozen-core approximation and geometries from Ref. Bomble  _et al._ , 2005
(HEAT). All calculations employed the double-zeta basis set (cc-pVDZ). We used
the QMCFort code to calculate Cholesky vectors, the reference Hartree-Fock
(RHF/UHF) wavefunctions, and all AFQMC energies. We used the PySCF package Sun
_et al._ (2018) to calculate the CCSD(T) energies. The Dice code Holmes,
Tubman, and Umrigar (2016); Sharma  _et al._ (2017) was used to calculate the
heat-bath CI energies. To eliminate possible sources of systematic errors in
the AFQMC, we (i) truncated Cholesky vectors at a conservative threshold of
$10^{-8}$ to ensure a reasonably accurate representation of the molecular
ERIs, (ii) chose a relatively small time step of $0.002\text{\,}\mathrm{\text{
$E_{\textup{ h }}$ }}^{-1}$ to exclude significant time-step errors, and (iii)
equilibrated the system for 40 000 time steps. To ensure good statistics, we
propagated 6 000 walkers until the standard error of the mean of the predicted
molecular energies dropped below $0.2\text{\,}\mathrm{m\text{ $E_{\textup{ h
}}$ }}$.

Table 2 compares the ph-AFQMC and ph∗-AFQMC molecular energies to the "gold-
standard" CCSD(T) and the more accurate coupled-cluster expansion including
variational triple, quadruple and pentuple excitations (CCSDTQP). We use the
latter as reference values since they are practically converged against the
FCI limit. The CCSDTQP value for the CO2 molecule is missing in Ref. Bomble
_et al._ , 2005, so we use the result of the heat-bath configuration
interaction (HCI).Holmes, Tubman, and Umrigar (2016) To verify the accuracy of
the HCI method, we selected a few other molecules (H2, H2O, CN, N2) and found
a maximal difference to the CCSDTQP values of $0.01\text{\,}\mathrm{m\text{
$E_{\textup{ h }}$ }}$, well below the statistical AFQMC errors. Comparing to
these reference values, we find a similar precision for the ph-AFQMC and
ph∗-AFQMC method with a root-mean-square deviation (RMSD) of
$2.89\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$ and
$2.73\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$, respectively. For
context, the CCSD(T) exhibits a RMSD of $1.65\text{\,}\mathrm{m\text{
$E_{\textup{ h }}$ }}$ for these molecules. The mean absolute deviations (MAD)
of $1.39\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$ (CCSD(T)),
$1.85\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$ (ph-AFQMC) and
$1.84\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$ (ph∗-AFQMC) show a
smaller difference between the methods indicating that some outliers taint the
overall performance of AFQMC. The mean signed deviation (MSD) of
$1.39\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$ (CCSD(T)) and
$1.17\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$ (ph∗-AFQMC) indicate
under-correlation in contrast to ph-AFQMC with an average over-correlation of
$-0.38\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$. The modified approach
(ph∗-AFMQC) is also closer to CCSD(T) with an RMSD of
$2.66\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$ compared to ph-AFQMC with
an RMSD of $3.67\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$.

Next, we present a more detailed analysis based on the graphical
representation of the results in Fig. 4. Let us focus on the more-widely-used
CCSD(T) method first. Fig. 4 illustrates that CCSD(T) systematically under-
correlates compared to the reference values. The most frequent deviation lies
around $1\text{\,}\mathrm{k}\mathrm{cal}\mathrm{/}\mathrm{mol}$ (i.e.
$1.59\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$). The worst case is the
CN molecule with a deviation of $4.1\text{\,}\mathrm{m\text{ $E_{\textup{ h
}}$ }}$. The reason for the larger deviation is probably the large spin
contamination in the UHF wavefunction.Bomble  _et al._ (2005)

The ph-AFQMC and ph∗-AFQMC results show similar trends. The O2 molecule is the
worst case for AFQMC with an undercorrelation of $8.6\text{\,}\mathrm{m\text{
$E_{\textup{ h }}$ }}$. Solving this problem requires a multi-determinant
trial wavefunction, for example from a complete active space self-consistent
field (CASSCF) calculation in the open-shell subspace. A similar problem in
open-shell atoms is solved by using multi-determinant CASSCF trial
wavefunctions or by using symmetry restoration techniques.Motta and Zhang
(2018) Excluding O2 from the statistics reduces the RMSD to 2.41 and
$2.11\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$ for the ph-AFQMC and
ph∗-AFQMC method, respectively.

In contrast to the O2 molecule, the ph-AFQMC values for CN, CO, and CO2
molecules show considerable over-correlation (up to $7\text{\,}\mathrm{m\text{
$E_{\textup{ h }}$ }}$ for the CN molecule). The origin of these deviations is
the fixed-node error. Our modified approach (ph∗-AFQMC) systematically
increases the energies compared to ph-AFQMC. It is noteworthy that the
differences are larger in the cases where the ph-AFQMC shows larger over-
correlation. The fixed-node errors are therefore significantly reduced for the
CN, CO, and CO2 molecules but the residual errors are still sizable. Better
trial wavefunctions could further reduce the fixed-node errors.

Borda _et. al._ Landinez Borda, Gomez, and Morales (2019) performed similar
benchmarks using the G1 test set that partially overlaps with the HEAT set.
They also used cc-pVDZ basis sets with the frozen-core approximation and
compared AFQMC total energies calculated with the QMCPACK package to the
CCSDTQ energies. While we used UHF trial wavefunctions, they used ROHF ones
for the AFQMC calculations. Fig. 5 shows the deviation of ph-AFQMC from the
reference energies. Overall, most energies agree within 1 m$E_{\rm h}$ between
the two ph-AFQMC calculations. For the remaining cases, our results appear to
be closer to the reference result except for the CN molecule. Possible reasons
for these differences include: (i) smaller statistical and systematic errors
in the present work, (ii) different molecular geometries, and (iii) different
trial wavefunctions. The latter matters particularly for the CN molecule
because the large spin contamination in the UHF wavefunction makes the ROHF
wavefunction a better choice for the trial wavefunction.

Figure 5: Absolute total energy difference $|\Delta E|$ between ph-AFMQC
(QMCPACK) and CCSDTQ (see Ref. Landinez Borda, Gomez, and Morales, 2019) as
well as ph-AFQMC (QMCFort) and CCSDTQP for the overlapping molecules in the G1
and HEAT set.

Next, we study the basis-set impact by comparing the total energy differences
between ph∗-AFQMC and CCSD(T) for the cc-pVDZ basis set and the roughly six-
times larger aug-cc-pVQZ one. The latter includes core electrons overcoming
the frozen-core approximation of the smaller basis. Fig. 6 shows remarkably
similar energy differences for both basis sets. Three notable exceptions are
F2, O2, and OH where the deviation is slightly larger than
$1\text{\,}\mathrm{k}\mathrm{cal}\mathrm{/}\mathrm{mol}$. This result is
important for two reasons: First, the comparison between ph∗-AFQMC and CCSD(T)
is sufficiently accurate using a cc-pVDZ basis set. Hardly any new conclusion
could be drawn from the complete basis-set limit. Secondly, the ph∗-AFQMC
precisely describes core-electron effects at the accuracy of the CCSD(T)
method.

Figure 6: Total-energy difference $\Delta E$ between ph∗-AFQMC and CCSD(T) for
the HEAT set. The cc-pVDZ calculations employ the frozen-core approximation;
the aug-cc-pVQZ ones include all electrons. Trial wavefunctions are single RHF
and UHF Slater determinants for closed- and open- shell molecules,
respectively. The solid black line represents chemical accuracy ($\Delta
E<$1\text{\,}\mathrm{k}\mathrm{cal}\mathrm{/}\mathrm{mol}$$).

###  IV.2 Benzene Molecule

Figure 7: Correlation energies for the benzene molecule using the cc-pVDZ
basis set and the frozen-core approximation. ph-AFQMC(RHF), ph∗-AFQMC(RHF)
correlation energies are obtained using QMCFort code, while the CCSD(T) energy
is calculated using PySCF code.Sun  _et al._ (2018) Other values are taken
from the Refs. Eriksen  _et al._ , 2020; Lee, Malone, and Reichman, 2020. The
red line represents the reference correlation energy of
$-863.0\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$, while the dashed black
lines represent the range of the chemical accuracy
($1\text{\,}\mathrm{k}\mathrm{cal}\mathrm{/}\mathrm{mol}$). All energy values
are given in $\mathrm{m\text{ $E_{\textup{ h }}$ }}$.

The ground state of the benzene molecule with double-zeta basis set (cc-pVDZ)
is an interesting system to benchmark state-of-the-art correlation-consistent
methods because it is among the largest systems that can be treated directly
using FCI diagonalization. The benzene molecule in the cc-pVDZ basis set
contains 108 orbitals and 30 electrons (in the frozen-core approximation).

Eriksen _et. al._ Eriksen  _et al._ (2020) performed a blind test on the
benzene molecule comparing ten different methods. They include the coupled-
cluster expansion methods, different selected CI methods, and quantum Monte
Carlo methods. The authors agreed that the full coupled-cluster reduction
(FCCR) and many-body FCI expansion (MBE-FCI) essentially yield the exact
correlation energy of $-863.0\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$.
Lee _et. al._ Lee, Malone, and Reichman (2020) supplemented the study with ph-
AFQMC results using two different trial wavefunctions: RHF wavefunction and
CASSCF(6,6) multiconfigurational wavefunction. The ph-AFQMC+RHF overestimates
the correlation energy by $3.1\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$,
while the ph-AFQMC+CAS(6,6) over-correlates by $1.2\text{\,}\mathrm{m\text{
$E_{\textup{ h }}$ }}$. For the sake of completeness, we augmented the study
with CCSD(T) results which are under-correlated by
$3.5\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$. Our ph-AFQMC+RHF result
is equivalent to the ph-AFQMC+RHF result calculated using the QMCPACK
code.Kent  _et al._ (2020); Kim  _et al._ (2018) This serves as an additional
validation of our AFQMC implementation in QMCFort. Similar to CN, CO, and CO2,
ph∗-AFQMC+RHF reduces the absolute value of the correlation energy
considerably and leads to an energy under-correlated by
$1.3\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$. In this case the quality
of the ph∗-AFQMC+RHF is seemingly similar to the quality of ph-AFQMC+CAS(6,6).
However, better trial wavefunctions promise improved accuracy and, most
importantly, more controlled results. All methods and respective correlation
energies are visualized in Fig. 7.

###  IV.3 Water Clusters

Great effort has been put into developing empirical models that faithfully
represent the properties of bulk water. It turns out that a good model must
also adequately describe the small water clusters present in the Earth’s
atmosphere. Temeslo _et. al._ collected representative structures for water
clusters of various sizes.Temelso, Archer, and Shields (2011) They
investigated which contributions are important to obtain accurate formation
energies. Motta _et al._ reported deviations of AFQMC and CCSD(T) larger than
the statistical fluctuations.Motta  _et al._ (2019) Here, we revisit these
clusters to scrutinize the accuracy of AFQMC for water.

We estimate the binding energy of the most stable water cluster with $n\leq 5$
H2O molecules

$E_{b}(n\text{H}_{2}\text{O})=E(n\text{H}_{2}\text{O})-nE(\text{H}_{2}\text{O})~{}.$
(29)

Here, $E(n\text{H}_{2}\text{O})$ is the total energy of the cluster, and
$E(\text{H}_{2}\text{O})$ is the ground-state energy of the water molecule.
The cluster geometries are taken from Ref. Temelso, Archer, and Shields, 2011.
We relaxed the single water molecule using MP2 at aug-cc-pVDZ basis set and
obtain a bond length of 0.96593 Å and a bond angle of 103.866$\degree$. We
used the heavy-augmented basis set (aug-cc-pVDZ for O atom and cc-pVDZ for H
atom) and all-electron wavefunctions to compare our results with Ref. Motta
_et al._ , 2019. The PySCF packageSun  _et al._ (2018) produced the CCSD(T)
correlation energies and QMCFort the RHF, MP2, and AFQMC ones. We truncated
the Cholesky decomposition of the ERIs with a threshold of
$10^{-6}$$\mathrm{\text{ $E_{\textup{ h }}$ }}$. We propagated 4 800 walkers
for 140 000 steps with a time-step of $0.01\text{\,}\mathrm{\text{
$E_{\textup{ h }}$ }}^{-1}$. The exchange energy was evaluated after every 100
steps. For the largest cluster, we used half as many walkers. This setup keeps
the systematic errors in the binding energy within
$0.05\text{\,}\mathrm{k}\mathrm{cal}\mathrm{/}\mathrm{mol}$.

Table 3 lists RHF, MP2, CCSD(T), ph-AFQMC, ph∗-AFQMC and AFQMC binding
energies from Ref. Motta  _et al._ , 2019. While the RHF binding energies
differ significantly from the correlation-consistent methods, all other
methods agree within chemical accuracy. Our AFQMC values are in better
agreement with MP2 and CCSD(T) values than the previous AFQMC values. The
large statistical errors in the reference data might partially explain the
discrepancy in AFQMC binding energies. However, the systematic under-
correlation of the reference AFQMC binding energies indicates the existence of
systematic errors, too. One possible source of the systematic errors could be
the looser threshold in the Cholesky decomposition of the ERIs
($10^{-4}$$\mathrm{\text{ $E_{\textup{ h }}$ }}$) in Ref. Motta  _et al._ ,
2019.

Table 3:  Binding energies of the four most stable water clusters containing up to 5 H2O molecules calculated using heavy-augmented cc-pVDZ basis set and all-electron wavefunctions. RHF, MP2, CCSD(T), and different ph-AFQMC values in $\text{kcal}/\text{mol}$ are reported. Cluster | RHF | MP2 | CCSD(T) | ph-AFQMC | ph∗-AFQMC | ph-AFQMCMotta  _et al._ (2019)
---|---|---|---|---|---|---
2Cs | $-3.82$ | $-5.22$ | $-5.18$ | $-5.17(5)$ | $-5.06(7)$ | $-5.11(31)$
3UUD | $-10.52$ | $-15.83$ | $-15.62$ | $-15.67(9)$ | $-15.68(9)$ | $-14.78(64)$
4S4 | $-19.00$ | $-28.36$ | $-27.87$ | $-28.12(10)$ | $-28.11(10)$ | $-26.49(46)$
5CYC | $-25.30$ | $-37.48$ | $-36.78$ | $-37.14(28)$ | $-37.31(28)$ | $-36.27(59)$

To test ph-AFQMC and ph∗-AFQMC for size consistency we dissociated the cluster
into individual molecules by shifting the second and third molecule in the
3UUD cluster by 8 Å in $x$ and $y$ direction, respectively. We computed the
total energy of the stretched cluster and of the individual water molecules
individually. We summarize our findings in Table 4. Within statistical
fluctuations, all fragments exhibit the same total energy which is equal to a
third of the stretched cluster. This demonstrates that ph-AFQMC and ph∗-AFQMC
are size-consistent.

Table 4:  The AFQMC energies of the stretched 3UUD cluster and its fragments, and the binding energies demonstrate that ph-AFQMC and ph∗-AFQMC are both size-consistent (all values are given in $\mathrm{\text{ $E_{\textup{ h }}$ }}$). Cluster | ph-AFQMC | ph∗-AFQMC
---|---|---
3UUD | $-228.82750(12)$ | $-228.82448(12)$
1-H2O | $-76.27576(4)$ | $-76.27475(4)$
2-H2O | $-76.27580(4)$ | $-76.27477(4)$
3-H2O | $-76.27585(4)$ | $-76.27481(4)$
$E_{b}$ | $0.00009(14)$ | $0.00015(14)$

##  V Conclusion

The phaseless auxiliary-field quantum Monte Carlo is increasingly popular due
to its high accuracy, the low polynomial scaling ($N^{3}-N^{4}$), and its
applicability to quantum chemistry and condensed matter physics. Using the DFT
or HF solutions as starting point, it can be considered as a natural extension
of these methods with similar scaling but higher accuracy.

We have presented a Fortran implementation of the AFQMC, QMCFort, that enables
efficient large-scale calculations on CPUs. The code is parallelized using
MPI, OpenMP, and parallel BLAS, and runs near peak performance for typical
systems (see Fig. 2). QMCFort can run AFQMC simulations independently or
obtain Cholesky vectors of the two-electron intermediates via interfaces to
VASP and PySCF.

Using QMCFort, we compared the accuracy of the ph-AFQMC method to the "gold-
standard" CCSD(T) and the more accurate CCSDTQP method. For this purpose, we
calculated the ph-AFQMC energies of the 26 molecules in the HEAT set, the
benzene molecule, and water clusters. For the HEAT set, the ph-AFQMC yields a
mean absolute deviation (MAD) of $1.85\text{\,}\mathrm{m\text{ $E_{\textup{ h
}}$ }}$ ($1.15\text{\,}\mathrm{k}\mathrm{cal}\mathrm{/}\mathrm{mol}$) compared
to CCSDTQP values. The poor performance for four molecules (CN, CO, CO2, O2)
is responsible for a large part of this deviation. Such outliers are clearly
troublesome as they are difficult to recognize in the absence of reference
calculations. Furthermore, CCSD(T) is certainly more robust for the HEAT set
and is highly accurate even if the groundstate wavefunction involves
significant double excitations. One possible explanation for the presence of
outliers and the failure of the AFQMC is that the method fails to account
accurately for strong double excitations. This suggests that the AFQMC with a
single determinantal trial wavefunction is only capable to treat weakly
correlated materials that lack strong static correlation effects even in the
space of double excitations.

For the remaining molecules in the HEAT set, ph-AFQMC performs similarly to or
better than the CCSD(T). We modified the phaseless approximation—ph∗-AFQMC—to
overcome at least partly the over-correlation problems often encountered with
ph-AFQMC. For the CN, CO, and CO2 molecules, where the over-correlation
effects are particularly pronounced, the modified approach indeed improves the
energies significantly. In the case of the benzene molecule, ph∗-AFQMC yields
a correlation energy $1.3\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$
higher than the FCI reference value. This is noticeably better than ph-AFQMC
and comparable to the accuracy of ph-AFQMC with a multi-determinant CAS(6,6)
trial wavefunction. Finally, the AFQMC binding energies of water clusters
agree with the CCSD(T) and MP2 binding energies to within
$0.5\text{\,}\mathrm{k}\mathrm{cal}\mathrm{/}\mathrm{mol}$. For the four water
clusters, our present results are generally closer to CCSD(T) results than
previous AFQMC results, substantiating our claim that the AFQMC method with a
single Slater determinant is particularly suitable for weakly correlated
molecules and materials.

Both the ph∗-AFQMC, as well as the ph-AFQMC, are rather ad hoc approaches to
deal with the fermionic sign problem and the exponential increase of the
noise. We feel that the present work does not conclusively show that the ph∗
approximation is superior to the ph approximation. However, the present work
clearly shows that even minor changes to the phaseless approximation can
affect the final results. Although we could not find an approximation that
more generally improves the results, we believe that further research in this
direction is well warranted.

In summary, the ph-AFQMC is a promising and reliable method with high accuracy
and reasonable computational cost. It provides nearly chemical accuracy for
the small molecules studied in this work if the molecules are weakly
correlated. We plan to further improve the accuracy of QMCFort through
algorithmic improvements and the use of more accurate trial wavefunctions.

##  VI Acknowledgments

Funding by the Austrian Science Foundation (FWF) within the project P 33440 is
gratefully acknowledged. All calculations were performed on the VSC4 / VSC5
(Vienna scientific cluster).

## Author declarations

### Conflict of Interest

The authors have no conflicts to disclose.

## Data availability

The data that support the findings of this study are available within the
article.

*

##  Appendix A Interface with VASP

Consider the set of mean-field orbitals $\\{\phi_{p}(\textbf{r})\\}$. The
single-particle Hamiltonian matrix elements are computed as

$h_{pq}=\int\text{d}^{3}(\textbf{r})\;\phi_{p}^{*}(\textbf{r})\hat{H}_{1}\phi_{q}(\textbf{r})~{},$
(30)

where $\hat{H}_{1}$ contains the kinetic energy and the Coulomb attraction
between an electron and the nuclei. To obtain objects equivalent to Cholesky
vectors $L_{g,pq}$, we introduce two-orbital densities

$\rho_{\textbf{G},pq}=\frac{1}{\sqrt{\Omega}}\int\text{d}^{3}(\textbf{r})\;e^{i\textbf{G}\textbf{r}}\;\phi_{p}^{*}(\textbf{r})\phi_{q}(\textbf{r})~{},$
(31)

where $\\{\textbf{G}\\}$ is the set of reciprocal lattice vectors defined by
the cutoff energy $\textbf{G}^{2}/2\leq E_{\text{cut}}$ and $\Omega$ is the
volume of the system. The relation between Cholesky vectors and two-orbital
densities is

$\sum_{g}^{N_{g}}L_{g,pq}L_{g,rs}=\sum_{\textbf{G}}^{N_{\textbf{G}}}\frac{4\pi}{\textbf{G}^{2}}\;\rho_{\textbf{G},pq}\rho_{\textbf{G},rs}~{},$
(32)

where $N_{g}\leq\text{min}(N_{\textbf{G}},N^{2})$. In Ref. Taheridehkordi,
2023, we introduced a truncated singular value decomposition to obtain the
Cholesky vectors

$L_{g,pq}=\text{SVD}\left(\frac{\sqrt{4\pi}}{|\textbf{G}|}\rho_{\textbf{G},pq}\right)~{}.$
(33)

## References

  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Kim  _et al._ (2015) M.-C. Kim, H. Park, S. Son, E. Sim, and K. Burke, J Phys. Chem. Lett. 6, 3802 (2015), pMID: 26722874\.
  * Cohen, Mori-Sánchez, and Yang (2008) A. J. Cohen, P. Mori-Sánchez, and W. Yang, Science 321, 792 (2008).
  * Perdew and Schmidt (2001) J. P. Perdew and K. Schmidt, AIP Conf. Proc. 577, 1 (2001), https://aip.scitation.org/doi/pdf/10.1063/1.1390175 .
  * Nagai, Akashi, and Sugino (2020) R. Nagai, R. Akashi, and O. Sugino, NPJ Comput. Mater. 6 (2020), 10.1038/s41524-020-0310-0.
  * Pederson, Kalita, and Burke (2022) R. Pederson, B. Kalita, and K. Burke, Nat. Rev. Phys. 4, 357 (2022).
  * Behler and Parrinello (2007) J. Behler and M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007).
  * Chen  _et al._ (0) M. S. Chen, J. Lee, H.-Z. Ye, T. C. Berkelbach, D. R. Reichman, and T. E. Markland, J. Chem. Theory Comput. 0, null (0), pMID: 36730728.
  * Szabo and Ostlund (1996) A. Szabo and N. S. Ostlund,  _Modern quantum chemistry: Introduction to advanced electronic structure theory_ (Dover Publications, New York, USA, 1996).
  * Helgaker, Jørgensen, and Olsen (2000) T. Helgaker, P. Jørgensen, and J. Olsen,  _Molecular Electronic Structure Theory_ (John Wiley & Sons, LTD, Chichester, 2000).
  * Handy (1980) N. Handy, Chem. Phys. Lett. 74, 280 (1980).
  * Holmes, Tubman, and Umrigar (2016) A. A. Holmes, N. M. Tubman, and C. J. Umrigar, J. Chem. Theory Comput. 12, 3674 (2016), pMID: 27428771\.
  * Sharma  _et al._ (2017) S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi, and C. J. Umrigar, J. Chem. Theory Comput. 13, 1595 (2017), pMID: 28263594\.
  * Dash  _et al._ (2018) M. Dash, S. Moroni, A. Scemama, and C. Filippi, J. Chem. Theory Comput. 14, 4176 (2018), pMID: 29953810\.
  * Booth, Thom, and Alavi (2009) G. H. Booth, A. J. Thom, and A. Alavi, J. Chem. Phys. 131, 054106 (2009).
  * Cleland, Booth, and Alavi (2010) D. Cleland, G. H. Booth, and A. Alavi, J. Chem. Phys. 132, 041103 (2010).
  * Ghanem, Guther, and Alavi (2020) K. Ghanem, K. Guther, and A. Alavi, J. Chem. Phys. 153, 224115 (2020).
  * Purvis and Bartlett (1982) G. D. Purvis and R. J. Bartlett, J. Chem. Phys. 76, 1910 (1982).
  * Bartlett and Musiał (2007) R. J. Bartlett and M. Musiał, Rev. Mod. Phys. 79, 291 (2007).
  * Tajti  _et al._ (2004) A. Tajti, P. G. Szalay, A. G. Császár, M. Kállay, J. Gauss, E. F. Valeev, B. A. Flowers, J. Vázquez, and J. F. Stanton, J. Chem. Phys. 121, 11599 (2004).
  * Bomble  _et al._ (2005) Y. J. Bomble, J. F. Stanton, M. Kállay, and J. Gauss, J. Chem. Phys. 123, 054101 (2005).
  * Harding  _et al._ (2008) M. E. Harding, J. Vázquez, B. Ruscic, A. K. Wilson, J. Gauss, and J. F. Stanton, J. Chem. Phys. 128, 114111 (2008).
  * White and Martin (1999) S. R. White and R. L. Martin, J. Chem. Phys. 110, 4127 (1999).
  * Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
  * Chan and Sharma (2011) G. K.-L. Chan and S. Sharma, Annu. Rev. Phys. Chem. 62, 465 (2011), pMID: 21219144.
  * McMillan (1965) W. L. McMillan, Physical Review 138, A442 (1965).
  * Anderson (1976) J. B. Anderson, J. Chem. Phys. 65, 4121 (1976).
  * Ceperley, Chester, and Kalos (1977) D. Ceperley, G. Chester, and M. Kalos, Phys. Rev. B 16, 3081 (1977).
  * Foulkes  _et al._ (2001) W. Foulkes, L. Mitas, R. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  * Nemec, Towler, and Needs (2010) N. Nemec, M. D. Towler, and R. J. Needs, J. Chem. Phys. 132, 034111 (2010).
  * Curtiss  _et al._ (1991) L. A. Curtiss, K. Raghavachari, G. W. Trucks, and J. A. Pople, J. Chem. Phys. 94, 7221 (1991).
  * Pfau  _et al._ (2020) D. Pfau, J. S. Spencer, A. G. D. G. Matthews, and W. M. C. Foulkes, Phys. Rev. Research 2, 033429 (2020).
  * Spencer  _et al._ (2020) J. S. Spencer, D. Pfau, A. Botev, and W. M. Foulkes, “Better, faster fermionic neural networks,”  (2020), arXiv:2011.07125 [physics.comp-ph] .
  * Gerard  _et al._ (2022) L. Gerard, M. Scherbela, P. Marquetand, and P. Grohs, “Gold-standard solutions to the schrödinger equation using deep learning: How much physics do we need?”  (2022).
  * Al-Saidi, Zhang, and Krakauer (2007) W. A. Al-Saidi, S. Zhang, and H. Krakauer, J. Chem. Phys. 127, 144101 (2007).
  * Purwanto, Zhang, and Krakauer (2009) W. Purwanto, S. Zhang, and H. Krakauer, J. Chem. Phys. 130, 094107 (2009).
  * Purwanto  _et al._ (2011) W. Purwanto, H. Krakauer, Y. Virgus, and S. Zhang, J. Chem. Phys. 135, 164105 (2011).
  * Virgus  _et al._ (2012) Y. Virgus, W. Purwanto, H. Krakauer, and S. Zhang, Phys. Rev. B 86, 241406 (2012).
  * Purwanto, Zhang, and Krakauer (2015) W. Purwanto, S. Zhang, and H. Krakauer, J. Chem. Phys. 142, 064302 (2015).
  * Motta  _et al._ (2017) M. Motta, D. M. Ceperley, G. K.-L. Chan, J. A. Gomez, E. Gull, S. Guo, C. A. Jiménez-Hoyos, T. N. Lan, J. Li, F. Ma, A. J. Millis, N. V. Prokof’ev, U. Ray, G. E. Scuseria, S. Sorella, E. M. Stoudenmire, Q. Sun, I. S. Tupitsyn, S. R. White, D. Zgid, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 7, 031059 (2017).
  * Lee, Malone, and Reichman (2020) J. Lee, F. D. Malone, and D. R. Reichman, J. Chem. Phys. 153, 126101 (2020).
  * Lee, Malone, and Morales (2020) J. Lee, F. D. Malone, and M. A. Morales, J. Chem. Theory Comput. 16, 3019 (2020), pMID: 32283932\.
  * Rudshteyn  _et al._ (2022) B. Rudshteyn, J. L. Weber, D. Coskun, P. A. Devlaminck, S. Zhang, D. R. Reichman, J. Shee, and R. A. Friesner, J. Chem. Theory Comput. 18, 2845 (2022), pMID: 35377642\.
  * Hirsch (1985) J. E. Hirsch, Phys. Rev. B 31, 4403 (1985).
  * Sugiyama and Koonin (1986) G. Sugiyama and S. Koonin, Ann. Phys. 168, 1 (1986).
  * White  _et al._ (1989) S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis, and R. T. Scalettar, Phys. Rev. B 40, 506 (1989).
  * Sorella  _et al._ (1989) S. Sorella, S. Baroni, R. Car, and M. Parrinello, Europhys. Lett. 8, 663 (1989).
  * Zhang, Carlson, and Gubernatis (1995) S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. Lett. 74, 3652 (1995).
  * Zhang, Carlson, and Gubernatis (1997) S. Zhang, J. Carlson, and J. Gubernatis, Phys. Rev. B 55, 7464 (1997).
  * Zhang and Krakauer (2003) S. Zhang and H. Krakauer, Phys. Rev. Lett. 90, 136401 (2003).
  * Suewattana  _et al._ (2007) M. Suewattana, W. Purwanto, S. Zhang, H. Krakauer, and E. J. Walter, Phys. Rev. B 75, 245123 (2007).
  * Motta  _et al._ (2019) M. Motta, J. Shee, S. Zhang, and G. K.-L. Chan, J. Chem. Theory Comput. 15, 3510 (2019), pMID: 31091103\.
  * Lee and Reichman (2020) J. Lee and D. R. Reichman, J. Chem. Phys. 153, 044131 (2020).
  * Malone, Zhang, and Morales (2019) F. D. Malone, S. Zhang, and M. A. Morales, J. Chem. Theory Comput. 15, 256 (2019).
  * Weber  _et al._ (2022) J. L. Weber, H. Vuong, P. A. Devlaminck, J. Shee, J. Lee, D. R. Reichman, and R. A. Friesner, J. Chem. Theory Comput. 18, 3447 (2022), pMID: 35507769\.
  * Zhang (2013) S. Zhang, Emergent Phenomena in Correlated Matter: Autumn School Organized by the Forschungszentrum Jülich and the German Research School for Simulation Sciences at Forschungszentrum Jülich 23-27 September 2013; Lecture Notes of the Autumn School Correlated Electrons 2013 3 (2013).
  * Motta and Zhang (2018) M. Motta and S. Zhang, WIRES Comput. Mol. Sci. 8, e1364 (2018), https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1364 .
  * Shi and Zhang (2013) H. Shi and S. Zhang, Phys. Rev. B 88, 125132 (2013).
  * Chang, Rubenstein, and Morales (2016) C.-C. Chang, B. M. Rubenstein, and M. A. Morales, Phys. Rev. B 94, 235144 (2016).
  * Mahajan and Sharma (2021) A. Mahajan and S. Sharma, J. Chem. Theory Comput. 17, 4786 (2021), pMID: 34232637.
  * Mahajan, Lee, and Sharma (2022) A. Mahajan, J. Lee, and S. Sharma, J. Chem. Phys. 156, 174111 (2022).
  * Huggins  _et al._ (2022) W. J. Huggins, B. A. O’Gorman, N. C. Rubin, D. R. Reichman, R. Babbush, and J. Lee, Nature 603, 416 (2022).
  * Landinez Borda, Gomez, and Morales (2019) E. J. Landinez Borda, J. Gomez, and M. A. Morales, J. Chem. Phys. 150, 074105 (2019).
  * Beebe and Linderberg (1977) N. H. F. Beebe and J. Linderberg, Int. J. Quantum Chem. 12, 683 (1977), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.560120408 .
  * Koch, Sánchez de Merás, and Pedersen (2003) H. Koch, A. Sánchez de Merás, and T. B. Pedersen, J. Chem. Phys. 118, 9481 (2003).
  * Hubbard (1959) J. Hubbard, Phys. Rev. Lett. 3, 77 (1959).
  * Stratonovich (1957) R. Stratonovich, in  _Sov. Phys. Dokl._ , Vol. 2 (1957) p. 416.
  * Balian and Brezin (1969) R. Balian and E. Brezin, Nuovo Cimento B 64, 37 (1969).
  * Note (1) QMCFort is available from the authors upon a reasonable request.
  * Sun  _et al._ (2018) Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, WIRES Comput. Mol. Sci. 8, e1340 (2018).
  * Trotter (1959) H. Trotter, 10, 545 (1959).
  * Suziki (1976) M. Suziki, 56, 1454 (1976).
  * Calandra Buonaura and Sorella (1998) M. Calandra Buonaura and S. Sorella, Phys. Rev. B 57, 11446 (1998).
  * Hummel, Tsatsoulis, and Grüneis (2017) F. Hummel, T. Tsatsoulis, and A. Grüneis, J. Chem. Phys. 146, 124105 (2017).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Taheridehkordi (2023) A. Taheridehkordi, “Phaseless auxiliary field quantum monte carlo with a projector augmented-wave method: application to solids,” (2023), manuscript in preparation.
  * Eriksen  _et al._ (2020) J. J. Eriksen, T. A. Anderson, J. E. Deustua, K. Ghanem, D. Hait, M. R. Hoffmann, S. Lee, D. S. Levine, I. Magoulas, J. Shen, N. M. Tubman, K. B. Whaley, E. Xu, Y. Yao, N. Zhang, A. Alavi, G. K.-L. Chan, M. Head-Gordon, W. Liu, P. Piecuch, S. Sharma, S. L. Ten-no, C. J. Umrigar, and J. Gauss, J. Phys. Chem. Lett. 11, 8922 (2020), pMID: 33022176\.
  * Kent  _et al._ (2020) P. R. C. Kent, A. Annaberdiyev, A. Benali, M. C. Bennett, E. J. Landinez Borda, P. Doak, H. Hao, K. D. Jordan, J. T. Krogel, I. Kylänpää, J. Lee, Y. Luo, F. D. Malone, C. A. Melton, L. Mitas, M. A. Morales, E. Neuscamman, F. A. Reboredo, B. Rubenstein, K. Saritas, S. Upadhyay, G. Wang, S. Zhang, and L. Zhao, J. Chem. Phys. 152, 174105 (2020).
  * Kim  _et al._ (2018) J. Kim, A. D. Baczewski, T. D. Beaudet, A. Benali, M. C. Bennett, M. A. Berrill, N. S. Blunt, E. J. L. Borda, M. Casula, D. M. Ceperley, S. Chiesa, B. K. Clark, R. C. Clay, K. T. Delaney, M. Dewing, K. P. Esler, H. Hao, O. Heinonen, P. R. C. Kent, J. T. Krogel, I. Kylänpää, Y. W. Li, M. G. Lopez, Y. Luo, F. D. Malone, R. M. Martin, A. Mathuriya, J. McMinis, C. A. Melton, L. Mitas, M. A. Morales, E. Neuscamman, W. D. Parker, S. D. P. Flores, N. A. Romero, B. M. Rubenstein, J. A. R. Shea, H. Shin, L. Shulenburger, A. F. Tillack, J. P. Townsend, N. M. Tubman, B. V. D. Goetz, J. E. Vincent, D. C. Yang, Y. Yang, S. Zhang, and L. Zhao, J. Phys. Cond. Matter 30, 195901 (2018).
  * Temelso, Archer, and Shields (2011) B. Temelso, K. A. Archer, and G. C. Shields, J. Phys. Chem. A 115, 12034 (2011), pMID: 21910428, https://doi.org/10.1021/jp2069489 .