ycode]python bgcolor=LightGray, fontsize=, breaklines=true

# Improved Modularity and New Features in ipie: Towards Even Larger AFQMC
Calculations on CPUs and GPUs at Zero and Finite Temperatures

Tong Jiang  Department of Chemistry and Chemical Biology, Harvard University,
Cambridge, MA, USA  Moritz K. A. Baumgarten  Department of Chemistry and
Chemical Biology, Harvard University, Cambridge, MA, USA  Pierre-François Loos
Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de
Toulouse, CNRS, UPS, France  Ankit Mahajan  Department of Chemistry, Columbia
University, New York, NY, USA  Anthony Scemama  Laboratoire de Chimie et
Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France  Shu
Fay Ung  Department of Chemistry, Columbia University, New York, NY, USA
Jinghong Zhang  Department of Chemistry and Chemical Biology, Harvard
University, Cambridge, MA, USA  Fionn D Malone  Google Research, Venice, CA
90291, United States  Joonho Lee<EMAIL_ADDRESS>Department of
Chemistry and Chemical Biology, Harvard University, Cambridge, MA, USA

###### Abstract

ipie is a Python-based auxiliary-field quantum Monte Carlo (AFQMC) package
that has undergone substantial improvements since its initial release [J.
Chem. Theory Comput., 2022, 19(1): 109-121]. This paper outlines the improved
modularity and new capabilities implemented in ipie. We highlight the ease of
incorporating different trial and walker types and the seamless integration of
ipie with external libraries. We enable distributed Hamiltonian simulations,
allowing for multi-GPU simulations of large systems. This development enabled
us to compute the interaction energy of a benzene dimer with 84 electrons and
1512 orbitals, which otherwise would not have fit on a single GPU. We also
support GPU-accelerated multi-slater determinant trial wavefunctions
[arXiv:2406.08314] to enable efficient and highly accurate simulations of
large-scale systems. This allows for near-exact ground state energies of
multi-reference clusters, [Cu2O2]2+ and [Fe2S2(SCH3)]2-. We also describe
implementations of free projection AFQMC, finite temperature AFQMC, AFQMC for
electron-phonon systems, and automatic differentiation in AFQMC for
calculating physical properties. These advancements position ipie as a leading
platform for AFQMC research in quantum chemistry, facilitating more complex
and ambitious computational method development and their applications.

##  I Introduction

Auxiliary-field quantum Monte Carlo (AFQMC) [1, 2] has become increasingly
popular in quantum chemistry [3, 4] and is furthermore being recognized as a
useful approach in the context of quantum algorithms. [5, 6, 7, 8, 9] A well-
maintained AFQMC program with flexibility and robust performance will play a
pivotal role at the intersection of many disciplines, including chemistry,
physics, materials science, and quantum information science.

Given the initial PAUXY development effort for rapid prototyping, while
lacking performance, [10, 11, 12, 13, 14] the Python-based AFQMC program ipie
was designed from scratch for high performance and ease of development. It was
officially introduced as a production-level package in Ref. 15. ipie was
optimized for high-performance computing architectures with both central and
graphical processing units (CPUs and GPUs). High performance was largely
achieved by integrating Numba’s JIT compilation [16] to fine-tune the
computational efficiency of specific kernels and MPI parallelism for effective
distributed computing. The utility of ipie was showcased through the
resolution of intricate quantum chemical challenges, notably the [Cu2O2]2+
torture track. [17, 15] Rigorous benchmarks on CPU and GPU platforms position
ipie competitively, displaying speed on par with–or surpassing–existing Python
and C++ codes. [15, 18, 19] Since its release, ipie has gained widespread
adoption in both quantum chemistry and quantum computing studies. [20, 21, 22,
3, 8, 7, 23]

In this article, we describe the recent development and current status of
ipie, introducing enhanced modularity, a suite of new features, interfaces to
external packages, and associated numerical examples. Below are the key
highlights:

##### High degree of modularity and customizability.

The AFQMC driver has been restructured to be fully modular, allowing for a
straightforward combination of features. This provides greater flexibility in
adapting to a wide range of user demands. Advancements in AFQMC algorithms
often focus on developing new trial wavefunctions to control better the
fermionic sign/phase problem. [5, 24, 25, 26, 23] The key routines in an AFQMC
calculation compute intermediates using walker and trial wavefunctions,
including overlap, force bias, Green’s function, and energy estimators.
Without altering the internal core code of ipie, users and developers can
customize all components making up an AFQMC simulation, including trial
wavefunctions, walkers, Hamiltonians, propagators, and estimators. The key
objects are all structured using object-oriented programming (OOP) principles,
facilitating straightforward customization through inheritance.

##### Development-friendly design.

This improvement standardizes the component interfaces and workflow processes
while offering a flexible system that adapts to various data types and user
requirements. The abstract base classes serve as a foundational blueprint,
ensuring all components adhere to a uniform structure and interact seamlessly.
Complementing this, the factory methods for common workflows simplify the
instantiation process, allowing users to set up standard calculations with
minimal effort and reduced potential for errors. Integrating type-based
dispatch through Plum [27] brings increased precision and efficiency in method
handling, ensuring every component downstream dynamically adjusts its
operations based on the specific trials and walkers requested.

Moreover, we have simplified the integration of ipie with external quantum
chemistry packages. While the necessary integrals and orbitals for running
AFQMC are most commonly obtained through an interface with PySCF, [28] our
simplified file format also ensures that other packages can be easily used.
For sophisticated trial wavefunctions, such as multiple Slater determinant
(MSD) trials derived from selected configuration interaction, interfaces with
PySCF, [28] Dice, [29, 19] and TrexIO[30] are available. An additional
interface with the Fermionic Quantum Emulator (FQE) [31, *rubin2021fermionic]
has also been introduced, facilitating the conversion between ipie’s MSD
wavefunction and quantum circuit wavefunctions, which facilitates AFQMC
applications in the quantum information science (QIS) community.

##### New improvements and features.

Several features have been added to ipie to improve memory management, handle
new problems, enhance integration testing, etc. To address the significant
memory demand associated with storing Cholesky vectors for large systems, ipie
offers options for shared memory across CPU processes or distributing among
CPU processes or GPU cards. Support for complex-valued Hamiltonians has also
been added.

While ipie is primarily developed for phaseless AFQMC (ph-AFQMC) targeting ab
initio systems, it supports other AFQMC methods as well, including free-
projection AFQMC, [24] finite temperature AFQMC, [12] AFQMC for coupled
electron-phonon systems, [13] and automatic differentiation within AFQMC for
calculating observables which do not commute with the Hamiltonian. [33]
Furthermore, ipie is equipped with improved integration testing, boosting the
package’s robustness and adaptability and significantly enhancing its
reliability and usability for end-users and developers alike.

The organization of this paper is as follows: Section II overviews the theory
of AFQMC; Section III details the components, software architecture, and
workflow of ipie, including examples to illustrate the framework’s
adaptability for AFQMC development; Section IV introduces new features in ipie
and provides corresponding examples; Section V outlines the interfaces to
external packages; and Section VI concludes with a summary and outlook.

##  II Theory of AFQMC

AFQMC is based on the following imaginary time evolution:

$\left|\Psi_{0}\right\rangle\propto\lim_{\tau\rightarrow\infty}\exp(-\tau\hat{H})\left|\Phi_{0}\right\rangle=\lim_{n\rightarrow\infty}(\exp(-\Delta\tau\hat{H}))^{n}|\Phi_{0}\rangle,$
(1)

where $\Delta\tau$ is an infinitesimal time step,
$\left|\Psi_{0}\right\rangle$ is the ground state wavefunction, and
$\left|\Phi_{0}\right\rangle$ is an initial state satisfying
$\langle\Phi_{0}|\Psi_{0}\rangle\neq 0$. While ipie supports some of the
prototypical model Hamiltonians, its development has focused on the simulation
of the ab initio Hamiltonian, which in second quantization is given by

$\hat{H}=\sum_{p,q=1}^{N}h_{pq}\hat{a}_{p}^{\dagger}\hat{a}_{q}+\frac{1}{2}\sum_{p,q,r,s=1}^{N}g_{psqr}\hat{a}_{p}^{\dagger}\hat{a}_{q}^{\dagger}\hat{a}_{r}\hat{a}_{s},$
(2)

where the two-electron repulsion integral (ERI) is factorized with the
Cholesky decomposition

$g_{psqr}=(ps|qr)=\sum_{\gamma=1}^{N_{\gamma}}L_{ps}^{\gamma}L_{qr}^{\gamma}.$
(3)

With this factorization, we have

$\hat{H}=\hat{v}_{0}-\frac{1}{2}\sum_{\gamma=1}^{N_{\gamma}}\hat{v}_{\gamma}^{2},$
(4)

where

$\hat{v}_{0}=\sum_{pq}\left[h_{pq}-\frac{1}{2}\sum_{r}(pr|rq)\right]a_{p}^{\dagger}a_{q}$
(5)
$\hat{v}_{\gamma}=\mathrm{i}\sum_{pq}L_{pq}^{\gamma}\hat{a}_{p}^{\dagger}\hat{a}_{q}.$
(6)

The short-time propagator with Trotter decomposition is written as

$\mathrm{e}^{-\Delta\tau\hat{H}}=\mathrm{e}^{-\frac{\Delta\tau}{2}\hat{v}_{0}}\mathrm{e}^{\frac{\Delta\tau}{2}\sum\hat{v}_{\gamma}^{2}}\mathrm{e}^{-\frac{\Delta\tau}{2}\hat{v}_{0}}+\mathcal{O}\left(\Delta\tau^{3}\right).$
(7)

Upon applying the Hubbard–Stratonovich transformation,[34, 35] our effective
propagator contains only one-body operators,

$\mathrm{e}^{-\Delta\tau\hat{H}}=\int\mathrm{d}\mathbf{x}\
p(\mathbf{x})\hat{B}(\mathbf{x},\Delta\tau)+\mathcal{O}\left(\Delta\tau^{2}\right),$
(8)

where $p(\mathbf{x})$ is the standard Gaussian distribution,
$\textbf{x}=(x_{1},x_{2},\cdots,x_{N_{\gamma}})$ are the auxiliary fields, and
the one-body propagator $\hat{B}$ is

$\hat{B}(\mathbf{x},\Delta\tau)=\mathrm{e}^{-\frac{\Delta\tau}{2}\hat{v}_{0}}\mathrm{e}^{-\sqrt{\Delta\tau}\mathbf{x}\cdot\hat{\mathbf{v}}}\mathrm{e}^{-\frac{\Delta\tau}{2}\hat{v}_{0}}.$
(9)

where
$\hat{\mathbf{v}}=(\hat{v}_{1},\hat{v}_{2},\cdots,\hat{v}_{N_{\gamma}})$.

In AFQMC, each walker samples auxiliary fields and represents a statistical
sample of the global wavefunction at imaginary time $\tau$, written as

$|\Psi(\tau)\rangle=\sum_{i}^{N_{\textrm{w}}}w_{i}(\tau)\frac{\left|\psi_{i}(\tau)\right\rangle}{\left\langle\Psi_{\textrm{T}}|\psi_{i}(\tau)\right\rangle},$
(10)

where $|\psi_{i}(\tau)\rangle$ is the wavefunction of the $i$-th walker at
time $\tau$ and $|\Psi_{\text{T}}\rangle$ is the trial wavefunction used for
importance sampling. The energy estimator is then

$E(\tau)=\sum_{i}^{N_{\textrm{w}}}w_{i}(\tau)E_{\textrm{loc},i}(\tau)=\sum_{i}^{N_{\textrm{w}}}w_{i}(\tau)\frac{\langle\Psi_{\textrm{T}}|\hat{H}|\psi_{i}(\tau)\rangle}{\left\langle\Psi_{\textrm{T}}|\psi_{i}(\tau)\right\rangle}.$
(11)

The walker state $\left|\psi_{i}(\tau)\right\rangle$ is updated by applying
the discrete-time propagator, and the weight $w_{i}(\tau)$ is updated
following the phaseless approximation: [2]

$\displaystyle\left|\psi_{i}(\tau+\Delta\tau)\right\rangle=\hat{B}\left(\mathbf{x}_{i}-\overline{\mathbf{x}}_{i},\Delta\tau\right)\left|\psi_{i}(\tau)\right\rangle$
(12) $\displaystyle
w_{i}(\tau+\Delta\tau)=I_{\mathrm{ph}}\left(\mathbf{x}_{i},\overline{\mathbf{x}}_{i},\tau,\Delta\tau\right)\times
w_{i}(\tau).$ (13)

We dynamically shift the distribution of auxiliary fields using the force bias
$\overline{\mathbf{x}}_{i}$ defined by

$\overline{\mathbf{x}}_{i}(\Delta\tau,\tau)=-\sqrt{\Delta\tau}\frac{\left\langle\Psi_{\textrm{T}}\left|\hat{\mathbf{v}}-\langle\hat{\mathbf{v}}\rangle_{\text{T}}\right|\psi_{i}(\tau)\right\rangle}{\left\langle\Psi_{\textrm{T}}|\psi_{i}(\tau)\right\rangle},$
(14)

where $\langle\hat{\mathbf{v}}\rangle_{\text{T}}$ is usually called the mean-
field shift in the AFQMC literature. [4] Furthermore, the phaseless importance
function [1, 2] used in the weight update (13) is given by

$I_{\mathrm{ph}}\left(\mathbf{x}_{i},\overline{\mathbf{x}}_{i},\tau,\Delta\tau\right)=\left|I\left(\mathbf{x}_{i},\overline{\mathbf{x}}_{i},\tau,\Delta\tau\right)\right|\times\max\left[0,\cos\left(\theta_{i}(\tau)\right)\right],$
(15)

which remains real and positive throughout the propagation. The so-called
hybrid importance function is given by

$I\left(\mathbf{x}_{i},\overline{\mathbf{x}}_{i},\tau,\Delta\tau\right)=S_{i}(\tau,\Delta\tau)\mathrm{e}^{\mathbf{x}_{i}\cdot\overline{\mathbf{x}}_{i}-\overline{\mathbf{x}}_{i}\cdot\overline{\mathbf{x}}_{i}/2},$
(16)

and the overlap ratio of the $i$-th walker is

$S_{i}(\tau,\Delta\tau)=\frac{\left\langle\Psi_{\textrm{T}}\left|\hat{B}\left(\Delta\tau,\mathbf{x}_{i}-\overline{\mathbf{x}}_{i}\right)\right|\psi_{i}(\tau)\right\rangle}{\left\langle\Psi_{\textrm{T}}|\psi_{i}(\tau)\right\rangle}.$
(17)

We define the phase of the overlap as

$\theta_{i}(\tau)=\arg\left[S_{i}(\tau,\Delta\tau)\right].$ (18)

The computation of local energies in Eq. 11, propagation of wavefunctions in
Eq. 12, and evaluation of force biases in Eq. 14 are the primary computational
hotspots in AFQMC calculations. We optimize these routines by employing
optimized algorithms and adopting high-performance computing techniques. We
provide an overview of the phaseless AFQMC algorithm described above in Fig.
1.

Figure 1: Overview of phaseless AFQMC in terms of computational hotspots.

##  III Software architecture and design principles

###  III.1 AFQMC driver

We overview the software architecture of ipie with emphasis on improved
modularity, as shown in Fig. 2.

Figure 2: The workflow of ipie.

ipie generates an AFQMC simulation by assembling the generic AFQMC driver from
problem-specific components and then running it to execute the simulation. The
driver requires several inputs to define the QMC simulation, including the
trial wavefunction, walkers, and the Hamiltonian, as well as more generic
parameters such as the number of blocks and time step. The following Python
Code Snippet shows how the driver can be instantiated:

⬇

class AFQMC(object):

def __init__(self, system, hamiltonian, trial,

walkers,propagator, mpi_handler, params:

QMCParams):

……

Figure 3: The AFQMC class within ipie.

The code takes several objects needed for a typical AFQMC driver:

  1. (1)

system contains information that defines the problem, including the number of
spin-up and spin-down electrons.

  2. (2)

hamiltonian stores molecular integrals as detailed in Sec. III.2.

  3. (3)

trial represents a trial wavefunction object, for which ipie provides many
options as detailed in Sec. III.3.

  4. (4)

walkers manages all information about walkers during the imaginary time
propagation. This includes the initial walker wavefunctions and the number of
walkers, among other details. ipie allows walkers to be explicitly passed into
the AFQMC driver or dynamically dispatched with default settings, as detailed
in Sec. III.4.

  5. (5)

propagator handles the update of walker weights and wavefunctions during the
imaginary time evolution. ipie’s AFQMC driver can take the propagator class as
an input or dispatch it internally with default settings.

  6. (6)

mpi_handler manages information related to the Message Passing Interface (MPI)
for parallel computing, including rank details, chunking, and groups (i.e., a
group is a collection of MPI processes). The parameter shared_comm in
mpi_handler refers to the MPI rank within a group where the chunking integrals
are distributed. The specifics of these are elaborated in Section IV.1.

  7. (7)

params encompasses other fundamental parameters for the QMC simulation, such
as the time step.

While users can directly provide all of the aforementioned inputs to construct
the AFQMC driver object, especially for development purposes, a factory method
is provided to simplify the process greatly. This is especially useful for
standard AFQMC calculations using single or multiple Slater determinant
trials. This method significantly reduces the number of inputs required to
construct the AFQMC object, as demonstrated in Code Snippet 4:

def build( num_elec: Tuple[int, int], hamiltonian, trial_wavefunction ) -¿
”AFQMC”: ……

Figure 4: Factory method for AFQMC driver.

As illustrated in the Code Snippet 4 and Fig. 2, the construction of the AFQMC
driver ultimately relies on two key inputs: the Hamiltonian and the trial
wavefunction. The Hamiltonian object is assembled with one-electron integrals
and the Cholesky decomposition (or density fitting) of two-electron integrals,
detailed in Sec. III.2. As for trial wavefunctions, ipie offers various
options that utilize interfaces with external packages such as PySCF, [28]
Dice, [19] TrexIO [30] to facilitate standardized workflows. ipie can also
accommodate user-customized trials, which offers flexibility in future
developments. Benefited from the flexibility of ipie, we are also able to
incorporate more complicated trial wavefunctions, such as density matrix
renormalization group (DMRG), and in Sec. III.3 we show the AFQMC calculation
with matrix product states (MPS) trials by incorporating ipie with the DMRG
package Renormalizer. [36, 37]

The build method in the AFQMC class is designed to streamline the setup of an
AFQMC calculation by requiring only three inputs: the number of electrons, the
Hamiltonian, and the trial wavefunction. Its main purpose is automating the
construction of key components such as walkers, propagators, and estimators,
etc.. Users need only specify the hamiltonian and trial inputs
trial_wavefunction, with the build method handling downstream instantiation of
the rest.

Once the driver is built, wavefunction propagation, stabilization, population
control, and estimator calculation are performed with the provided QMC
parameters. Simulation outputs in both text file and hdf5 file formats are
saved and updated after each block, which enables real-time reblocking
analysis [38] via built-in tools.

###  III.2 Hamiltonian

The construction of a Hamiltonian object can be achieved in several ways. One
can construct the object directly by inputting the one-electron integrals and
the Cholesky decomposition (Eq. (3)) of either 8-fold (for real symmetric
integrals) or 4-fold (for complex hermitian integrals) symmetric two-body
integrals:

from ipie.hamiltonians.generic import Generic as HamGeneric ham =
HamGeneric(h1e, chol, ecore)

Figure 5: Constructing the Hamiltonian object with provided electron integrals

where h1e, chol and ecore are the one-electron integrals, the Cholesky
decomposition (Eq. 6) of the two-electron integrals within the desired
orbitals in Eq. 2, and the nuclear repulsion energy, respectively.

ipie also offers built-in functions for assembling the Hamiltonian object from
PySCF calculations. One can provide: (i) a PySCF check file; (ii) a PySCF mol
object together with the molecular orbital (MO) coefficient matrix mo_coeff;
or (iii) a mol object, mo_coeff, and a basis transformation matrix X that
transforms MOs to orthogonal atomic orbitals (OAO), natural orbitals and so
on, depending on the single particle basis used in the AFQMC calculation:

from ipie.utils.from_pyscf import generate_hamiltonian,
generate_hamiltonian_from_chk

# OAO basis ham = generate_hamiltonian_from_chk( ’scf.chk’, use_mcscf=False,
chol_cut=1e-5, num_frozen_core=0, ortho_ao=True )

# MO basis ham = generate_hamiltonian(mol, mo_coeff, h1e, mo_coeff,
chol_cut=1e-5, num_frozen_core=0) # basis with transform matrix X ham =
generate_hamiltonian(mol, mo_coeff, h1e, X, chol_cut=1e-5)

Figure 6: Constructing the Hamiltonian object from PySCF.

The argument chol_cut specifies the cutoff for the Cholesky decomposition;
use_mcscf specifies whether to use the multi-configurational self-consistent
field (MCSCF) MO coefficients; and num_frozen_core specifies the number of
frozen cores for the subsequent AFQMC calculations.

Given that the size of Cholesky vectors scales as $\mathcal{O}(N^{3})$, with
$N$ as the number of orbitals, large systems with large basis sets can
generate electron integrals of considerable size. For such systems, one cannot
store copies of integrals for each MPI process and must resort to other
strategies. Utilizing shared memory across different MPI processes, ipie
offers an approach that allows processes on the same node to access Cholesky
vectors in the same memory block. The following code snippet can be adopted to
ensure the use of shared memory Hamiltonian objects:

from ipie.utils.mpi import get_shared_comm from ipie.hamiltonians.utils import
get_hamiltonian from mpi4py import MPI shared_comm =
get_shared_comm(MPI.COMM_WORLD) ham =
get_hamiltonian(”ham.h5”,comm=shared_comm, pack_chol=True)

Figure 7: Constructing the Hamiltonian object with shared memory.

This reads the integrals from the ham.h5 file and generates the Hamiltonian
object with shared integrals across all MPI processes. pack_chol refers to
another strategy used by ipie to reduce memory usage by a factor of two and
speed up the propagation. Namely, the permutation symmetry of Cholesky vectors
is exploited to utilize only the upper-triangular part of the Cholesky
vectors. Propagation using the packed matrix contracts only the symmetric part
in Eq. (12) and is thus more efficient.

However, the shared memory strategy often creates a memory access overhead and
reduces the MPI parallel efficiency when many cores are used. We also enable a
chunking (or distributed-memory) strategy. The Cholesky vectors are divided
into segments, or ‘chunked,’ and distributed across different MPI processes to
optimize memory usage and computational efficiency. The Cholesky vectors are
divided into several chunks, each of which will be allocated to separate MPI
processes. The specifics of this chunking process and the associated MPI
communication strategy will be thoroughly explained in Sec. IV.1. Given enough
nodes, our approach completely removes the memory bottleneck in AFQMC and can
work on both CPU and GPU nodes. We especially recommend our distributed-memory
implementation for GPUs, where memory capacity is much more limited than that
of CPUs.

###  III.3 Trial wavefunctions

Any new trial wavefunction may inherit from the TrialWavefunctionBase class,
which contains methods to calculate essential quantities including the overlap
(Eq. (17)), force bias (Eq. (14)), etc., as shown in Fig. 8. Some trials, such
as single/multiple Slater determinants, also require computation of Green’s
function and half rotation of electron repulsion integrals. [11]

Trial (WavefunctionBase) calc_overlap() calc_force_bias()
calc_greens_function()$*$ half_rotate()$*$
$\langle\Psi_{\textrm{T}}|\psi_{i}\rangle$ Eq. (14)
$\frac{\langle\Psi_{\textrm{T}}\left|\hat{a}_{p}^{\dagger}\hat{a}_{q}\right|\psi_{i}\rangle}{\langle\Psi_{\textrm{T}}\mid\psi_{i}\rangle}$
see Ref. 11 Figure 8: The functions to be defined for a trial object. $*$: The
function is optional depending on the trial type, used in certain trials to
computations like overlap, force bias, or local energy.

There is flexibility to either inherit these functions directly or overwrite
them with customized implementations, ensuring both consistency in fundamental
operations and adaptability for specialized needs.

Single Slater determinant (SD) trials can be generated with

from ipie.trial_wavefunction.single_det import SingleDet trial =
SingleDet(np.hstack([orbs_a, orbs_b]), nelec, num_basis)

Figure 9: The single determinant trial object.

where orbs_a and orbs_b are the spin-up and spin-down coefficient matrices for
the SD trials. Similarly, MSD trials can be generated with:

from ipie.trial_wavefunction.particle_hole import ParticleHole trial =
ParticleHole(wfn, nelec, num_basis)

Figure 10: The multiple Slater determinant trial object.

with wfn containing the occupation indices in the basis used in the
configuration interaction (CI) expansion of the MSD and their coefficients.

A new customized trial can be defined straightforwardly in ipie. For instance,
we consider a “noisy” SD trial that can be inherited from our existing SD
trial class. The noisy trial adds Gaussian noise to the overlap evaluation and
keeps everything else the same as for regular SD trials. The following code
snippet achieves this:

class NoisySingleDet(SingleDet): def __init__(self, wavefunction, num_elec,
num_basis, noise_level=1e-12): super().__init__(wavefunction, num_elec,
num_basis) self._noise_level = noise_level def calc_overlap(self, walkers) -¿
np.ndarray: ovlp = super().calc_overlap(walkers) noise = np.random.normal(
scale=self._noise_level, size=ovlp.size ) return ovlp * (1 + noise)

Figure 11: Example of constructing a custom noisy trial class.

Recently, we have also implemented a more sophisticated extension using a
matrix product state (MPS) trial wavefunction. [23] MPS is the variational
ansatz used in the density matrix renormalization group (DMRG) algorithm, and
we obtain the MPS solution from external DMRG packages such as Renormalizer.
[36, *ren2022time] The theory explaining the structure of the MPS trial and
associated functions is detailed in Ref. 23. A representative calculation is
shown in Fig. 12.

Figure 12: Application of MPS-AFQMC to two-dimensional hydrogen lattice with
$r=4.2a_{0}$. Energy deviation from exact results using MPS-AFQMC, UHF-AFQMC,
and DMRG with different bond dimensions. [23]

###  III.4 Dispatchers of walkers, propagators and estimators

As the program grows, additional complexity arises from the need to dispatch
different features for arbitrary combinations of smaller building blocks.
While there are numerous ways to achieve this, ipie previously relied on
conditional if statements at the computational subroutines, which–despite its
simplicity–will result in many intractable conditions as the number of
features increases. The latest release of ipie thus utilizes the multiple
dispatch library Plum [27] to dispatch features for different combinations of
objects in a compact and flexible manner. This can be especially useful in
AFQMC, where the behavior of many functions varies based on the type of trial,
walkers, and Hamiltonians. In Fig. 13, we illustrate the dispatcher for
walkers depending on the type of trial. As for ipie’s native trial
wavefunctions, including both single and multiple (orthogonal and non-
orthogonal) Slater determinants, once the trial and Hamiltonian objects are in
place, the AFQMC driver can be formed. In instances where the walkers object
is not explicitly provided, the AFQMC driver constructs it inside the build
function based on the type of the trial wavefunction via Plum. The appropriate
energy estimators and propagators are also built according to the object type
of the trial, walkers, and Hamiltonian, completing the necessary
initialization for AFQMC calculations.

Plum
DispatcherBaseWalkersGHFWalkersUHFWalkersUHFWalkersParticleHoleUHFWalkersNOCIParticleHoleSingleDetNOCISingleDetGHFGHFWalkersTrialUHFWalkersTrialWavefunctionBase

Figure 13: Dispatcher for different walkers depending on the trial. Lines with
arrows indicate class inheritance, while dashed lines signify the use of the
Plum dispatcher for dispatching.

###  III.5 Estimators

ipie supports well-optimized energy estimators for wavefunctions, including
SD, [11] MSD, [39] and non-orthogonal configuration interaction (NOCI) trials.
It supports real and complex Hamiltonians on both CPUs and GPUs. As mentioned
in Section. III.4, dispatchers can automatically determine a specific
implementation of the estimator required for given user inputs. Users can also
compose customized estimators for new types of trial wavefunctions.
Considering the example in Code Snippet 11, we may define the corresponding
energy estimator as

class NoisyEnergyEstimator(EnergyEstimator): def __init__( self, system, ham,
trial): super().__init__(system=system, ham=ham, trial=trial) def
compute_estimator(self, system, walkers, hamiltonian, trial, istep=1):
trial.calc_greens_function(walkers) energy = local_energy_batch(system,
hamiltonian, walkers, trial) self._data[”ENumer”] = xp.sum(walkers.weight *
energy[:, 0].real) self._data[”EDenom”] = xp.sum(walkers.weight)
self._data[”E1Body”] = xp.sum(walkers.weight * energy[:, 1].real)
self._data[”E2Body”] = xp.sum(walkers.weight * energy[:, 2].real) return
self.data

Figure 14: Customized noisy estimator.

where xp serves as an abstract linear algebra backend, facilitating seamless
integration of either NumPy for pure CPU-based computations or CuPy for GPU-
accelerated tasks. In this estimator, we still evaluate the local energy as we
do for any SD trial. This customized estimator can be passed to the AFQMC
driver simply with

add_est = ”energy”: NoisyEnergyEstimator(Generic(mol.nelec), ham, trial)
afqmc.run(additional_estimators=add_est)

Figure 15: Adding the customized noisy estimator to the AFQMC driver.

For more sophisticated cases such as MPS trials, in addition to the usual
energy estimator, we can furthermore define other estimators such as one to
measure the bond dimension of walkers when converted into an MPS. [23]

##  IV New Features

Many new features have been added to ipie since our first release paper; [15]
we highlight some of the more important ones in this section.

###  IV.1 Distributed Hamiltonian for limited memory

When using GPUs, we routinely encounter cases where the Cholesky vectors
cannot entirely fit into a single GPU card’s memory. Or one could imagine a
large system where a single CPU node cannot fit the Cholesky vectors in memory
anymore. In such circumstances, the Cholesky vectors are distributed over
several GPU cards or CPU nodes that constitute a group. Each group stores the
full Cholesky tensor as shown in Fig. 16(a). Whenever the tensor is required
by the computational subroutines, for example, in computing the force bias or
local energy, the walkers in each member of the group are communicated to
other members in the same group, as illustrated in Fig. 16(b). We employ a
cyclic data passing scheme of walkers among MPI processes to apply all
existing Cholesky chunks to given walkers. With this strategy, one can
accumulate locally contracted quantities before moving to the next QMC time
step.

Figure 16: (a) Distribution of walkers and the chunked Hamiltonian over
multiple groups. (b) MPI communication between MPI processes within a group
with the distributed Hamiltonian and walker batches.

We performed timing benchmarks as an application by computing the parallel-
displaced benzene dimer, a Van der Waals complex taken from the S22 dataset.
[40] Disk usage for the benzene dimer calculations varies with the choice of
basis set: 11GB for the aug-cc-pVTZ basis and 67GB for the aug-cc-pVQZ basis
with a Cholesky decomposition threshold of $10^{-3}$. The threshold is set
based on the Hartree-Fock energy error around typical density fitting errors
[41] ($<50\mu E_{\textrm{h}}$ per atom here). We used our shared memory
framework for CPU calculations and distributed Cholesky vectors and half-
rotated integrals [11] across multiple cards were used for GPU computations.
For the benzene dimer analysis employing an aug-cc-pVTZ basis set, a single
A100 GPU card achieves performance around tenfold faster than 36 CPU cores
(see Fig. 17). Increasing the GPU count for non-chunked Hamiltonian
configurations (nmembers=1) results in a nearly linear acceleration
proportional to the number of GPUs deployed, as shown in Fig. 17(a).
Distributing the Hamiltonian across additional GPUs slightly slows down the
performance, predominantly during the VHS step, which involves the
construction of two-body propagators (where HS stands for Hubbard-
Stratonovich) during propagation.

Figure 17: (a) Timing benchmarks for AFQMC calculations on the benzene dimer
with the aug-cc-pVTZ basis employing 36 CPUs with shared memory and multiple
A100 GPUs with the distributed Hamiltonian. The data were obtained by
averaging 100 blocks, each containing 25 steps with a time step of
0.005$E_{\textrm{h}}^{-1}$. The interaction energy of the dimer computed from
AFQMC is compared against other methods using the aug-cc-pVTZ basis in (b) and
in the complete basis set limit obtained from a two-point basis extrapolation
[42] from aug-cc-pVTZ and aug-cc-pVQZ calculations in (c). The MP2, CCSD, and
CCSD(T) data using the aug-cc-pVTZ basis in (b) were obtained from Refs. 43
and 44. The CCSD(T) and FN-DMC data in (c) were extracted from Ref. 45. All
AFQMC results were obtained with 1224 total walkers.

Using our GPU implementation, we calculated the interaction energy of the
benzene dimer with counterpoise correction. [46] We compared AFQMC results
with those from other methodologies, as shown in Fig. 17(b)(c). While coupled
cluster with singles and doubles (CCSD) underestimates the binding energy,
second-order Møller-Plesset perturbation theory overestimates. Crucially, we
observed that RHF-AFQMC is on par with CCSD with perturbative triples
(CCSD(T)).

###  IV.2 Support of complex Cholesky vectors

Most calculations in quantum chemistry deal with real-valued integrals.
However, in some cases, especially in the presence of external magnetic
fields, when considering spin-orbit coupling or employing Bloch orbitals in
solid-state calculations, the wavefunctions (and orbitals) can become complex-
valued. [47] This adds a layer of mathematical and computational complexity.
To handle these cases, ipie supports complex Cholesky vectors; the number of
auxiliary fields is hence twice as many as the number of Cholesky vectors. The
implementation details are available in Ref. 48.

###  IV.3 GPU accelerated MSD-AFQMC and timing benchmarks

ph-AFQMC with MSD trials is a powerful tool for systems with multi-reference
characteristics, especially when the phaseless bias can be converged away.
Building on top of our successful CPU implementation [15] based on Wick’s
theorem, [39] As detailed in Ref. 49, ipie has now enabled GPU-accelerated
MSD-AFQMC implementations. Here, we summarize some of the important features
of our implementation and discuss applications of this feature. Interested
readers are referred to Ref. 49 for more details.

Similar to our single determinant GPU implementation, ipie employs cupy.einsum
via the cuTENSOR library to enhance matrix multiplication performance.
Additionally, Wick’s algorithm is adapted for GPU execution using custom CUDA
kernels, effectively speeding up the computational hotspots in MSD-AFQMC, such
as Green’s function, overlap, and local energy computations.

The GPU implementation of MSD-AFQMC is benchmarked on the example [Cu2O2]2+
from ipie’s first release paper. [15] In addition to employing a total of 10
walkers to directly compare with our previous results in Ref. 15, calculations
using a total of 640 walkers were also performed to obtain a more realistic
comparison. With 10 walkers, the GPU implementation on a single NVIDIA A100
card is six times faster than using a single CPU core when the number of
determinants in the trial is less than $10^{3}$. As the number of determinants
increases, the performance improvement is even more pronounced: tenfold with
$10^{4}$ determinants and approximately 100-fold with $10^{6}$ determinants,
as illustrated in Fig. 18(a). With 640 walkers, the GPU code on a single A100
card is compared against with 32 CPU cores. As shown in Fig. 18(b), a fourfold
speedup is observed with less than $10^{2}$ determinants that increase to
around sixfold with more determinants. The determinants were divided into
chunks and computed sequentially to handle a larger number of determinants
within the memory limits. This approach introduces an overhead that mostly
affects the energy estimator, as reported in Ref. 49.

Figure 18: (a) Comparisons of the time cost per block for AFQMC calculations
on [Cu2O2]2+ using the BS1 basis [15] with 10 walkers on a single CPU core and
a single A100 GPU. (b) Time per block with 640 walkers on a single A100
against 32 CPU cores. Figure adapted from Ref. 49.

The timing and absolute energy benchmarks are further benchmarked on the (20o,
30e) active space of the [Fe2S2(SCH3)]2- cluster, [49] as shown in Fig. 19. In
contrast to the 108-orbital calculation of [Cu2O2]2+, a merely two-fold
speedup is observed using a single A100 card compared to 32 CPUs. This speedup
increases to almost tenfold when the number of determinants exceeds $10^{4}$.
As shown in Fig. 19(b), another observation is using natural orbitals for
AFQMC calculation leads to chemical accuracy with around 3$\times 10^{5}$
determinants in the trial. In contrast, reaching chemical accuracy is more
challenging when using localized orbitals employed for DMRG calculations. [50]
These observations highlight the critical need to select an appropriate set of
orbitals for AFQMC calculations involving strongly correlated systems.
Additionally, it is noteworthy that the energies obtained using localized
orbitals with a limited number of determinants are initially lower than the
FCI energy, and they tend to converge upward toward the FCI reference as the
number of determinants increases.

Figure 19: Time cost and absolute energy benchmarks on the [Fe2S2(SCH3)]2-
cluster. (a) The time cost per block with 640 walkers on 32 CPU cores and a
single A100 GPU. (b) Absolute energies are derived using localized atomic
orbitals and natural orbitals. Figure adapted from Ref. 49.

###  IV.4 Free projection AFQMC

Free projection (fp-) AFQMC is a numerically exact method for calculating the
eigenvalues of a Hamiltonian. While the more commonly used ph-AFQMC variant
employs the phaseless constraint to control the sign problem, fp-AFQMC
attempts to sample the ground state energy brute-force with exponential-
scaling sample complexity. Despite this scaling, it is possible to perform
relatively large active space calculations in practice using accurate trial
states.

In fp-AFQMC, the ground-state energy is estimated by

$E(\tau)=\frac{\langle\psi_{l}|\hat{H}\mathrm{e}^{-\tau\hat{H}}|\psi_{r}\rangle}{\langle\psi_{l}|\mathrm{e}^{-\tau\hat{H}}|\psi_{r}\rangle}=\frac{\int
d\mathbf{x}\
p(\mathbf{x})\langle\psi_{l}|\hat{H}\hat{B}(\mathbf{x},\Delta\tau)|\psi_{r}\rangle}{\int
d\mathbf{x}\
p(\mathbf{x})\langle\psi_{l}|\hat{B}(\mathbf{x},\Delta\tau)|\psi_{r}\rangle},$
(19)

where $\ket{\psi_{l}}$ and $\ket{\psi_{r}}$ are left and right trial states,
respectively; $\hat{B}(\mathbf{x},\Delta\tau)$ is defined in Eq. 9; and $\tau$
is the imaginary time. Analogously to the ph-AFQMC propagator sampling, the
auxiliary fields $\mathbf{x}$ are sampled from the Gaussian distribution
$p(\mathbf{x})$. Unlike ph-AFQMC, we do not use the force bias to perform
importance sampling. Instead, only mean-field subtraction is used. We also do
not employ population control.

A judicious choice of $\ket{\psi_{l}}$ and $\ket{\psi_{r}}$ offers two
advantages. First, more accurate trial states (eg., selected CI) reduce the
imaginary time required to project the ground state energy to a given
accuracy, thereby reducing the noise in the energy estimate. [24, 39]
Integrating the GPU-accelerated MSD code [49] will also significantly
accelerate the fp-AFQMC calculations with MSD trials. Second, since energies
are measured using the state $\ket{\psi_{l}}$, the closer this state is to the
ground state, the smaller the variance in the energy estimate due to the zero
variance principle. [51] Furthermore, one can use a CCSD wavefunction as the
initial state $\ket{\psi_{r}}$, which reduces the projection time. CCSD
wavefunction is employed by performing a Hubbard–Stratonovich transformation
on the exponential of the cluster operator. [24] Using trial states belonging
to specific symmetry sectors also allows one to target the lowest energy
states in the corresponding sectors.

As an illustrative example, we consider the $D_{4h}$ symmetric transition
state of cyclobutadiene. This state has a biradical character, which makes it
a challenging problem for single reference electronic structure methods,
including spin-restricted and spin-unrestricted CCSD(T) (RCCSD(T) and
UCCSD(T)). We performed fp-AFQMC calculations on this system using an MSD
trial, obtained via heat-bath configuration interaction (HCI), as
$|\psi_{l}\rangle$ and the spin-restricted CCSD (RCCSD) state as
$|\psi_{r}\rangle$. The HCI state was obtained from an HCI calculation with a
crude $\epsilon_{1}=10^{-4}$ in the full space except for four frozen HF
orbitals. These orbitals were kept frozen in all correlated calculations. The
geometry was taken from Ref. 24, which also reported fp-AFQMC energies for
this system. As additional validation of our results, we converged the ph-
AFQMC energy with respect to the number of determinants in the trial HCI
state; we expect this energy to be nearly exact. Results are shown in Fig. 20.

The fp-AFQMC energy converges to the converged ph-AFQMC energy within
statistical error bars. Our results suggest that hybrid (H)-AFQMC energies
reported in Ref. 52 are likely biased as they are too much lower than our
converged ph-AFQMC and fp-AFQMC energies. RCCSD(T) is about 13.5
m$E_{\textrm{h}}$ higher than ph-AFQMC, while UCCSD(T) is closer but still
about 4.5 m$E_{\textrm{h}}$ higher. The differences in the barrier height
largely come from transition state energies because all these methods work
well for the equilibrium $D_{2h}$ geometry.

Figure 20: Convergence of the fp-AFQMC energy for the cyclobutadiene
transition state in the cc-pVDZ basis set with projection time. The inset
shows the convergence of ph-AFQMC energy with respect to the number of
determinants in the trial.

###  IV.5 AFQMC beyond the ground-state electronic structure energy

The preceding discussions concentrate on the ground-state ab initio electronic
structure energy calculations. Apart from this, ipie now also accommodates
other types of calculations within the framework of AFQMC, such as property
calculations, [33] finite temperature calculations [53, 12] and ground-state
calculations of electron-phonon coupled model systems. [13] The code
structures for these different features mirror those of the zero temperature
ab initio electronic structure implementation.

####  IV.5.1 Electronic structure at finite temperatures

An extension of the phaseless AFQMC method at zero temperature, finite-
temperature AFQMC (FT-AFQMC), was developed to study systems at finite
temperatures. [12] It is customary (though not necessary [54]) to work in the
grand canonical ensemble described by the temperature $T$, volume $V$, and
chemical potential $\mu$. The central quantities of interest are thermal
expectation values computed from the partition function

$\displaystyle\Xi=\mathrm{Tr}\left[\mathrm{e}^{-\beta(\hat{H}-\mu\hat{N})}\right],$
(20)

where $\hat{N}$ is the total number operator and $\mu$ is the chemical
potential. The Boltzmann factor in Eq. (20) can be interpreted as a propagator
in imaginary time $\tau=\beta$.

Analogously to the zero temperature case, $\tau$ is first discretized into $l$
intervals of length $\Delta\tau=\tau/l$,

$\displaystyle\Xi=\mathrm{Tr}\left[\mathrm{e}^{-\beta(\hat{H}-\mu\hat{N})}\right]=\mathrm{Tr}\left[\lim_{l\to\infty}\prod_{k=1}^{k=l}\mathrm{e}^{-\Delta\tau(\hat{H}-\mu\hat{N})}\right],$
(21)

the short-time propagator is then Trotter decomposed as

$\displaystyle\mathrm{e}^{-\Delta\tau(\hat{H}-\mu\hat{N})}\simeq\mathrm{e}^{-\frac{\Delta\tau}{2}(\hat{v}_{0}-\mu\hat{N})}\mathrm{e}^{\frac{\Delta\tau}{2}\sum_{\gamma}\hat{v}_{\gamma}^{2}}\mathrm{e}^{-\frac{\Delta\tau}{2}(\hat{v}_{0}-\mu\hat{N})},$
(22)

and the application of the Hubbard–Stratonovich transformation gives

$\displaystyle\mathrm{e}^{-\Delta\tau(\hat{H}-\mu\hat{N})}$
$\displaystyle\simeq\int d\mathbf{x}\
p(\mathbf{x})\hat{B}(\mathbf{x},\Delta\tau,\mu),$ (23)

where the one-body propagator $\hat{B}$ is now also function of $\mu$:

$\hat{B}(\mathbf{x},\Delta\tau,\mu)=\mathrm{e}^{-\frac{\Delta\tau}{2}(\hat{v}_{0}-\mu\hat{N})}\mathrm{e}^{-\sqrt{\Delta\tau}\mathbf{x}\cdot\hat{\mathbf{v}}}\mathrm{e}^{-\frac{\Delta\tau}{2}(\hat{v}_{0}-\mu\hat{N})}.$
(24)

The grand canonical partition function is thus evaluated as

$\displaystyle\Xi$
$\displaystyle=\mathrm{Tr}\left[\mathrm{e}^{-\beta(\hat{H}-\mu\hat{N})}\right]$
$\displaystyle=\int d\mathbf{x}_{1}\cdots d\mathbf{x}_{l}\
\underbrace{p(\mathbf{x}_{1})\cdots
p(\mathbf{x}_{l})}_{p(\mathbf{x}_{1},\dots,\mathbf{x}_{l})}\
\mathrm{Tr}\left[\prod_{k=1}^{l}\hat{B}(\mathbf{x}_{k},\Delta\tau,\mu)\right],$
(25)

with $p(\mathbf{x}_{1},\dots,\mathbf{x}_{l})$ being the probability of
sampling a specific path designated by auxiliary fields
$\mathbf{x}_{1},\dots,\mathbf{x}_{l}$. Furthermore, the trace in Eq. (25) can
be written in terms of a determinant, [55, 56, 57] which finally yields

$\displaystyle\Xi$ $\displaystyle=\int d\mathbf{x}_{1}\cdots d\mathbf{x}_{l}\
p(\mathbf{x}_{1},\dots,\mathbf{x}_{l})\
\mathrm{det}\left[\mathbf{I}+\prod_{k=1}^{l}\mathbf{B}(\mathbf{x}_{k},\Delta\tau,\mu)\right],$
(26)

where $\mathbf{I}$ is the identity matrix and $\mathbf{B}$ is a matrix
representation of $\hat{B}$ in a single-particle basis. It should be noted
that the partition function itself is not explicitly calculated–expectation
values derived from it are. For some generic observable $\hat{O}$, its
expectation value is

$\displaystyle\braket{\hat{O}}$
$\displaystyle=\frac{1}{\Xi}\mathrm{Tr}\left[\mathrm{e}^{-\beta(\hat{H}-\mu\hat{N})}\hat{O}\right]$
(27) $\displaystyle=\int d\mathbf{X}\ \frac{p(\mathbf{X})\
\mathrm{Tr}\left[A(\mathbf{X},\Delta\tau,\mu)\right]}{\int d\mathbf{Y}\
p(\mathbf{Y})\ \mathrm{Tr}\left[A(\mathbf{Y},\Delta\tau,\mu)\right]}\times$
(28)
$\displaystyle\frac{\mathrm{Tr}\left[A(\mathbf{X},\Delta\tau,\mu)\hat{O}\right]}{\mathrm{Tr}\left[A(\mathbf{X},\Delta\tau,\mu)\right]},$
(29)

where we introduced the shorthand $\mathbf{X}$ for the set of auxiliary fields
along an imaginary path
$\left\\{\mathbf{x}_{1},\cdots,\mathbf{x}_{l}\right\\}$, and

$\displaystyle A(\mathbf{X},\Delta\tau,\mu)$
$\displaystyle=\prod_{k=1}^{l}\hat{B}(\mathbf{x}_{k},\Delta\tau,\mu).$ (30)

Rewriting Eq. (27) in the form Eq. (29) allows $\braket{\hat{O}}$ to be
estimated through an importance sampling procedure. The field configurations
$\mathbf{X}$ are obtained via Monte Carlo sampling from the modified
probability distribution $\tilde{p}(\mathbf{X})=p(\mathbf{X})\
\mathrm{Tr}\left[A(\mathbf{X},\Delta\tau,\mu)\right]$, while the computed
random variables are the local expectation values

$\displaystyle
O_{L}(\mathbf{X},\Delta\tau,\mu)=\frac{\mathrm{Tr}\left[A(\mathbf{X},\Delta\tau,\mu)\hat{O}\right]}{\mathrm{Tr}\left[A(\mathbf{X},\Delta\tau,\mu)\right]}.$
(31)

Similar to the zero temperature case, the importance sampling is implemented
by initializing a set of walkers in the space of auxiliary fields with weights
$w_{i}$ and propagating them in imaginary time. We therefore evaluate Eq. (29)
in practice as

$\displaystyle\braket{\hat{O}}=\frac{\sum_{i}w_{i}O_{L}(\mathbf{X}_{i})}{\sum_{i}w_{i}}.$
(32)

An example of the imaginary time trace for observables is provided in Fig. 21,
which reproduces Fig. 3 in Ref. 12.

Figure 21: Imaginary time trace from FT-AFQMC of the uniform electron gas
model with 2 electrons. The top panel depicts the average electron number at
imaginary time $\tau$, $\braket{\hat{N}}_{\tau}$, while the bottom panel shows
the error in the total energy compared to FCI, i.e. $E_{\text{FT-
AFQMC}}(\tau)-E_{\text{FCI}}$. These results reproduce Fig. 3 in Ref. 12.

####  IV.5.2 Electrons coupled to phonons

Utilizing the tools provided in ipie, one may write a new projector Monte
Carlo method to calculate the ground state of electrons coupled to phonons
(i.e., lattice vibration).[13]

$\hat{H}=\underbrace{\sum_{pq}h_{pq}a_{p}^{\dagger}a_{q}}_{\hat{H}_{\text{el}}}+\underbrace{\sum_{\nu}\omega_{\nu}b_{\nu}^{\dagger}b_{\nu}}_{\hat{H}_{\text{ph}}}+\underbrace{\sum_{pq\nu}g_{pq\nu}a_{p}^{\dagger}a_{q}(b_{\nu}+b_{\nu}^{\dagger})}_{\hat{H}_{\text{el-
ph}}}$ (33)

where $b^{(\dagger)}$ represents the bosonic annihilation (and creation)
operator, the first term represents the electronic band term, the second term
represents the phonon band term, and the third term represents the coupling
between electrons and phonons.

The corresponding ground-state projector is obtained from trotterizing the
imaginary time propagator corresponding to Eq. (33),

$\displaystyle\mathrm{e}^{-\Delta\tau\hat{H}}$
$\displaystyle=\mathrm{e}^{-\frac{\Delta\tau}{2}\hat{H}_{\mathrm{el}}}\mathrm{e}^{-\frac{\Delta\tau}{2}\hat{H}_{\mathrm{ph}}}\mathrm{e}^{-\Delta\tau\hat{H}_{\mathrm{el-
ph}}}\mathrm{e}^{-\frac{\Delta\tau}{2}\hat{H}_{\mathrm{ph}}}\mathrm{e}^{-\frac{\Delta\tau}{2}\hat{H}_{\mathrm{el}}}$
$\displaystyle+\mathcal{O}(\Delta\tau^{3})$ (34)

In our method, we choose walkers of the form
$|\psi_{\mathrm{w}}\rangle\otimes|\mathbf{X}_{\mathrm{w}}\rangle\equiv|\psi_{\mathrm{w}}(\tau),\mathbf{X}_{\mathrm{w}}(\tau)\rangle$,
with $|\psi_{\mathrm{w}}\rangle$ being a single determinant and
$|\mathbf{X}_{\mathrm{w}}\rangle$ being the coordinates for phonon
displacements. Using the Monte Carlo sampling, we can work in the position
space of phonons without invoking any boson number truncation. This strategy
differs from other standard Monte Carlo approaches in this area. [58, 59]

With importance sampling via a trial wavefunction, $|\Psi_{\text{T}}\rangle$,
our global wavefunction is a weighted linear combination of walker vibronic
wavefunctions,

$|\Psi(\tau)\rangle=\sum_{\mathrm{w}=1}^{N_{\text{walkers}}}w_{\mathrm{w}}(\tau)\>\frac{|\psi_{\mathrm{w}}(\tau),\mathbf{X}_{\mathrm{w}}(\tau)\rangle}{\langle\Psi_{\text{T}}|\psi_{\mathrm{w}}(\tau),\mathbf{X}_{\mathrm{w}}(\tau)\rangle}.$
(35)

The propagation of the bosonic degrees of freedom under
$\hat{H}_{\mathrm{ph}}$ is performed via a diffusion Monte Carlo algorithm.
[60] With the walker representation, the propagation with
$\hat{H}_{\text{ph}}$ and $\hat{H}_{\text{el-ph}}$ is straightforward by
exploiting the Thouless theorem. More details can be found in Ref. 13. As an
example calculation, we picked the one-dimensional Holstein model under a
periodic boundary condition, [61] which involves $h_{pq}=-t\delta_{p+1,q}$,
$\omega_{\nu}=\omega$, and $g_{pq\nu}=g\delta_{pq}\delta_{p\nu}$. For 20-site
at half-filling, we considered four different unitless electron-phonon
coupling strengths $\lambda=\frac{g^{2}}{2t\omega}$, which are displayed in
Fig. 22.

Figure 22: Reproduction of Fig. 3 in Ref. 13. Results are shown for a 20-site
Holstein model at half-filling. We used a coherent state trial to obtain the
AFQMC/S results. CSMP2 data points correspond to perturbation theory-based
computations based on molecular orbitals obtained from SCF calculations for a
Holstein Fock operator.[13] Panels (a)-(d) show scans along phonon frequencies
$\omega$ for $\lambda$ equals (a) 0.1, (b) 0.3, (c) 0.8, and (d) 2.

####  IV.5.3 Automatic differentiatiable AFQMC

The computation of observables that do not commute with $\hat{H}$ poses
additional challenges in projector Monte Carlo, such as AFQMC. A recently
proposed approach by Mahajan et. al.[33] aims to compute the response
estimator (Eq. (36)), which, though still not exact due to the discontinuity
in the distribution, is accurate enough to give reliable results:

$\langle\hat{O}\rangle_{\text{response}}=\frac{\mathrm{d}E_{\text{AFQMC}}(\lambda)}{\mathrm{d}\lambda}\Bigg{|}_{\lambda=0}.$
(36)

Here $\lambda$ is a parameter for the perturbed Hamiltonian:
$\hat{H}(\lambda)=\hat{H}+\lambda\hat{O}$, and AFQMC energy
$E_{\text{AFQMC}}(\lambda)$ is given by the Monte Carlo estimator

$E_{\text{AFQMC}}(\lambda)=\frac{\sum_{i}w_{i}(\lambda)\frac{\langle\Psi_{T}(\lambda)|\hat{H}(\lambda)|\psi_{i}(\lambda)\rangle}{\langle\Psi_{T}(\lambda)|\psi_{i}(\lambda)\rangle}}{\sum_{i}w_{i}(\lambda)},$
(37)

where $w_{i}(\lambda)$ is the weight for the $i$-th walker, and
$|\psi_{i}(\lambda)\rangle$ is the $i$-th walker state at coupling strength
$\lambda$.

We implemented this scheme for computing observables using the Automatic
Differentiation (AD) functionality of PyTorch[62] as an add-on within ipie. It
is noteworthy that storing the computation graph for an entire AFQMC run
requires substantial memory. In practice, we use the concept of AD blocks.[33]
We only track the computation graph within an AD block; there is no connection
between different AD blocks. Using this approach, we manage memory costs by
adjusting block size. We also use the gradient checkpointing technique[63] to
reduce the memory cost further.

For relatively large systems, the AFQMC calculation is parallelized using MPI
to distribute walkers over MPI tasks. The differentiation of the AFQMC global
energy estimator is not embarrassingly parallel since walkers will mix between
different MPI ranks via population control. The AD implementation in ipie does
not support communication between MPI tasks. Therefore, we perform “local”
population control within an MPI rank. Because the local population control
does not mix the walkers between different ranks, it is valid to regard those
estimates as independent samples. We thus perform block analysis on those
samples to obtain the final result. Effectively, this amounts to running very
low walker population simulations (50 walkers per task here), which is only
practically possible for small system sizes where population control is not a
concern. For intermediate system sizes, one could imagine using large memory
nodes and OpenMP threading or using single GPUs with a sufficient amount of
memory.

Here, we present the results of AD-AFQMC calculations on various molecular
systems. The molecular integrals are obtained by PySCF,[28] and the modified
Cholesky decomposition is performed with ipie, with a threshold of $10^{-5}$.
We used a time step of 0.01 a.u., and periodic reorthogonalization of walkers
is performed every 5 steps for all calculations. We used the restricted
Hartree Fock (RHF) trial for all systems. For accurate statistical analysis,
we perform block analysis for $\geq 200$ gradient samples for each
calculation. All AFQMC calculations are performed using the frozen-core
approximation.

We benchmarked our implementation on various small molecules using aug-cc-pVTZ
[64] basis set and compared the results to Ref. 33 in Table 1. All dipole
moments align strongly with the AD-AFQMC results reported in Ref. 33.
Furthermore, except for CO, AD-AFQMC matches experimental values better than
RCCSD and exhibits accuracy comparable to RCCSD(T). This supports the widely
held view that AFQMC’s accuracy falls between CCSD and CCSD(T).

Table 1:  Comparison of AD-AFQMC dipole moment (in a.u.) of various molecules at equilibrium geometry with the implementation in Ref. 33 and other quantum chemistry methods. The data using RCCSD and RCCSD(T) are also extracted from Ref. 33. Molecule | This work | AD-AFQMC in Ref. 33 | RCCSD | RCCSD(T) | Experiment
---|---|---|---|---|---
H2O | 0.723(2) | 0.720(2) | 0.7335 | 0.7247 | 0.730[65]
NH3 | 0.592(2) | 0.592(2) | 0.6015 | 0.5938 | 0.581(1)[66]
CO | 0.022(4) | 0.019(4) | 0.0199 | 0.0429 | 0.048(1)[67]
HCl | 0.428(1) | 0.429(1) | 0.4318 | 0.4273 | 0.430[68]
HBr | 0.332(2) | 0.329(2) | 0.3289 | 0.3245 | 0.325[68]

###  IV.6 Enhanced integration testing and no-MPI mode

As ipie expands its functionalities, robust testing workflows become crucial.
ipie supports a new integration testing framework that enhances the package’s
robustness and adaptability and significantly contributes to its reliability
and ease of use for the end-users. The GitHub continuous integration (CI)
workflow automatically tests new pull requests, ensuring that every change to
the codebase does not break existing functionalities. The workflow encompasses
a comprehensive test suite and executes various linting and code formatting
checks before running serial and parallel unit tests and integration tests.

Recognizing the diverse computational environments in which ipie might be
deployed, we introduced the no-MPI mode. This mode is designed for situations
where MPI is unavailable, or its use is not desired, offering greater
flexibility for users and developers. The CI workflow includes a job that
tests ipie’s functionality without the MPI dependency. This mode is
particularly beneficial for users who wish to perform quantum Monte Carlo
simulations on personal computing setups or in environments where setting up
MPI is challenging.

##  V Interfaces to external packages

###  V.1 Dice interface and SHCI-AFQMC

ipie includes utilities for converting the output from Dice, [19] a package
that employs semistochastic heat bath configuration interaction (SHCI) [29,
19] as the complete active space (CAS) solver.

python -u /path-to-ipie/tools/extract_dice.py –dice-wfn /path-to-Dice-
output/dets.bin –sort –verbose

Figure 23: Convert Dice output to the MSD trial in ipie.

which produces the wfn.h5 file containing the coefficients and the indices of
occupied orthogonal orbitals for spin-up and spin-down sectors that follows
the block-formatted ($\alpha\alpha\cdots\beta\beta\cdots$) orbital convention
in ipie. The SHCI calculation generates a good MSD trial for challenging
systems. [15, 39] However, it is not a black-box approach and often requires
careful handling, such as selecting CAS and using natural orbitals, as
discussed in Section IV.3.

These detailed discussions are beyond the scope of this article, and we
briefly mention the procedures in the semi-blackbox example provided within
ipie: [15]

  1. 1.

Initiate a preliminary rough SHCI calculation in an extensive active space for
the system under study.

  2. 2.

Derive the SHCI one-electron reduced density matrix (1-RDM), extracting the
resultant natural orbital occupation number (NOON) and natural orbitals.

  3. 3.

Define an active space criterion based on a predetermined NOON threshold.

  4. 4.

Adjust the orbitals through rotation, aligning them with the unitary
transformation specified by the natural orbitals.

  5. 5.

Execute a refined SHCI self-consistent field calculation within the new active
space determined in the previous step.

  6. 6.

Use this MSD trial in AFQMC.

This strategy ensures the trial wavefunction encapsulates static correlations
within the active space via the MSD trial. At the same time, AFQMC
incorporates the residual dynamic correlations. This procedure is folded into
the factory utility method ipie.utils.from_dice.build_driver_from_shciscf.

###  V.2 TREXIO support and CIPSI-AFQMC

The TrexIO library and file format have been developed to offer a robust and
efficient solution for storing and exchanging wavefunction parameters and
matrix elements.[30] This library supports bindings in several programming
languages, including Python, and can be conveniently installed via the pip
package manager.

The compatibility of ipie with the TrexIO format facilitates its integration
with various software packages. Specifically, it allows ipie to utilize trial
wavefunctions produced by Quantum Package,[69] along with the associated one-
electron integrals and Cholesky-decomposed electron repulsion integrals. This
interface allowed us to check if AFQMC could complement configuration
interaction using perturbative selection done iteratively (CIPSI) calculations
to improve full configuration interaction (FCI) energy estimates of large
systems.

Typically, to estimate the FCI energy from a CIPSI calculation, one
extrapolates to zero the variational energy, $E_{\text{trial}}$, as a function
of the renormalized second-order perturbative correction, $E_{\text{rPT2}}$.
[69] However, the AFQMC energy, $E_{\text{AFQMC}}$, is anticipated to provide
a closer approximation to the FCI energy than the sum of $E_{\text{trial}}$
and $E_{\text{rPT2}}$. This expectation is particularly relevant for systems
with large $E_{\text{rPT2}}$ corrections.

Figure 24: Energy of the trial wavefunction ($E_{\text{trial}}$) as a function of $\Delta E$, where $\Delta E=E_{\text{rPT2}}$ for CIPSI calculations and $\Delta E=E_{\text{AFQMC}}-E_{\text{trial}}$ for AFQMC. exFCI is the extrapolated FCI energy obtained from CIPSI calculations. Table 2:  Renormalized PT2 correction $E_{\text{rPT2}}$, Variational energy $E_{\text{var}}$, and AFQMC energy $E_{\text{AFQMC}}$ for the nitroxyl and benzene molecules as functions of the number of determinants in the variational space. System, Basis | $N_{\text{det}}$ | $E_{\text{trial}}$ | $E_{\text{trial}}+E_{\text{rPT2}}$ | $E_{\text{AFQMC}}$
---|---|---|---|---
HNO, 6-31G | $1$ | $-129.712554$ | $-129.987361$ | $-129.971(1)$
| $2$ | $-129.752113$ | $-129.987158$ | $-129.974(1)$
| $41$ | $-129.781149$ | $-129.981219$ | $-129.9727(6)$
| $748$ | $-129.894775$ | $-129.973405$ | $-129.9724(2)$
| $2838$ | $-129.942391$ | $-129.975057$ | $-129.9726(1)$
| $14201$ | $-129.955368$ | $-129.973620$ | $-129.9727(1)$
C6H6, cc-pVDZ | $1$ | $-230.7204904$ | $-231.39330$ | $-231.5887(7)$
| $23$ | $-230.7615816$ | $-231.42167$ | $-231.5866(7)$
| $828$ | $-230.8687072$ | $-231.44602$ | $-231.5864(6)$
| $15690$ | $-231.1300087$ | $-231.49330$ | $-231.5848(4)$
| $109869$ | $-231.3152399$ | $-231.53022$ | $-231.5842(4)$

We conducted AFQMC calculations on the nitroxyl and benzene molecules using
CIPSI trial wavefunctions of increasing sizes. The results of these
calculations are detailed in Table 2. Supporting our expectation, Fig. 24
demonstrates that the AFQMC corrections align the data points along a straight
line, [70] validating the hypothesis that AFQMC energies are more reliable for
the extrapolation towards the FCI value, especially in cases with large rPT2
corrections. Employing calculations with comparatively small wavefunctions for
benzene, as shown in Table 2, a three-point linear extrapolation based on the
rPT2 correction yields a correlation energy of $-858.6\text{\,}\mathrm{m\text{
$E_{\textup{ h }}$ }}$. In contrast, extrapolation using AFQMC energies
results in a correlation energy of $-862.2\text{\,}\mathrm{m\text{
$E_{\textup{ h }}$ }}$. This latter value is significantly closer to the
correlation energy of $-863.4\text{\,}\mathrm{m\text{ $E_{\textup{ h }}$ }}$
achieved through CIPSI with 167 million determinants.[71]

These preliminary calculations enabled by the TrexIO interface illustrate that
the integration of AFQMC with CIPSI emerges as a promising methodology for
estimating the FCI energy of systems larger than those currently feasible.

###  V.3 FQE interfaces

The Fermionic quantum emulator (FQE) [31, *rubin2021fermionic] is a
lightweight fermionic circuit simulator, which is particularly useful in
quantum computing where it aids in the development and testing of quantum
algorithms tailored for fermionic systems. ipie provides the conversion
between the ipie’s MSD wavefunction and the FQE wavefunction.

##  VI Conclusions and outlooks

This paper summarized the improvements and new features added in ipie since
its original release.[15] These improvements enhance modularity and
computational efficiency and offer intuitive user-end APIs. New features and
interfaces aim to expand a broader spectrum of AFQMC calculations in quantum
chemistry.

We summarize the key features we highlighted in this manuscript:

  1. 1.

Distributed Hamiltonians to remove the memory bottleneck. We demonstrated
ipie’s capacity for studying large systems deploying GPUs with significantly
higher efficiency than CPU-based implementations, exemplified in our case
study assessing the interaction energies in a benzene dimer.

  2. 2.

GPU support for MSD trial wavefunctions. With customized CUDA kernels, we
enabled an efficient realization of Wick’s theorem. Timing benchmarks for
[Cu2O2]2+ and [Fe2S2(SCH3)]2- were shown to achieve more than an order of
magnitude speedup compared to our CPU implementation for a large MSD trial.
[49]

  3. 3.

Support for complex-valued Cholesky vectors. ipie can handle complex-valued
Cholesky vectors that may arise when the underlying basis functions are
complex-valued.

  4. 4.

Free-projection AFQMC. A numerically exact AFQMC approach can be used to study
small strongly correlated systems.

  5. 5.

Finite-temperature AFQMC. A finite-temperature AFQMC algorithm based on the
grand canonical ensemble was added.

  6. 6.

Electron-phonon QMC. A QMC algorithm that computes the ground state of
electron-phonon problems was added.

  7. 7.

Automatic differentiable AFQMC. We offer AFQMC property calculations via
automatic differentiation.

  8. 8.

External package interfaces. ipie is now interfaced with PySCF, Dice, TrexIO,
and FQE.

We hope that ipie will serve as a community code base for developing ab initio
AFQMC methods and their applications. Furthermore, as ipie is written mainly
in Python, we anticipate its use in machine learning and quantum computing
communities will also grow.

##  VII Acknowledgements

Work by T.J. and J.L. were supported by Harvard University’s startup funds,
the quantum information science project from the U.S. Department of Energy
(DOE), and Wellcome Leap as part of the Quantum for Bio Program. Computations
were carried out partly on the FASRC cluster supported by the FAS Division of
Science Research Computing Group at Harvard University. This work also used
the Delta system at the National Center for Supercomputing Applications
through allocation CHE230032, CHE230088, and PHY230192 from the Advanced
Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS)
program, which is supported by National Science Foundation grants #2138259,
#2138286, #2138307, #2137603, and #2138296. P.F.L. and A.S. have received
financial support from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (Grant agreement no.
863481). A.S. was also supported by the European Centre of Excellence in
Exascale Computing (TREX) which has received funding from the European Union’s
Horizon 2020 — Research and Innovation program — under grant agreement no.
952165. S.F.U. thanks David Reichman for his support. We thank Yifei Huang and
Dingshun Lv for their GPU-MSD code contributions. [49]

##  VIII Data Availability

The data that support the findings of this study are openly available in our
Zenodo repository at https://doi.org/10.5281/zenodo.12522916.

## References

  * Zhang, Carlson, and Gubernatis [1995] S. Zhang, J. Carlson, and J. E. Gubernatis, “Constrained Path Quantum Monte Carlo Method for Fermion Ground States,” Phys. Rev. Lett. 74, 3652–3655 (1995).
  * Zhang and Krakauer [2003] S. Zhang and H. Krakauer, “Quantum Monte Carlo Method using Phase-Free Random Walks with Slater Determinants,” Phys. Rev. Lett. 90, 136401 (2003).
  * Lee, Pham, and Reichman [2022] J. Lee, H. Q. Pham, and D. R. Reichman, “Twenty Years of Auxiliary-Field Quantum Monte Carlo in Quantum Chemistry: An Overview and Assessment on Main Group Chemistry and Bond-Breaking,” J. Chem. Theory Comput. 18, 7024–7042 (2022).
  * Motta and Zhang [2018] M. Motta and S. Zhang, “Ab initio computations of molecular systems by the auxiliary-field quantum Monte Carlo method,” WIREs Computational Molecular Science 8, e1364 (2018).
  * Huggins  _et al._ [2022] W. J. Huggins, B. A. O’Gorman, N. C. Rubin, D. R. Reichman, R. Babbush, and J. Lee, “Unbiasing fermionic quantum Monte Carlo with a quantum computer,” Nature 603, 416–420 (2022).
  * Wan  _et al._ [2023] K. Wan, W. J. Huggins, J. Lee, and R. Babbush, “Matchgate shadows for fermionic quantum simulation,” Commun. Math. Phys 404, 629–700 (2023).
  * Amsler  _et al._ [2023] M. Amsler, P. Deglmann, M. Degroote, M. P. Kaicher, M. Kiser, M. Kühn, C. Kumar, A. Maier, G. Samsonidze, A. Schroeder, _et al._ , “Classical and quantum trial wave functions in auxiliary-field quantum monte carlo applied to oxygen allotropes and a cubr2 model system,” J. Chem. Phys. 159 (2023).
  * Kiser  _et al._ [2023] M. Kiser, A. Schroeder, G.-L. R. Anselmetti, C. Kumar, N. Moll, M. Streif, and D. Vodola, “Classical and quantum cost of measurement strategies for quantum-enhanced auxiliary field quantum monte carlo,” New J. Phys. (2023).
  * Huang  _et al._ [2024a] B. Huang, Y.-T. Chen, B. Gupt, M. Suchara, A. Tran, S. McArdle, and G. Galli, “Evaluating a quantum-classical quantum Monte Carlo algorithm with Matchgate shadows,” arXiv (2024a), 2404.18303 .
  * Lee, Malone, and Morales [2019] J. Lee, F. D. Malone, and M. A. Morales, “An auxiliary-field quantum monte carlo perspective on the ground state of the dense uniform electron gas: An investigation with hartree-fock trial wavefunctions,” J. Chem. Phys. 151 (2019).
  * Lee and Reichman [2020] J. Lee and D. R. Reichman, “Stochastic resolution-of-the-identity auxiliary-field quantum monte carlo: Scaling reduction without overhead,” J. Chem. Phys. 153 (2020), 10.1063/5.0015077.
  * Lee, Morales, and Malone [2021] J. Lee, M. A. Morales, and F. D. Malone, “A phaseless auxiliary-field quantum monte carlo perspective on the uniform electron gas at finite temperatures: Issues, observations, and benchmark study,” J. Chem. Phys. 154 (2021).
  * Lee, Zhang, and Reichman [2021] J. Lee, S. Zhang, and D. R. Reichman, “Constrained-path auxiliary-field quantum monte carlo for coupled electrons and phonons,” Phys. Rev. B 103, 115123 (2021).
  * Lee  _et al._ [2021] J. Lee, F. D. Malone, M. A. Morales, and D. R. Reichman, “Spectral functions from auxiliary-field quantum monte carlo without analytic continuation: The extended koopmans’ theorem approach,” J. Chem. Theory Comput. 17, 3372–3387 (2021).
  * Malone  _et al._ [2022] F. D. Malone, A. Mahajan, J. S. Spencer, and J. Lee, “Ipie: A python-based auxiliary-field quantum monte carlo program with flexibility and efficiency on cpus and gpus,” J. Chem. Theory Comput. 19, 109–121 (2022).
  * Lam, Pitrou, and Seibert [2015] S. K. Lam, A. Pitrou, and S. Seibert, “Proceedings of the second workshop on the llvm compiler infrastructure in hpc. llvm’15,” (2015).
  * Cramer  _et al._ [2006] C. J. Cramer, M. Włoch, P. Piecuch, C. Puzzarini, and L. Gagliardi, “Theoretical models on the cu2o2 torture track: Mechanistic implications for oxytyrosinase and small-molecule analogues,” J. Phys. Chem. A 110, 1991–2004 (2006).
  * Kent  _et al._ [2020] P. R. 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ää, _et al._ , “Qmcpack: Advances in the development, efficiency, and application of auxiliary field and real-space variational and diffusion quantum monte carlo,” J. Chem. Phys. 152 (2020).
  * Sharma  _et al._ [2017] S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi, and C. J. Umrigar, “Semistochastic heat-bath configuration interaction method: Selected configuration interaction with semistochastic perturbation theory,” J. Chem. Theory Comput. 13, 1595–1604 (2017).
  * Hehn, Deglmann, and Kühn [2024] L. Hehn, P. Deglmann, and M. Kühn, “Chelate complexes of 3d transition metal ions–a challenge for electronic-structure methods?” J. Chem. Theory Comput. (2024).
  * Vysotskiy, Filippi, and Ryde [2024] V. P. Vysotskiy, C. Filippi, and U. Ryde, “Scalar relativistic all-electron and pseudopotential ab initio study of a minimal nitrogenase [fe (sh) 4h]- model employing coupled-cluster and auxiliary-field quantum monte carlo many-body methods,” J. Phys. Chem. A (2024).
  * Chen  _et al._ [2023a] M. S. Chen, J. Lee, H.-Z. Ye, T. C. Berkelbach, D. R. Reichman, and T. E. Markland, “Data-efficient machine learning potentials from transfer learning of periodic correlated electronic structure methods: Liquid water at afqmc, ccsd, and ccsd (t) accuracy,” J. Chem. Theory Comput. 19, 4510–4519 (2023a).
  * Jiang  _et al._ [2024] T. Jiang, B. O’Gorman, A. Mahajan, and J. Lee, “Unbiasing fermionic auxiliary-field quantum monte carlo with matrix product state trial wavefunctions,” arXiv (2024), 2405.05440 .
  * Mahajan and Sharma [2021] A. Mahajan and S. Sharma, “Taming the sign problem in auxiliary-field quantum monte carlo using accurate wave functions,” J. Chem. Theory Comput. 17, 4786–4798 (2021).
  * Vitali, Rosenberg, and Zhang [2019] E. Vitali, P. Rosenberg, and S. Zhang, “Calculating ground state properties of correlated fermionic systems with BCS trial wave functions in Slater determinant path-integral approaches,” Phys. Rev. A 100, 023621 (2019).
  * Chang, Rubenstein, and Morales [2016] C.-C. Chang, B. M. Rubenstein, and M. A. Morales, “Auxiliary-field-based trial wave functions in quantum Monte Carlo calculations,” Phys. Rev. B 94, 235144 (2016).
  * [27] “Plum: Multiple dispatch in python,” https://github.com/beartype/plum .
  * 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, _et al._ , “Pyscf: the python-based simulations of chemistry framework,” Wiley Interdiscip Rev Comput Mol Sci 8, e1340 (2018).
  * Holmes, Tubman, and Umrigar [2016] A. A. Holmes, N. M. Tubman, and C. Umrigar, “Heat-bath configuration interaction: An efficient selected configuration interaction algorithm inspired by heat-bath sampling,” J. Chem. Theory Comput. 12, 3674–3680 (2016).
  * Posenitskiy  _et al._ [2023] E. Posenitskiy, V. G. Chilkuri, A. Ammar, M. Hapka, K. Pernal, R. Shinde, E. J. Landinez Borda, C. Filippi, K. Nakano, O. Kohulák, S. Sorella, P. de Oliveira Castro, W. Jalby, P. L. Ríos, A. Alavi, and A. Scemama, “TREXIO: A file format and library for quantum chemistry,” J. Chem. Phys. 158 (2023), 10.1063/5.0148161, 174801.
  * Rubin  _et al._ [2021a] N. C. Rubin, K. Gunst, A. White, L. Freitag, K. Throssell, G. K.-L. Chan, R. Babbush, and T. Shiozaki, “The fermionic quantum emulator, version 0.2.0,”  (2021a).
  * Rubin  _et al._ [2021b] N. C. Rubin, K. Gunst, A. White, L. Freitag, K. Throssell, G. K.-L. Chan, R. Babbush, and T. Shiozaki, “The fermionic quantum emulator,” Quantum 5, 568 (2021b).
  * Mahajan  _et al._ [2023] A. Mahajan, J. S. Kurian, J. Lee, D. R. Reichman, and S. Sharma, “Response properties in phaseless auxiliary field quantum Monte Carlo,” J. Chem. Phys. 159, 184101 (2023).
  * Hubbard [1959] J. Hubbard, “Calculation of partition functions,” Phys. Rev. Lett. 3, 77 (1959).
  * Stratonovich [1958] R. Stratonovich, Sov. Phys. Dokl. 2, 416 (1958).
  * Ren  _et al._ [2021] J. Ren, W. Li, T. Jiang, Y. Wang, and Z. Shuai, “The renormalizer package,” https://github.com/shuaigroup/renormalizer (2021).
  * Ren  _et al._ [2022] J. Ren, W. Li, T. Jiang, Y. Wang, and Z. Shuai, “Time-dependent density matrix renormalization group method for quantum dynamics in complex systems,” Wiley Interdiscip Rev Comput Mol Sci 12, e1614 (2022).
  * Flyvbjerg and Petersen [1989] H. Flyvbjerg and H. G. Petersen, “Error estimates on averages of correlated data,” J. Chem. Phys. 91, 461–466 (1989).
  * Mahajan, Lee, and Sharma [2022] A. Mahajan, J. Lee, and S. Sharma, “Selected configuration interaction wave functions in phaseless auxiliary field quantum monte carlo,” J. Chem. Phys. 156 (2022), https://doi.org/10.1063/5.0087047.
  * Jurečka  _et al._ [2006] P. Jurečka, J. Šponer, J. Černỳ, and P. Hobza, “Benchmark database of accurate (mp2 and ccsd (t) complete basis set limit) interaction energies of small model complexes, dna base pairs, and amino acid pairs,” Phys. Chem. Chem. Phys. 8, 1985–1993 (2006).
  * Weigend, Kattannek, and Ahlrichs [2009] F. Weigend, M. Kattannek, and R. Ahlrichs, “Approximated electron repulsion integrals: Cholesky decomposition versus resolution of the identity methods,” J. Chem. Phys. 130 (2009).
  * Halkier  _et al._ [1998] A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, H. Koch, J. Olsen, and A. K. Wilson, “Basis-set convergence in correlated calculations on ne, n2, and h2o,” Chem. Phys. Lett. 286, 243–252 (1998).
  * gat [2024] “Benchmark data for S22 reaction set,”  (2024).
  * Burns, Marshall, and Sherrill [2014] L. A. Burns, M. S. Marshall, and C. D. Sherrill, “Appointing silver and bronze standards for noncovalent interactions: A comparison of spin-component-scaled (scs), explicitly correlated (f12), and specialized wavefunction approaches,” J. Chem. Phys. 141 (2014), https://doi.org/10.1063/1.4903765.
  * Al-Hamdani  _et al._ [2021] Y. S. Al-Hamdani, P. R. Nagy, A. Zen, D. Barton, M. Kállay, J. G. Brandenburg, and A. Tkatchenko, “Interactions between large molecules pose a puzzle for reference quantum mechanical methods,” Nat. Commun. 12, 3927 (2021).
  * Boys and Bernardi [1970] S. F. Boys and F. Bernardi, “The calculation of small molecular interactions by the differences of separate total energies. some procedures with reduced errors,” Mol. Phys. 19, 553–566 (1970).
  * Blaschke and Stopkowicz [2022] S. Blaschke and S. Stopkowicz, “Cholesky decomposition of complex two-electron integrals over giaos: Efficient mp2 computations for large molecules in strong magnetic fields,” J. Chem. Phys. 156 (2022), https://doi.org/10.1063/5.0076588.
  * Suewattana  _et al._ [2007] M. Suewattana, W. Purwanto, S. Zhang, H. Krakauer, and E. J. Walter, “Phaseless auxiliary-field quantum monte carlo calculations with plane waves and pseudopotentials: Applications to atoms and molecules,” Phys. Rev. B 75, 245123 (2007).
  * Huang  _et al._ [2024b] Y. Huang, Z. Guo, H. Q. Pham, and D. Lv, “Gpu-accelerated auxiliary-field quantum monte carlo with multi-slater determinant trial states,” arXiv (2024b), 2406.08314 .
  * Li and Chan [2017] Z. Li and G. K.-L. Chan, “Spin-projected matrix product states: Versatile tool for strongly correlated systems,” J. Chem. Theory Comput. 13, 2681–2695 (2017).
  * Assaraf and Caffarel [1999] R. Assaraf and M. Caffarel, “Zero-Variance Principle for Monte Carlo Algorithms,” Phys. Rev. Lett. 83, 4682–4685 (1999).
  * Chen  _et al._ [2023b] Y. Chen, L. Zhang, W. E, and R. Car, “Hybrid auxiliary field quantum monte carlo for molecular systems,” J. Chem. Theory Comput. 19, 4484–4493 (2023b).
  * He  _et al._ [2019] Y.-Y. He, M. Qin, H. Shi, Z.-Y. Lu, and S. Zhang, “Finite-temperature auxiliary-field quantum monte carlo: Self-consistent constraint and systematic approach to low temperatures,” Phys. Rev. B 99, 045108 (2019).
  * Shen  _et al._ [2020] T. Shen, Y. Liu, Y. Yu, and B. M. Rubenstein, “Finite temperature auxiliary field quantum Monte Carlo in the canonical ensemble,” J. Chem. Phys. 153 (2020), 10.1063/5.0026606.
  * Blankenbecler, Scalapino, and Sugar [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, “Monte carlo calculations of coupled boson-fermion systems. i,” Phys. Rev. D 24, 2278–2286 (1981).
  * Hirsch [1985] J. E. Hirsch, “Two-dimensional hubbard model: Numerical simulation study,” Phys. Rev. B 31, 4403–4419 (1985).
  * Santos [2003] R. R. d. Santos, “Introduction to quantum monte carlo simulations for fermionic systems,” Braz. J. Phys. 33, 36–54 (2003).
  * Macridin [2003] A. Macridin, “Phonons, Charge and Spin in Correlated Systems,” PhD Thesis (2003).
  * Marchand, Stamp, and Berciu [2017] D. J. J. Marchand, P. C. E. Stamp, and M. Berciu, “Dual coupling effective band model for polarons,” Phys. Rev. B 95, 035117 (2017).
  * Kalos, Levesque, and Verlet [1974] M. H. Kalos, D. Levesque, and L. Verlet, “Helium at zero temperature with hard-sphere and other forces,” Phys. Rev. A 9, 2178–2195 (1974).
  * Holstein [1959] T. Holstein, “Studies of polaron motion : Part I. The molecular-crystal model,” Ann. Phys. 8, 325–342 (1959).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” in  _Advances in Neural Information Processing Systems 32_ (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Chen  _et al._ [2016] T. Chen, B. Xu, C. Zhang, and C. Guestrin, “Training deep nets with sublinear memory cost,” arXiv:1604.06174 (2016), https://doi.org/10.48550/arXiv.1604.06174.
  * Dunning Jr [1989] T. H. Dunning Jr, “Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen,” J. Chem. Phys. 90, 1007–1023 (1989).
  * Shostak, Ebenstein, and Muenter [1991] S. L. Shostak, W. L. Ebenstein, and J. S. Muenter, “The dipole moment of water. i. dipole moments and hyperfine properties of h2o and hdo in the ground and excited vibrational states,” J. Chem. Phys. 94, 5875–5882 (1991).
  * Shimizu [1970] F. Shimizu, “Stark spectroscopy of nh3 $\nu$ 2 band by 10-$\mu$ co2 and n2o lasers,” J. Chem. Phys. 52, 3572–3576 (1970).
  * Muenter [1975] J. Muenter, “Electric dipole moment of carbon monoxide,” J. Mol. Spectrosc. 55, 490–491 (1975).
  * Lovas  _et al._ [2003] F. Lovas, E. Tiemann, J. Coursey, S. Kotochigova, J. Chang, K. Olsen, and R. Dragoset, “Diatomic spectral database,” (2003).
  * Garniron  _et al._ [2019] Y. Garniron, T. Applencourt, K. Gasperich, A. Benali, A. Ferté, J. Paquier, B. Pradines, R. Assaraf, P. Reinhardt, J. Toulouse, P. Barbaresco, N. Renon, G. David, J.-P. Malrieu, M. Véril, M. Caffarel, P.-F. Loos, E. Giner, and A. Scemama, “Quantum package 2.0: An open-source determinant-driven suite of programs,” J. Chem. Theory Comput. 15, 3591–3609 (2019).
  * Burton and Loos [2024] H. G. A. Burton and P.-F. Loos, “Rationale for the extrapolation procedure in selected configuration interaction,” J. Chem. Phys. 160, 104102 (2024).
  * Loos, Damour, and Scemama [2020] P.-F. Loos, Y. Damour, and A. Scemama, “The performance of CIPSI on the ground state electronic energy of benzene,” J. Chem. Phys. 153 (2020), 10.1063/5.0027617.