# Selected configuration interaction wave functions in phaseless auxiliary
field quantum Monte Carlo

Ankit Mahajan<EMAIL_ADDRESS>Department of Chemistry, University of
Colorado, Boulder, CO 80302, USA  Joonho Lee  Department of Chemistry,
Columbia University, New York, New York 10027, USA  Sandeep Sharma
<EMAIL_ADDRESS>Department of Chemistry, University of Colorado, Boulder, CO
80302, USA

###### Abstract

We present efficient algorithms for using selected configuration interaction
(sCI) trial wave functions in phaseless auxiliary field quantum Monte Carlo
(ph-AFQMC). These advancements, geared towards optimizing computational
performance for longer CI expansions, allow us to use up to a million
configurations feasibly in the trial state for ph-AFQMC. In one example, we
found the cost of ph-AFQMC per sample to increase only by a factor of about 3
for a calculation with $10^{4}$ configurations compared to that with a single
one, demonstrating the tiny computational overhead due to a longer expansion.
This favorable scaling allows us to study the systematic convergence of the
phaseless bias in AFQMC calculations with an increasing number of
configurations and provides a means to gauge the accuracy of ph-AFQMC with
other trial states. We also show how the scalability issues of sCI trial
states for large system sizes could be mitigated by restricting them to a
moderately sized orbital active space and leveraging the near-cancellation of
out of active space phaseless errors.

##  1 Introduction

Quantum Monte Carlo (QMC) is a powerful tool in our arsenal to tackle the
quantum many-body problem.Kalos et al. (1974); Ceperley et al. (1977);
Ceperley and Alder (1980); Nightingale and Umrigar (1998); Foulkes et al.
(2001); Booth et al. (2013); Becca and Sorella (2017) Among various QMC
approaches, phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC)Zhang and
Krakauer (2003) has emerged as an accurate and efficient method. While
originally from the condensed matter community (therein usually referred to as
constrained-path AFQMC),Fahy and Hamann (1990); Zhang et al. (1995, 1997) ph-
AFQMC has gained popularity in chemistry in recent years.Al-Saidi et al.
(2006, 2007); Suewattana et al. (2007); Purwanto et al. (2015); Motta and
Zhang (2018); Hao et al. (2018); Motta et al. (2019); Shee et al. (2019,
2019); Williams et al. (2020); Lee and Reichman (2020); Lee et al. (2020); Shi
and Zhang (2021) The accuracy and scalability of ph-AFQMC are determined in
large part by the choice of the trial wave function. The use of a trial wave
function becomes necessary for retaining statistical efficiency (in sample
complexity) to control the fermionic phase (or sign) problemLoh Jr et al.
(1990); Troyer and Wiese (2005). The constraint imposed to control the phase
problem is called the phaseless approximationZhang and Krakauer (2003). When
the trial wave function approaches the exact ground state, the corresponding
ph-AFQMC energy tends to the ground state energy. The commonly used trial wave
function for ph-AFQMC is the broken-symmetry Hartree-Fock (HF) wave function.
It scales as $O(N^{3})$ to obtain the trial wave function and $O(N^{5})$ to
perform the ph-AFQMC calculation with the trial to obtain energy for a fixed
statistical error, where $N$ is the system size. Beyond HF, one may try to use
a single determinant trial using approximate Bruckner orbitalsLee et al.
(2020), which essentially keeps the same computational scaling for ph-AFQMC.
While single determinant trial wave functions are attractive due to their
scalability, their accuracy can be limited and questionable in many
examplesAl-Saidi et al. (2007); Williams et al. (2020); Mahajan and Sharma
(2021).

On the other hand, there has been a flurry of developments in selected
configuration interaction (sCI) methods in the last few years.Giner et al.
(2013); Evangelista (2014); Holmes et al. (2016); Tubman et al. (2016) Despite
their steep (formally exponential) scaling with system size, sCI methods have
increasingly been employed to study moderately sized systems and perform large
active space correlated calculations. These methods have the capability of
generating systematically more accurate approximations to the state of
interest by increasing the number of configurations judiciously. This
attractive feature makes sCI wave functions natural candidates for trial
states in projection QMC methods. Based on significant advances in real-space
QMC algorithms with such trial states,Filippi et al. (2016); Assaraf et al.
(2017) several papers have already reported calculations with long sCI
expansions in diffusion Monte Carlo.Dash et al. (2018); Pineda Flores and
Neuscamman (2019); Benali et al. (2020) There has also been some recent
progress for using them in ph-AFQMC,Shee et al. (2018); Williams et al.
(2020); Shi and Zhang (2021) but the commonly used algorithm based on the
Sheman-Morrison-Woodbury identity has an inherent computational bottleneck
with the local energy evaluation scaling as $O(N_{c}N^{4})$, where $N_{c}$ is
the number of determinants. This steeper scaling compared to real-space
methods is due to the Gaussian orbital representation of the electron
repulsion interaction. Furthermore, the force bias evaluation scales as
$O(N_{c}N^{3})$, which can be expensive since it is performed once every time
step. Due to these, it has been challenging to scale up the ph-AFQMC
calculations beyond 100-1000 configurations.

In this work, we propose to use one of the sCI methods, heat-bath CI
(HCI)Holmes et al. (2016); Sharma et al. (2017); Smith et al. (2017), to
generate a large trial with up to $10^{6}$ configurations. We then combine
several algorithmic advances made by two of usMahajan and Sharma (2020, 2021)
using the generalized Wick’s theorem and make further improvements to
accelerate ph-AFQMC calculations with a large trial. In particular, in our new
algorithm, the local energy evaluation scales as $O(N_{c}N+N^{4})$, and the
force bias evaluation scales as $O(N_{c}+N^{3})$. To alleviate the scaling
issues of sCI wave functions, we study the efficacy of active space trial
states in calculating energy differences and properties. The active space size
is another variable to converge the phaseless error systematically in larger
systems while keeping the calculation cost manageable. We use the new
algorithm for ph-AFQMC with HCI trials to investigate the behavior of the
phaseless bias as a function of the number of configurations in the trial for
several challenging systems.

This paper is organized as follows: we first review the ph-AFQMC algorithm
(Section 2), we then discuss how once can drastically speed up the evaluation
of force bias and local energy of sCI trials using the generalized Wick’s
theorem (Section 3), we present results that show how phaseless errors in the
ground state energy and dipole moments change as a function of the number of
determinants in hydrogen chains, transition metal oxides, and a few small
molecules (Section 4), and we finally conclude (Section 5).

##  2 Phaseless auxiliary field quantum Monte Carlo

We briefly summarize the procedure for phaseless AFQMC here and refer the
reader to reference Motta and Zhang, 2018 for more details. Consider the ab
initio Hamiltonian given by

$\hat{H}=\sum_{ij}h_{ij}\hat{a}_{i}^{\dagger}\hat{a}_{j}+\frac{1}{2}\sum_{\gamma}\left(\sum_{ij}L^{\gamma}_{ij}\hat{a}_{i}^{\dagger}\hat{a}_{j}\right)^{2},$
(1)

where $h_{ij}$ are one-electron integrals and $L^{\gamma}_{ij}$ are Cholesky
decomposed two-electron integrals in an orthonormal orbital basis. We will use
letters $N$, $M$, and $X$ to denote the number of electrons, number of
orbitals, and number of Cholesky vectors, respectively. We note that in
chemical systems, empirically, $X\sim O(M)$ with a proportionality constant
usually smaller than 10. AFQMC uses the exponential form of the projector to
converge to the ground state of the Hamiltonian as

$e^{-\tau\hat{H}}|\psi_{I}\rangle\xrightarrow{\tau\rightarrow\infty}|\Psi_{0}\rangle,$
(2)

where $\tau$ is the imaginary time, $|\Psi_{0}\rangle$ is the ground state,
and $|\psi_{I}\rangle$ is an initial state such that
$\langle\psi_{I}|\psi_{0}\rangle\neq 0$. Using the Hubbard-Stratonovic
transform, the short-time exponential projector can be written as

$e^{-\Delta\tau\hat{H}}=\int d\mathbf{x}p(\mathbf{x})\hat{B}(\mathbf{x}),$ (3)

where $\mathbf{x}$ is the vector of auxiliary fields (one scalar field per
Cholesky component), $p(\mathbf{x})$ is the standard normal Gaussian
distribution of the auxiliary fields, and $\hat{B}(\mathbf{x})$ is a complex
propagator given by the exponential of a one-body operator. Due to Thouless’
theorem, $\hat{B}(\mathbf{x})$ acts on a Slater determinant $|\phi\rangle$ as

$\hat{B}(\mathbf{x})|\phi\rangle=|\phi(\mathbf{x})\rangle,$ (4)

where $|\phi(\mathbf{x})\rangle$ is another Slater determinant obtained by a
linear transformation of the orbitals in $|\phi\rangle$. By importance
sampling the auxiliary fields in equation 3 and applying the short time
propagator sufficiently many times, we get a stochastic representation of the
ground state wave function in the long time limit as

$|\Psi_{0}\rangle\propto\sum_{i}w_{i}\frac{|\phi_{i}\rangle}{\langle\psi_{T}|\phi_{i}\rangle},$
(5)

where $w_{i}$ are weights, $|\phi_{i}\rangle$ are Slater determinants with
complex orbitals obtained from the orbital transformations sampled during
propagation, and $|\psi_{T}\rangle$ is the trial state used for importance
sampling. More accurate trial states lead to less noisy weights. We use the
hybrid approximation in this paper, which avoids the expensive calculation of
local energy at each propagation step. Importance sampling in AFQMC involves
shifting the sampled auxiliary fields by the force bias given as

$\bar{x}_{\gamma}=-\sqrt{\Delta\tau}\frac{\langle\psi_{T}|\sum_{ij}L^{\gamma}_{ij}\hat{a}_{i}^{\dagger}\hat{a}_{j}|\phi\rangle}{\langle\psi_{T}|\phi\rangle},$
(6)

where $|\phi\rangle$ is the walker. The force bias can be thought of as a
dynamic correction to the mean-field contour shiftRom et al. (1997) that
vanishes as $\Delta\tau\rightarrow 0$. Force bias is not enough, by itself, to
control the large fluctuations stemming from the phase problem. The phaseless
constraint can be used to overcome the phase problem at the expense of a
systematic bias in the sampled wave function. The size of this phaseless bias
is dictated by the accuracy of the trial state.

After an equilibration time, the energy of the sampled wave function can be
measured as

$E\approx\frac{\sum_{i}w_{i}E_{L}(\phi_{i})}{\sum_{i}w_{i}},$ (7)

where

$E_{L}(\phi_{i})=\frac{\langle\psi_{T}|\hat{H}|\phi_{i}\rangle}{\langle\psi_{T}|\phi_{i}\rangle}$
(8)

is the local energy of the walker $|\phi_{i}\rangle$. We note that this
estimator has a zero-variance property i.e. in the limit that the trial state
is the exact, ground state the variance of the estimator vanishes. In
practice, this means that more accurate trial states lead to less noisy local
energy values.

It is evident that more accurate trial states lead to less noisy simulations
and smaller phaseless biases. But, in general, the use of more sophisticated
trial states comes at the expense of greater computational cost. In the
following section, we present efficient algorithms to calculate force bias and
local energy of sCI trial states.

##  3 Selected configuration interaction trial states

A selected CI wave function is given by

$|\psi_{T}\rangle=\sum_{n}^{N_{c}}c_{n}|\psi_{n}\rangle=\sum_{n}^{N_{c}}c_{n}\prod_{\mu}^{k_{n}}\hat{a}^{\dagger}_{t_{n_{\mu}}}\hat{a}_{p_{n_{\mu}}}|\psi_{0}\rangle,$
(9)

where $|\psi_{n}\rangle$ are configurations obtained by particle-hole
excitations from the reference configuration $|\psi_{0}\rangle$, $N_{c}$ is
the number of configurations in the expansion, $c_{n}$ are real expansion
coefficients, and $k_{n}$ are the excitation ranks. We will use the indices
$p_{\mu}$ and $t_{\mu}$ to denote occupied and virtual orbitals, respectively,
whereas indices $i,j,\dots$ will be used for general orbitals. The orthonormal
orbital basis set used in the expansion can be chosen to be natural orbitals
obtained from an HCI calculation or can be optimized using a self-consistent
procedure. In the following, we use the same basis set to express the
Hamiltonian in equation 1 and, in general, all second-quantized operators
refer to this basis. In our prior work, we have proposed efficient algorithms
for using selected CI states in variational Monte CarloMahajan and Sharma
(2020) and free projection AFQMC.Mahajan and Sharma (2021) These focused on
the calculation of local energy using the generalized Wick’s theorem. Here, we
discuss an algorithm for calculating the force bias required in ph-AFQMC and
briefly summarize the calculation of local energy.

###  3.1 Force bias

Calculation of the force bias for importance sampling is one of the
computationally intensive parts of propagation in ph-AFMQC. Substituting the
sCI trial state into the force bias expression in equation 6, we get

$\bar{x}_{\gamma}=-\sqrt{\Delta\tau}\frac{\sum_{n}\sum_{ij}c_{n}L^{\gamma}_{ij}\langle\psi_{n}|\hat{a}_{i}^{\dagger}\hat{a}_{j}|\phi\rangle}{\sum_{n}c_{n}\langle\psi_{n}|\phi\rangle}.$
(10)

The conventional way of evaluating the force bias proceeds by separately
calculating the contribution of each configuration in the CI expansion using
the Sherman-Morrison-Woodbury identity, resulting in a cost scaling of
$O(N_{c}NM+XNM)$.Shee et al. (2018) Here we present an algorithm with cost
scaling as $O(N_{c}+XM^{2})$ by essentially separating the sum over CI
excitations from that over Cholesky integral indices in equation 10.

We define the Green’s function matrix for the reference configuration
$|\psi_{0}\rangle$ and the walker configuration $|\phi\rangle$ asLöwdin
(1955); Balian and Brezin (1969)

$G^{i}_{j}=\frac{\langle\psi_{0}|\hat{a}_{i}^{\dagger}\hat{a}_{j}|\phi\rangle}{\langle\psi_{0}|\phi\rangle}=\left[\phi(\psi_{0}^{\dagger}\phi)^{-1}\psi_{0}^{\dagger}\right]^{j}_{i},$
(11)

where $\psi_{0}$ and $\phi$ are the orbital coefficient matrices of the
corresponding configurations, and superscripts and subscripts denote row and
column indices, respectively. Note that since we work in the orbital basis of
the CI expansion, $\psi_{0}$ has a particularly simple form with all its
columns being unit vectors. Thus the cost of calculating the Green’s function
in this basis scales as $O(N^{2}M)$. For convenience, we also define the
related quantity

$\mathcal{G}^{i}_{j}=-\frac{\langle\psi_{0}|\hat{a}_{j}\hat{a}_{i}^{\dagger}|\phi\rangle}{\langle\psi_{0}|\phi\rangle}=G^{i}_{j}-\delta^{i}_{j}.$
(12)

According to the generalized Wick’s theoremBalian and Brezin (1969), we have

$\frac{\langle\psi_{0}|\prod_{\mu}^{k}\hat{a}^{\dagger}_{p_{\mu}}\hat{a}_{t_{\mu}}|\phi\rangle}{\langle\psi_{0}|\phi\rangle}=\det\left(G^{\left\\{p_{\mu}\right\\}}_{\left\\{t_{\mu}\right\\}}\right),$
(13)

where the sets of indices $\left\\{p_{\mu}\right\\}$ and
$\left\\{t_{\mu}\right\\}$ denote the $k\times k$ slice of the $G$ matrix, $k$
being the rank of the excitation. This expression is obtained by taking
pairwise contractions of the operators in the string of excitations according
to the generalized Wick’s theorem, with the determinant structure arising due
to fermionic permutation parity factors. Therefore, the denominator in
equation 10 can be expressed as

$\frac{\langle\psi_{T}|\phi\rangle}{\langle\psi_{0}|\phi\rangle}=\sum_{n}c_{n}\det\left(G^{\left\\{p_{n_{\mu}}\right\\}}_{\left\\{t_{n_{\mu}}\right\\}}\right).$
(14)

The computational cost scaling of this overlap ratio is thus $O(N_{c})$ once
the Green’s function has been calculated. Since the Wick’s theorem applies
naturally to overlap ratios, we will find it convenient to express matrix
elements as ratios with the reference overlap $\langle\psi_{0}|\phi\rangle$.

To evaluate the numerator of equation 10, consider the matrix element ratio
given by

$\frac{\langle\psi_{n}|\hat{a}_{i}^{\dagger}\hat{a}_{j}|\phi\rangle}{\langle\psi_{0}|\phi\rangle}=\frac{\langle\psi_{0}|\left(\prod_{\mu}^{k}\hat{a}^{\dagger}_{p_{\mu}}\hat{a}_{t_{\mu}}\right)\hat{a}_{i}^{\dagger}\hat{a}_{j}|\phi\rangle}{\langle\psi_{0}|\phi\rangle},$
(15)

where we have dropped the configuration index subscript ($n$) on the
excitations for clarity. This ratio can also be evaluated using the
generalized Wick’s theorem and by choosing the appropriate order of operations
one can achieve a significant reduction in the cost of calculating the
numerator. To this end, we note that the pairwise contractions in the string
of excitation operators can be divided into two groups: $i$) those without any
contractions between CI excitations and $\hat{a}_{i}^{\dagger}\hat{a}_{j}$,
$ii$) those containing contractions between two of the CI operators and
$\hat{a}_{i}^{\dagger}\hat{a}_{j}$. Algebraically, the two terms are given as

$\begin{split}&\frac{\langle\psi_{0}|\left(\prod_{\mu}^{k}a^{\dagger}_{p_{\mu}}a_{t_{\mu}}\right)\hat{a}_{i}^{\dagger}\hat{a}_{j}|\phi\rangle}{\langle\psi_{0}|\phi\rangle}=\det\begin{pmatrix}G^{i}_{j}&\mathcal{G}^{i}_{\left\\{t_{\mu}\right\\}}\\\\[10.00002pt]
G^{\left\\{p_{\mu}\right\\}}_{j}&G^{\left\\{p_{\mu}\right\\}}_{\left\\{t_{\mu}\right\\}}\\\
\end{pmatrix}\\\ &\
=G^{i}_{j}\det\left(G^{\left\\{p_{\mu}\right\\}}_{\left\\{t_{\mu}\right\\}}\right)\\\
&\qquad\qquad\
+\sum_{\nu,\lambda}^{k}(-1)^{1+\nu+\lambda}G^{p_{\nu}}_{j}\mathcal{G}^{i}_{t_{\lambda}}\det\left(G^{\left\\{p_{\mu}\right\\}\char
92\relax p_{\nu}}_{\left\\{t_{\mu}\right\\}\char 92\relax
t_{\lambda}}\right),\end{split}$ (16)

where on the first line the matrix of size $(k+1)\times(k+1)$ is written in a
block form, and the notation $\left\\{p_{\mu}\right\\}\char 92\relax p_{\nu}$
indicates the set of indices $\left\\{p_{\mu}\right\\}$ excluding $p_{\nu}$.
Using this grouping of terms, we get for the numerator in equation 10

$\begin{split}\sum_{n}&\sum_{ij}c_{n}L_{ij}^{\gamma}\frac{\langle\psi_{n}|\hat{a}_{i}^{\dagger}\hat{a}_{j}|\phi\rangle}{\langle\psi_{0}|\phi\rangle}\\\
&=\left(\sum_{n}c_{n}\det\left(G^{\left\\{p_{\mu_{n}}\right\\}}_{\left\\{t_{\mu_{n}}\right\\}}\right)\right)\left(\sum_{ij}L^{\gamma}_{ij}G^{i}_{j}\right)\\\
&\qquad\qquad\qquad\qquad+\sum_{pt}I^{t}_{p}\left(\sum_{ij}L^{\gamma}_{ij}G^{p}_{j}\mathcal{G}^{i}_{t}\right),\end{split}$
(17)

where

$I^{t}_{p}=\sum_{n}c_{n}\sum_{\nu,\lambda}^{k_{n}}\delta_{p_{n_{\nu}},p}\delta_{t_{n_{\lambda}},t}(-1)^{1+\nu+\lambda}\det\left(G^{\left\\{p_{n_{\mu}}\right\\}\char
92\relax p_{n_{\nu}}}_{\left\\{t_{n_{\mu}}\right\\}\char 92\relax
t_{n_{\lambda}}}\right)$ (18)

is an intermediate formed by summing over the CI configurations. The two terms
in equation 17 are shown as tensor contraction diagrams in table 1 for the
case of doubly excited CI configurations. Note that a selected CI state does
not necessarily include all excitations of a given rank, but since it is not
well suited for a tensor representation we use all doubly excited
configurations in the table. In practice, the sum over CI excitations in table
1 is only performed over the configurations in the selected CI state. The
first term has a cost scaling of $O(N_{c}+XNM)$, while the second one using
the intermediate has a scaling of $O(N_{c}+XM^{2})$. If the selected CI state
is restricted to an active space of size $A$, the overall cost can be reduced
to $O(N_{c}+XNM+XAM)$.

Table 1: Summary of computational cost scaling of force bias and two-body local energy calculations (given the reference Green’s function matrix $G$). SMW refers to overall cost scaling of the algorithm that makes use of the Sherman-Morrison-Woodbury formula as described in reference Shee et al., 2018. The last columns shows tensor contraction diagrams for all distinct types of pairwise Wick’s contractions using all doubly excited configurations $\sum c_{ptqu}\hat{a}_{t}^{\dagger}\hat{a}_{p}\hat{a}_{u}^{\dagger}\hat{a}_{q}|\psi_{0}\rangle$ as an example. Combining all possible contractions of a given type along with the fermionic parity signs (not shown here) leads to the determinant expressions in equations 17 and 20. The same types of terms arise for higher than doubly excited configurations and the adjacent column shows the optimal cost scaling for calculating each type for a general CI trial with $N_{c}$ configurations.  Quantity |  | SMW |  | This work
---|---|---|---|---
|  |  |  | Scaling |  | Tensor diagrams for doubly excited configurations
Force bias |  | $NM(N_{c}+X)$ |  | $N_{c}+XNM$ |  |
|  | $N_{c}+XM^{2}$ |  |
Local energy |  | $N_{c}XN^{2}M$ |  | $N_{c}+XNM$ |  |
|  | $N_{c}+XN^{2}M$ |  |
|  | $N_{c}+XM^{2}$ |  |
|  | $N_{c}+XNM^{2}$ |  |
|  | $X(N_{c}+NM^{2})$ |  |
Key for notation
$N_{c}$: number of configurations, $N$: number of electrons, $M$: number of
orbitals, $X$: number of Cholesky vectors
$c$: CI coefficients, $L$: Cholesky vectors, $G$ and $\mathcal{G}$: Green’s
functions
$i,j,k,l$: Hamiltonian indices; $p,q,t,u$: CI excitation indices

###  3.2 Local energy

The calculation of local energy was discussed in detail in reference Mahajan
and Sharma, 2021. Here, we give a brief summary of the algorithm. Consider the
two-body part of the local energy (see equation 8), expressed as

$E^{2}_{L}\left[\phi\right]=\frac{\sum_{n}\sum_{\gamma
ijkl}c_{n}L^{\gamma}_{ij}L^{\gamma}_{kl}\frac{\langle\psi_{n}|\hat{a}_{i}^{\dagger}\hat{a}_{k}^{\dagger}\hat{a}_{l}\hat{a}_{j}|\phi\rangle}{\langle\psi_{0}|\phi\rangle}}{\sum_{n}c_{n}\frac{\langle\psi_{n}|\phi\rangle}{\langle\psi_{0}|\phi\rangle}}.$
(19)

The calculation of the overlap ratios in the denominator is described in the
last section for the force bias. The matrix element in the numerator for the
$n$th configuration is given by using the Wick’s theorem as

$\frac{\langle\psi_{0}|\left(\prod_{\mu}^{k}\hat{a}^{\dagger}_{p_{\mu}}\hat{a}_{t_{\mu}}\right)\hat{a}_{i}^{\dagger}\hat{a}_{k}^{\dagger}\hat{a}_{l}\hat{a}_{j}|\phi\rangle}{\langle\psi_{0}|\phi\rangle}=\det\begin{pmatrix}G^{\left\\{i,k\right\\}}_{\left\\{j,l\right\\}}&\mathcal{G}^{\left\\{i,k\right\\}}_{\left\\{t_{\mu}\right\\}}\\\\[10.00002pt]
G^{\left\\{p_{\mu}\right\\}}_{\left\\{j,l\right\\}}&G^{\left\\{p_{\mu}\right\\}}_{\left\\{t_{\mu}\right\\}}\\\
\end{pmatrix},$ (20)

where we have again dropped the configuration index subscript. Similar to the
force bias calculation, the terms in this determinant can be split into two
groups based on the kinds of pairwise contractions: those that do not involve
cross contractions between CI and Hamiltonian excitation operators and those
that involve at least one such cross contraction. The second group containing
cross contractions can be further split into two groups depending on whether
one or both of the Hamiltonian excitations are cross contracted with CI
excitations. The resulting terms are shown as a tensor diagram in table 1.
Again, we have used all double excited CI configurations only for
representational convenience. Different orders of performing the tensor
contractions lead to different cost scalings, and using the one described in
reference Mahajan and Sharma, 2021 leads to a scaling of $O(N_{c}X+XNM^{2})$.
For a selected CI state restricted to an active space of size $A$, this cost
is reduced to $O(N_{c}X+XNAM+XN^{2}M)$.

##  4 Results

In this section, we present the results of our ph-AFQMC/HCI calculations and
analyze the convergence of the phaseless error for several systems. In
particular, we consider the utility of active space trial states for obtaining
accurate energy differences and dipole moments. We used PySCFSun et al. (2018)
to obtain molecular integrals and to perform all quantum chemistry wave
function calculations. The SHCI code Dice was used to obtain the trial states.
The code used to perform ph-AFQMC calculations is available in a public
repository.dqm  Input and output files for all calculations can also be
accessed from a public repository.afq  Sources of systematic errors in ph-
AFQMC calculations, besides the phaseless bias, include Trotter error,
population control bias, errors in the Cholesky decomposition of the electron
repulsion integrals, matrix exponential Taylor series truncation errors, and
bias due to filtering of rare large fluctuations. We used a conservative time
step of 0.005 a.u in all ph-AFQMC calculations. For population control, we
used the reconfiguration procedure described in reference Buonaura and
Sorella, 1998. Cholesky decompositions were calculated up to a threshold error
of $10^{-5}$. Matrix exponentials were calculated by keeping terms up to the
sixth power in the Taylor expansion. We filtered large fluctuations by capping
weights and local energies.Motta and Zhang (2018) Although it is difficult to
get an accurate estimate of the systematic errors due to all these factors in
general, they can be controlled systematically. We estimate these errors to be
much smaller than the statistical errors in all the results presented here.

###  4.1 Hydrogen chains

Figure 1: Convergence of the ph-AFQMC/HCI phaseless energy error and scaling
of the cost of computation with increasing number of configurations for H50
using the sto-6g basis at $d=1.6$ Bohr. The phaseless error is calculated with
respect to the near-exact DMRG energy, and the cost factor is calculated as
the ratio of computational costs with respect to a single configuration
calculation.

Hydrogen chain systems are convenient for benchmarking due to their simplicity
and availability of accurate density matrix renormalization group (DMRG)
ground state energies.Motta et al. (2017) First, we consider a chain of 50
equidistant Hydrogen atoms with an interatomic distance $d=1.6$ Bohr in the
minimal sto-6g basis. We performed ph-AFQMC/HCI calculations on this system
using progressively larger HCI expansion. These trial states were constructed
as truncations of the HCI state obtained with $\epsilon_{1}=10^{-4}$ in the
canonical RHF basis. Figure 1 shows the convergence of the phaseless error
(calculated as a difference to the DMRG energy reported in reference ) with
the number of configurations in the trial state. We see a nearly monotonic
decrease in the phaseless error and the error is less than 2 mH for the trial
with $10^{6}$ configurations. The figure also shows the scaling of the cost of
the ph-AFQMC calculation against the number of configurations with respect to
the cost of a single configuration calculation. Remarkably, the cost factor
for up to $10^{4}$ configurations is less than 10 and increases roughly
linearly with the number of configurations thereafter. This observation is
consistent with the scaling relations discussed in section 3, with the
$O(N_{c}X)$ factor in the cost scaling of local energy evaluation dominating
for large HCI expansions. We found changing the value of $\epsilon_{1}$ or
slightly changing the selection criterion for choosing the configurations to
only make a difference in ph-AFQMC energies for smaller number of
configurations, leading to similar convergence behaviors for longer
expansions.

Figure 2: ph-AFQMC/HCI phaseless errors for H10 using the cc-pVDZ basis set at
different bond lengths. The trial states with different number of
configurations are constructed by truncating a (10e, 10o) CASSCF wave function
in each case. The hollow symbols show ph-AFQMC/UHF phaseless errors. (They are
almost zero for $d=1.6$ Bohr and $d=2.4$ Bohr.)

We also studied the H10 chain in the cc-pVDZ basis at different bond lengths
to gauge the performance of active space trial states. We first performed a
(10e, 10o) CASSCF (complete active space self consistent field) calculation
and generated trial states by truncating the active space wave function.
Figure 2 shows the convergence of the phaseless error with the size of the HCI
trial. For all three bond lengths, we see a systematic convergence to an
almost vanishing error with the active space trial. The shortest bond length
near equilibrium converges rapidly, while the stretched bond length requires
almost the full active space wave function to reach convergence. The out of
active space phaseless error, arising due to not correlating any orbitals
outside the active space in the trial state, is tiny in this case and, more
importantly, is very similar for all three bond lengths. We also show ph-
AFQMC/UHF phaseless errors in the figure for comparison. ph-AFQMC/UHF is known
to be very accurate for hydrogen chains except at stretched geometries. At
stretched geometries, a CI expansion based on an unrestricted HF reference has
been employed in the past as a trial stateMotta et al. (2017) and could be a
handy tool in other problems as well. The algorithms described in this work
can be straightforwardly extended to work with UHF-based CI expansions.

###  4.2 Transition metal oxides

A recent study of transition metal oxide molecules reported benchmark ground
state energies using various accurate many-body methods.Williams et al. (2020)
Here, we use ph-AFQMC/HCI to study convergence of the phaseless errors for
active space trials in TiO, CrO, MnO, and FeO. The Trails-Needs
pseudopotential and the corresponding dz basis setsTrail and Needs (2017) were
used at equilibrium geometries used in reference Williams et al., 2020.

Figure 3 shows the convergence of phaseless errors with the number of
configurations in the trial state at equilibrium geometries used in reference
Williams et al., 2020. We first performed CASSCF calculations with active
spaces consisting of the metal $3d$ and $4d$ orbitals and oxygen $2p$
orbitals. The HCI trial states were then generated in a space including the
CASSCF active orbitals and the core orbitals, including all electrons.
$\epsilon_{1}=10^{-5}$ was used for the HCI calculations. Interestingly, for
FeO, CrO, and MnO the phaseless error first increases before converging to a
fixed value as the number of configurations in the trial is increased. This
peculiar behavior highlights the fact that a lower energy trial state does not
guarantee a smaller phaseless error. We note that the trends for smaller
expansions are dependent on the particular schemes used for choosing the
orbital spaces and for truncating the HCI states. But results for a large
number of configurations are largely independent of such somewhat subjective
choices. For TiO, using more configurations makes little difference to the
phaseless error. We also performed ph-AFQMC/UHF calculations for comparison
and found the phaseless errors to be 2(1), 5.3(8), -10(1), and 30.6(9) mH for
TiO, CrO, MnO, and FeO, respectively.

Figure 3: Convergence of ph-AFQMC/HCI phaseless energy errors for four
transition metal oxide molecules using the Trail-Needs dz pseudopotential and
basis set. SHCI reference energies were used as references for calculating the
phaseless errors. The trial states were constructed by truncating active space
HCI wave functions. Hollow symbols show ph-AFQMC/UHF errors.

The residual phaseless errors can be attributed solely to the out of HCI
active space correlation. These errors can be systematically reduced by
increasing the size of the active space, but getting vanishing errors this way
is unlikely to be practical for larger systems. On the other hand, as noted in
section 4.1, it is reasonable to expect near-cancellation of this out of
active space error when energy differences at different geometries are
calculated using the same active space. Similar trends have been observed in
multireference perturbation theories. To gauge the efficacy of this heuristic,
we calculated ph-AFQMC/HCI energies of CrO at three bond lengths using the
procedure described above. Table 2 shows that the converged phaseless errors
are almost identical at the three points. We note that it is crucial to
converge the phaseless error of the active space trial, and different
geometries may require trial states with different numbers of configurations
to achieve convergence. This evidence suggests that in large problems,
scalability issues of HCI trial states can be mitigated using active spaces.

Table 2: CrO ground state energies and out of active space phaseless errors at different bond lengths. Bond lengths ($d$) are in Bohr, ph-AFQMC/HCI and SHCI energies are in H, the phaseless errors are in mH. $d$ |  | ph-AFQMC/HCI |  | SHCI |  | Phaseless error
---|---|---|---|---|---|---
2.65 |  | -102.4993(5) |  | -102.5040(5) |  | 4.6(7)
3.06 |  | -102.5531(5) |  | -102.5584Williams et al. (2020) |  | 5.3(5)
3.59 |  | -102.5269(7) |  | -102.5325(5) |  | 5.6(8)

###  4.3 Dipole moments

The mixed estimator of an observable $\hat{O}$ in the ground state is given by

$\langle
O\rangle_{\text{mixed}}=\frac{\langle\psi_{T}|\hat{O}|\Psi_{0}\rangle}{\langle\psi_{T}|\Psi_{0}\rangle},$
(21)

where $|\Psi_{0}\rangle$ is the ground state. In ph-AFQMC, one has access to
an approximate sampling of the ground state, making the calculation of this
mixed estimator convenient:

$\langle
O\rangle_{\text{mixed}}\approx\frac{\sum_{i}w_{i}O_{L}(\phi_{i})}{\sum_{i}w_{i}},$
(22)

where

$O_{L}(\phi_{i})=\frac{\langle\psi_{T}|\hat{O}|\phi_{i}\rangle}{\langle\psi_{T}|\phi_{i}\rangle}$
(23)

is the local observable value.

For observables $\hat{O}$ that do not commute with the system Hamiltonian, the
mixed estimator has a systematic bias due to the trial state used to measure
the observable. This bias can be reduced by using the extrapolated estimator,
often employed in diffusion Monte Carlo,Whitlock et al. (1979); Ceperley and
Kalos (1986); Rothstein (2013) given by

$\langle\hat{O}\rangle_{\text{extrapolated}}=2\frac{\langle\psi_{T}|\hat{O}|\Psi_{0}\rangle}{\langle\psi_{T}|\Psi_{0}\rangle}-\frac{\langle\psi_{T}|\hat{O}|\psi_{T}\rangle}{\langle\psi_{T}|\psi_{T}\rangle},$
(24)

where the variational estimator of the trial state is used to correct the
mixed estimator. The accuracy of this extrapolation depends critically on the
quality of the trial state. Since ph-AFQMC calculations on molecular systems
are typically performed using crude single determinant trials, extrapolation
is not a viable option, and backpropagation is used to evaluate variational
estimators of the sampled state instead.Motta and Zhang (2017) Here, we study
the behavior of the mixed and extrapolated estimators of dipole moments of a
few small molecules. Evaluation of the local dipole moment for an HCI trial
state is similar to the force bias calculation outlined in section 3.1 since
they both involve matrix elements of one-electron operators.

Figure 4: Convergence of different estimators for the dipole moment of NH3
using the 6-31g basis. Variational and mixed estimators are calculated using
HCI and ph-AFQMC/HCI, respectively.

We first look at the exactly solvable problem of NH3 in the 6-31g basis
(geometry in the SI). We fully diagonalize this (10e, 14o) problem and
calculate its ground state dipole moment. Mixed and extrapolated ph-AFQMC
estimators were evaluated with trial states consisting of an increasing number
of leading configurations from the ground state in the RHF canonical orbital
basis. Figure 4 shows the convergence of the estimators with the trial state.
In this case, we see a systematic convergence of the variational and mixed
estimators to the exact value. The extrapolation works reasonably well, with
consistently smaller errors than the mixed estimator.

Table 3: Dipole moments (in a.u.) of three molecules using the aug-cc-pVQZ basis set. The reported ph-AFQMC/HCI values are converged with the number of configurations in the 50 orbital active space trial state. Molecule |  | ph-AFQMC/HCI |  | RHF |  | MP2 |  | CCSD |  | CCSD(T) |  | Experiment
---|---|---|---|---|---|---|---|---|---|---|---|---
|  | Mixed | Extrapolated |  |  |  |  |  |  |  |  |  |
CO |  | 0.070(3) | 0.044(3) |  | -0.104 |  | 0.108 |  | 0.024 |  | 0.048 |  | 0.048Muenter (1975)
BF |  | 0.356(4) | 0.321(4) |  | 0.333 |  | 0.377 |  | 0.314 |  | 0.325 |  | -
H2O |  | 0.726(2) | 0.734(2) |  | 0.780 |  | 0.733 |  | 0.738 |  | 0.729 |  | 0.730Lide (2004)

Table 3 shows the calculated dipole moments for three molecules at equilibrium
geometries in the aug-cc-pVQZ basis. The large basis set precludes exact
evaluation, but CCSD(T) moments are known to be very accurate for the
molecules considered here and agree well with the experimental values. We
computed orbital relaxed coupled cluster and MP2 dipole moments by calculating
energies at two different electric fields and evaluating the numerical
derivative. For obtaining ph-AFQMC/HCI estimators, we first performed a full-
valence CASSCF calculations. The trial states were then constructed using HCI
to correlate all electrons in a space of 50 CASSCF orbitals. In all cases,
both mixed and extrapolated estimators were nearly converged with the number
of configurations in the active space (details can be found in the SI).
Consider the case of the challenging CO molecule, which has a small dipole
moment, and RHF predicts the wrong direction polarity in this case. The
convergence of the dipole moment estimators with number of configurations in
the trial state is shown in figure 5. For a small number of configurations
(fewer than 100), the variational estimator is qualitatively incorrect similar
to RHF. Adding more configurations flips the sign, and the variational
estimator seems to eventually converge to the asymptotic value for the
truncated space of 50 orbitals. The mixed estimator also exhibits a non-
monotonic convergence with its value for a single determinant close to the
experimental value because of a fortuitous cancellation of mixed estimator and
phaseless biases. The extrapolated estimator for small number of
configurations in the trial is very poor due to the erroneous variational
estimates. For more than $10^{3}$ configurations, while the mixed estimator
seems to be converging to a substantially biased dipole moment value, the
extrapolation rectifies this bias reasonably well. The extrapolated estimators
are in good agreement with the CCSD(T) and experimental values for all three
molecules. Evidently, extrapolation effectively corrects the out of active
space bias in the mixed estimator for CO and BF.

Figure 5: Convergence of different estimators for the dipole moment of CO
using the aug-cc-pVQZ basis. Variational and mixed estimators are calculated
using HCI and ph-AFQMC/HCI, respectively. Trial states are obtained by
truncating a (14e, 50o) active space HCI wave function.

##  5 Conclusion

In this work, we presented efficient algorithms for using selected CI trial
states in ph-AFQMC. We demonstrated their favorable scaling by showing that
simulations with long sCI expansions incur little overhead. In our analysis of
the convergence of phaseless error in ground-state energies and dipole moments
as a function of the number of configurations, we found it to be non-monotonic
in some cases, highlighting the importance of using a sufficient number of
determinants in these systems. The use of states restricted to moderately
sized active spaces yielded excellent results for energy differences and
dipole moments. Our numerical experiments suggest that this may be a practical
way to tackle larger systems, where the generation and handling of sCI states
are challenging. Finally, we showed that the extrapolated estimator for dipole
moments is accurate if enough configurations are present in the trial state.
If this behavior turns out to be true in general, extrapolated estimators
could be a viable alternative to backpropagation for calculating ground state
properties.

Our technique can be employed in different ways for studying challenging
systems. It can be used to validate the use of simpler and less expensive
trial states in ph-AFQMC for large calculations. ph-AFQMC/HCI can be adopted
as an accurate solver in embedding schemes for describing the embedded
cluster. It would be interesting to study the feasibility of the extrapolated
estimator in larger systems and for calculating two-body properties. Lastly,
we note that this technique can be utilized in the quantum-classical hybrid
AFQMC (QC-AFQMC) in prototypical and practical quantum computations.Huggins et
al. (2021)

## Data availability

The data that support the findings of this study are available in a public
repository at afq, .

## Acknowledgements

AM and SS were supported by the National Science Foundation through the grant
CHE-1800584. SS was also partly supported through the Sloan research
fellowship. JL thanks David Reichman for support and encouragement. All
calculations were performed on the Blanca and Summit clusters at CU Boulder.

## References

  * Kalos et al. (1974) Kalos, M. H.; Levesque, D.; Verlet, L. Helium at zero temperature with hard-sphere and other forces. _Physical Review A_ 1974, _9_ , 2178.
  * Ceperley et al. (1977) Ceperley, D.; Chester, G. V.; Kalos, M. H. Monte Carlo simulation of a many-fermion study. _Physical Review B_ 1977, _16_ , 3081\.
  * Ceperley and Alder (1980) Ceperley, D. M.; Alder, B. J. Ground state of the electron gas by a stochastic method. _Physical review letters_ 1980, _45_ , 566.
  * Nightingale and Umrigar (1998) Nightingale, M. P.; Umrigar, C. J. _Quantum Monte Carlo methods in physics and chemistry_ ; Springer Science & Business Media, 1998.
  * Foulkes et al. (2001) Foulkes, W.; Mitas, L.; Needs, R.; Rajagopal, G. Quantum Monte Carlo simulations of solids. _Reviews of Modern Physics_ 2001, _73_ , 33.
  * Booth et al. (2013) Booth, G. H.; Grüneis, A.; Kresse, G.; Alavi, A. Towards an exact description of electronic wavefunctions in real solids. _Nature_ 2013, _493_ , 365.
  * Becca and Sorella (2017) Becca, F.; Sorella, S. _Quantum Monte Carlo approaches for correlated systems_ ; Cambridge University Press, 2017.
  * Zhang and Krakauer (2003) Zhang, S.; Krakauer, H. Quantum Monte Carlo method using phase-free random walks with Slater determinants. _Physical review letters_ 2003, _90_ , 136401.
  * Fahy and Hamann (1990) Fahy, S.; Hamann, D. Positive-projection Monte Carlo simulation: A new variational approach to strongly interacting fermion systems. _Physical review letters_ 1990, _65_ , 3437.
  * Zhang et al. (1995) Zhang, S.; Carlson, J.; Gubernatis, J. E. Constrained path quantum Monte Carlo method for fermion ground states. _Physical review letters_ 1995, _74_ , 3652.
  * Zhang et al. (1997) Zhang, S.; Carlson, J.; Gubernatis, J. E. Constrained path Monte Carlo method for fermion ground states. _Physical Review B_ 1997, _55_ , 7464\.
  * Al-Saidi et al. (2006) Al-Saidi, W.; Zhang, S.; Krakauer, H. Auxiliary-field quantum Monte Carlo calculations of molecular systems with a Gaussian basis. _The Journal of chemical physics_ 2006, _124_ , 224101.
  * Al-Saidi et al. (2007) Al-Saidi, W. A.; Zhang, S.; Krakauer, H. Bond breaking with auxiliary-field quantum Monte Carlo. _The Journal of chemical physics_ 2007, _127_ , 144101.
  * Suewattana et al. (2007) Suewattana, M.; Purwanto, W.; Zhang, S.; Krakauer, H.; Walter, E. J. Phaseless auxiliary-field quantum Monte Carlo calculations with plane waves and pseudopotentials: Applications to atoms and molecules. _Physical Review B_ 2007, _75_ , 245123.
  * Purwanto et al. (2015) Purwanto, W.; Zhang, S.; Krakauer, H. An auxiliary-field quantum Monte Carlo study of the chromium dimer. _J. Chem. Phys._ 2015, _142_ , 064302\.
  * Motta and Zhang (2018) Motta, M.; Zhang, S. Ab initio computations of molecular systems by the auxiliary-field quantum Monte Carlo method. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2018, _8_ , e1364.
  * Hao et al. (2018) Hao, H.; Shee, J.; Upadhyay, S.; Ataca, C.; Jordan, K. D.; Rubenstein, B. M. Accurate predictions of electron binding energies of dipole-bound anions via quantum Monte Carlo methods. _The journal of physical chemistry letters_ 2018, _9_ , 6185–6190.
  * Motta et al. (2019) Motta, M.; Shee, J.; Zhang, S.; Chan, G. K.-L. Efficient ab initio auxiliary-field quantum Monte Carlo calculations in Gaussian bases via low-rank tensor decomposition. _Journal of chemical theory and computation_ 2019, _15_ , 3510–3521.
  * Shee et al. (2019) Shee, J.; Arthur, E. J.; Zhang, S.; Reichman, D. R.; Friesner, R. A. Singlet–triplet energy gaps of organic biradicals and polyacenes with auxiliary-field quantum Monte Carlo. _Journal of chemical theory and computation_ 2019, _15_ , 4924–4932.
  * Shee et al. (2019) Shee, J.; Rudshteyn, B.; Arthur, E. J.; Zhang, S.; Reichman, D. R.; Friesner, R. A. On Achieving High Accuracy in Quantum Chemical Calculations of 3 d Transition Metal-Containing Systems: A Comparison of Auxiliary-Field Quantum Monte Carlo with Coupled Cluster, Density Functional Theory, and Experiment for Diatomic Molecules. _Journal of chemical theory and computation_ 2019, _15_ , 2346–2358.
  * Williams et al. (2020) Williams, K. T.; Yao, Y.; Li, J.; Chen, L.; Shi, H.; Motta, M.; Niu, C.; Ray, U.; Guo, S.; Anderson, R. J., et al. Direct comparison of many-body methods for realistic electronic Hamiltonians. _Physical Review X_ 2020, _10_ , 011041.
  * Lee and Reichman (2020) Lee, J.; Reichman, D. R. Stochastic resolution-of-the-identity auxiliary-field quantum Monte Carlo: Scaling reduction without overhead. _The Journal of Chemical Physics_ 2020, _153_ , 044131.
  * Lee et al. (2020) Lee, J.; Malone, F. D.; Morales, M. A. Utilizing essential symmetry breaking in auxiliary-field quantum Monte Carlo: Application to the spin gaps of the C36 fullerene and an iron porphyrin model complex. _Journal of chemical theory and computation_ 2020, _16_ , 3019–3027.
  * Shi and Zhang (2021) Shi, H.; Zhang, S. Some recent developments in auxiliary-field quantum Monte Carlo for real materials. _The Journal of Chemical Physics_ 2021, _154_ , 024107.
  * Loh Jr et al. (1990) Loh Jr, E.; Gubernatis, J.; Scalettar, R.; White, S.; Scalapino, D.; Sugar, R. Sign problem in the numerical simulation of many-electron systems. _Physical Review B_ 1990, _41_ , 9301.
  * Troyer and Wiese (2005) Troyer, M.; Wiese, U.-J. Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations. _Physical review letters_ 2005, _94_ , 170201.
  * Mahajan and Sharma (2021) Mahajan, A.; Sharma, S. Taming the Sign Problem in Auxiliary-Field Quantum Monte Carlo Using Accurate Wave Functions. _Journal of Chemical Theory and Computation_ 2021, _17_ , 4786–4798.
  * Giner et al. (2013) Giner, E.; Scemama, A.; Caffarel, M. Using perturbatively selected configuration interaction in quantum Monte Carlo calculations. _Canadian Journal of Chemistry_ 2013, _91_ , 879–885.
  * Evangelista (2014) Evangelista, F. A. Adaptive multiconfigurational wave functions. _The Journal of Chemical Physics_ 2014, _140_ , 124114.
  * Holmes et al. (2016) Holmes, A. A.; Tubman, N. M.; Umrigar, C. J. Heat-Bath Configuration Interaction: An Efficient Selected Configuration Interaction Algorithm Inspired by Heat-Bath Sampling. _J. Chem. Theory Comput._ 2016, _12_ , 3674–3680, PMID: 27428771.
  * Tubman et al. (2016) Tubman, N. M.; Lee, J.; Takeshita, T. Y.; Head-Gordon, M.; Whaley, K. B. A deterministic alternative to the full configuration interaction quantum Monte Carlo method. _The Journal of chemical physics_ 2016, _145_ , 044112.
  * Filippi et al. (2016) Filippi, C.; Assaraf, R.; Moroni, S. Simple formalism for efficient derivatives and multi-determinant expansions in quantum Monte Carlo. _The Journal of chemical physics_ 2016, _144_ , 194105.
  * Assaraf et al. (2017) Assaraf, R.; Moroni, S.; Filippi, C. Optimizing the Energy with Quantum Monte Carlo: A Lower Numerical Scaling for Jastrow–Slater Expansions. _J. Chem. Theory Comput._ 2017, _13_ , 5273–5281, PMID: 28873307\.
  * Dash et al. (2018) Dash, M.; Moroni, S.; Scemama, A.; Filippi, C. Perturbatively selected configuration-interaction wave functions for efficient geometry optimization in quantum Monte Carlo. _Journal of chemical theory and computation_ 2018, _14_ , 4176–4182.
  * Pineda Flores and Neuscamman (2019) Pineda Flores, S. D.; Neuscamman, E. Excited state specific multi-Slater Jastrow wave functions. _The Journal of Physical Chemistry A_ 2019, _123_ , 1487–1497.
  * Benali et al. (2020) Benali, A.; Gasperich, K.; Jordan, K. D.; Applencourt, T.; Luo, Y.; Bennett, M. C.; Krogel, J. T.; Shulenburger, L.; Kent, P. R.; Loos, P.-F., et al. Toward a systematic improvement of the fixed-node approximation in diffusion Monte Carlo for solids—A case study in diamond. _The Journal of Chemical Physics_ 2020, _153_ , 184111\.
  * Shee et al. (2018) Shee, J.; Arthur, E. J.; Zhang, S.; Reichman, D. R.; Friesner, R. A. Phaseless auxiliary-field quantum Monte Carlo on graphical processing units. _Journal of chemical theory and computation_ 2018, _14_ , 4109–4121.
  * Sharma et al. (2017) Sharma, S.; Holmes, A. A.; Jeanmairet, G.; Alavi, A.; Umrigar, C. J. Semistochastic Heat-Bath Configuration Interaction Method: Selected Configuration Interaction with Semistochastic Perturbation Theory. _J. Chem. Theory Comput._ 2017, _13_ , 1595–1604, PMID: 28263594\.
  * Smith et al. (2017) Smith, J. E.; Mussard, B.; Holmes, A. A.; Sharma, S. Cheap and near exact CASSCF with large active spaces. _J. Chem. Theory Comput._ 2017, _13_ , 5468–5478.
  * Mahajan and Sharma (2020) Mahajan, A.; Sharma, S. Efficient local energy evaluation for multi-Slater wave functions in orbital space quantum Monte Carlo. _The Journal of Chemical Physics_ 2020, _153_ , 194108.
  * Rom et al. (1997) Rom, N.; Charutz, D.; Neuhauser, D. Shifted-contour auxiliary-field Monte Carlo: circumventing the sign difficulty for electronic-structure calculations. _Chemical physics letters_ 1997, _270_ , 382–386.
  * Löwdin (1955) Löwdin, P.-O. Quantum theory of many-particle systems. I. Physical interpretations by means of density matrices, natural spin-orbitals, and convergence problems in the method of configurational interaction. _Physical Review_ 1955, _97_ , 1474.
  * Balian and Brezin (1969) Balian, R.; Brezin, E. Nonunitary Bogoliubov transformations and extension of Wick’s theorem. _Il Nuovo Cimento B (1965-1970)_ 1969, _64_ , 37–55.
  * Sun et al. (2018) Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S.; Wouters, S.; Chan, K.-L. G. PySCF: the Python-based simulations of chemistry framework. _WIREs Comput. Mol. Sci._ 2018, _8_ , e1340.
  * (45) https://github.com/sanshar/VMC/.
  * (46) https://github.com/ankit76/ph_afqmc.
  * Buonaura and Sorella (1998) Buonaura, M. C.; Sorella, S. Numerical study of the two-dimensional Heisenberg model using a Green function Monte Carlo technique with a fixed number of walkers. _Physical Review B_ 1998, _57_ , 11446.
  * Motta et al. (2017) Motta, M.; Ceperley, D. M.; Chan, G. K.-L.; Gomez, J. A.; Gull, E.; Guo, S.; Jiménez-Hoyos, C. A.; Lan, T. N.; Li, J.; Ma, F., et al. Towards the solution of the many-electron problem in real materials: Equation of state of the hydrogen chain with state-of-the-art many-body methods. _Physical Review X_ 2017, _7_ , 031059.
  * Trail and Needs (2017) Trail, J. R.; Needs, R. J. Shape and energy consistent pseudopotentials for correlated electron systems. _The Journal of chemical physics_ 2017, _146_ , 204107.
  * Whitlock et al. (1979) Whitlock, P.; Ceperley, D.; Chester, G.; Kalos, M. Properties of liquid and solid He 4. _Physical Review B_ 1979, _19_ , 5598.
  * Ceperley and Kalos (1986) Ceperley, D. M.; Kalos, M. _Monte Carlo methods in statistical physics_ ; Springer, 1986; pp 145–194.
  * Rothstein (2013) Rothstein, S. M. A survey on pure sampling in quantum Monte Carlo methods. _Canadian Journal of Chemistry_ 2013, _91_ , 505–510.
  * Motta and Zhang (2017) Motta, M.; Zhang, S. Computation of ground-state properties in molecular systems: Back-propagation with auxiliary-field quantum Monte Carlo. _Journal of chemical theory and computation_ 2017, _13_ , 5367–5378.
  * Muenter (1975) Muenter, J. Electric dipole moment of carbon monoxide. _Journal of Molecular Spectroscopy_ 1975, _55_ , 490–491.
  * Lide (2004) Lide, D. R. _CRC handbook of chemistry and physics_ ; CRC press, 2004; Vol. 85.
  * Huggins et al. (2021) Huggins, W. J.; O’Gorman, B. A.; Rubin, N. C.; Reichman, D. R.; Babbush, R.; Lee, J. Unbiasing Fermionic Quantum Monte Carlo with a Quantum Computer. _arXiv_ 2021,