# Selected non-orthogonal configuration interaction with compressed single and
double excitations

Chong Sun  Department of Chemistry, Rice University, Houston, TX 77005-1892
Fei Gao  Department of Physics and Astronomy, Rice University, Houston, TX
77005-1892  Gustavo E. Scuseria  Department of Chemistry, Rice University,
Houston, TX 77005-1892  Department of Physics and Astronomy, Rice University,
Houston, TX 77005-1892<EMAIL_ADDRESS>
###### Abstract

Addressing both dynamic and static correlation accurately is a primary goal in
electronic structure theory. Non-orthogonal configuration interaction (NOCI)
is a versatile tool for treating static correlation, offering chemical
insights by combining diverse reference states. Nevertheless, achieving
quantitative accuracy requires the inclusion of missing dynamic correlation.
This work introduces a framework for compressing orthogonal single and double
excitations into an NOCI of a much smaller dimension. This compression is
repeated with each Slater determinant in a reference NOCI, resulting in
another NOCI that includes all its single and double excitations (NOCISD),
effectively recovering the missing dynamic correlations from the reference.
This compressed NOCISD is further refined through a selection process using
metric and energy tests (SNOCISD). We validate the effectiveness of SNOCISD
through its application to the dissociation of the nitrogen molecule and the
hole-doped two-dimensional Hubbard model at various interaction strengths.

††preprint: AIP/123-QED

##  I Introduction

Modern quantum chemistry is increasingly focused on accurate simulations of
strongly correlated systems. These systems, where Hartree-Fock falls short,
include challenges like excited states, nonequilibrium geometries, and many
transition-metal compounds, among others. Post-Hartree-Fock (post-HF) methods,
such as Møller–Plesset perturbation theory (MP), truncated configuration
interaction (CI), and coupled cluster theory (CC), provide improvement over
the HF reference by incorporating dynamic (or weak) correlations from
particle-hole excitations. However, these methods encounter significant
challenges when strong (static or non-dynamic) correlations become
significant. Addressing static correlation often involves considering multiple
reference configurations. The Complete Active Space Self-Consistent Field
(CASSCF) method, Roos (1980) incorporating all possible configurations within
a chosen active space, was developed to tackle static correlation. However, as
the active space size increases, exact solvers such as full configuration
interaction (FCI) become intractable, hindering large-scale simulations.

On the other hand, non-orthogonal configuration interaction (NOCI)Thom and
Head-Gordon (2009); Jiménez-Hoyos, Rodríguez-Guzmán, and Scuseria (2013);
Jimenez-Hoyos (2013) offers a compact and versatile alternative to CASSCF
without sacrificing accuracy. In NOCI, each electronic configuration, defined
by its unique set of molecular orbitals (MOs), is not necessarily orthogonal
to others. This "different orbitals for different configurations" approach
endows NOCI with great flexibility, requiring fewer configurations to describe
sets of orthogonal excitations. Its diabatic nature makes NOCI particularly
useful in chemically intuitive studies of phenomena like avoided crossings and
conical intersections. Barca, Gilbert, and Gill (2014); Jensen  _et al._
(2018); Mahler and Thompson (2021) Methods that can generate an NOCI
wavefunction include symmetry breaking and restoration in projected Hartree-
Fock (PHF) theory, Smeyers and Doreste-Suarez (1973); Smeyers and Delgado-
Barrio (1974); Jiménez-Hoyos  _et al._ (2012); Garza, Jiménez-Hoyos, and
Scuseria (2013); Tsuchimochi and Ten-no (2016) or combining non-unitary
Thouless rotations Thouless (1960) over a reference Slater determinant
(henceforth, determinant). The Thouless parameters defining the latter
strategy can be optimized with either the resonating HF (ResHF) approach
Fukutome (1988) or the few-determinant (FED) method. Schmid  _et al._ (1989);
Schmid (2004) ResHF optimizes all determinants simultaneously, whereas FED
only optimizes determinants added to the expansion, keeping the previous ones
frozen. ResHF often yields better energies, but it requires a more involved
optimization process not always worth the effort. Jiménez-Hoyos, Rodríguez-
Guzmán, and Scuseria (2013) As more determinants are included in the NOCI
wavefunction, FCI accuracy is asymptotically approached both with ResHF and
FED. However, the energy improvement usually reaches a plateau, Jiménez-Hoyos,
Rodríguez-Guzmán, and Scuseria (2013) highlighting a saturation point where
dynamic correlation is added in a less-than-efficient manner. This observation
underscores the challenge we want to address in this work: how can we
effectively incorporate dynamic correlation into an NOCI whose static
correlation is predominantly included by the initial terms of the expansion?
We here propose an innovative strategy.

Dynamic correlation missing from a given multireference wavefunction can be
effectively recovered by including low-rank excitations. However, applying
these concepts to NOCI wavefunctions, which use different MOs for each
configuration, is not straightforward. Efforts to tackle this NOCI challenge
include Ten-no’s addition of CI and CC, Ten-no (1998) Burton and Thom’s
second-order perturbation theory, Burton and Thom (2020) and Nite and Jimenez-
Hoyos’ implementation of the CISD expansion for each reference configuration.
Nite and Jiménez-Hoyos (2019) The latter necessitates a significant truncation
to manage computational demands.

This paper introduces a framework for compressing non-orthogonal configuration
interaction with single and double-electron excitations (NOCISD) from each
term, without brute-force truncation. The key component of our method is a
strategy for compressing linear combinations of particle-hole excitations from
each determinant into a significantly smaller set of non-orthogonal Slater
determinants (NOSD). The method is applicable across all excitation orders.
Its efficiency is demonstrated by showing that in the spin-orbital basis, only
$2$ NOSD are required for reproducing an expansion with all singly-excited
states (CIS), and at most $(2N_{\mathrm{vir}}N_{\mathrm{occ}}+1)$ NOSD are
needed for all doubly-excited states (CID), achieving a square-root reduction
in expansion length. After compressing the excited state set, we employ a
selection process based on metric and energy criteria to further reduce the
NOCI expansion without compromising accuracy. The resulting method is named
selected NOCISD (SNOCISD); it provides a versatile and compact approach for
incorporating residual dynamic correlation into an NOCI wavefunction. The
performance of SNOCISD is here benchmarked on the N2 molecule dissociation
curve and the ground state energy of the hole-doped two-dimensional Hubbard
model at different interaction regimes. In addition, we analyze the SNOCISD
spin-spin correlation function of the two-dimensional case, providing insights
into its structure.

##  II Theory

In this section, we outline the SNOCISD algorithm, which begins by generating
a small set of non-orthogonal reference determinants through ResHF or FED
optimization. Subsequently, a CISD procedure is applied to each reference
determinant, followed by their compression. The following step involves
selecting the most important determinants from these compressed sets, which
are then incorporated into an expanded NOCI. Finally, a generalized eigenvalue
problem is solved with respect to the expanded NOCI. Detailed formulae and
derivations describing the key steps necessary for implementing the SNOCISD
algorithm are presented below.

###  II.1 Non-orthogonal configuration interaction

Mutually non-orthogonal Slater determinants (NOSD) can be constructed through
Thouless rotations applied to a reference determinant $|\Phi_{0}\rangle$. A
new determinant is defined by

$\begin{split}|\Phi(\mathbf{Z})\rangle=e^{Z}|\Phi_{0}\rangle,\quad
Z=\sum_{a=1}^{N_{\rm vir}}\sum_{i=1}^{N_{\rm occ}}Z_{ai}E_{a}^{i},\end{split}$
(1)

where $Z$ represents the Thouless rotation operator, characterized by
$N_{\mathrm{vir}}N_{\mathrm{occ}}$ parameters $\mathbf{Z}$. The excitation
operator $E_{a}^{i}=b^{\dagger}_{a}b_{i}$ transfers an electron from the $i$th
occupied orbital to the $a$th virtual orbital. We use the subscripts "occ" and
"vir" to represent occupied and virtual spaces, respectively. Assuming
$|\Phi_{0}\rangle$ is normalized, $|\Phi(\mathbf{Z})\rangle$ maintains
intermediate normalization, meaning
$\langle\Phi(\mathbf{Z})|\Phi_{0}\rangle=1$. In the limiting case where one
parameter $Z_{ai}\gg 1$ and the others are zero,
$|\Phi(\mathbf{Z})\rangle=|\Phi_{0}\rangle+Z_{ai}E_{i}^{a}|\Phi_{0}\rangle$
becomes a singly excited state from $|\Phi_{0}\rangle$. Higher-order
excitations can be similarly obtained when two or more parameters are large
and the rest are zero. This is the construction yielding the traditional FCI
using a linear combination of NOSD.

In our case, we build the many-body wavefunction as

$\begin{split}|\Psi\rangle=\sum_{\mu=0}^{L-1}w_{\mu}|\Phi_{\mu}\rangle,\end{split}$
(2)

where $\mathbf{w}=[w_{0},\cdots,w_{L-1}]$ are the coefficients of each
$|\Phi_{\mu}\rangle=|\Phi(\mathbf{Z_{\mu}})\rangle$. Typically,
$|\Phi_{0}\rangle$ is the HF solution. Each $|\Phi_{\mu}\rangle$ is based on a
distinct set of molecular orbitals (MOs) differing from those derived from HF.

The coefficients $\mathbf{w}$ are obtained by solving a generalized eigenvalue
equation, using the generalized Wick’s theoremLöwdin (1955); Ripka, Blaizot,
and Ripka (1986); Balian and Brezin (1969); Chen and Scuseria (2023) to
compute the matrix elements. The NOCI energy is minimized by optimizing the
Thouless parameters $\\{\mathbf{Z}_{\mu}\\}$ with the aforementioned ResHF or
FED methods.

###  II.2 Compression of excited configurations

We next discuss how a small number of NOSD can represent a linear combination
of excited configurations like single and double excitations. We work with
spin-orbital excitations, while applications to more specific spin cases are
straightforward. We start by adding a real parameter $t$ to the Thouless
rotation in Eq. (1),

$\begin{split}|\Phi(\mathbf{Z},t)\rangle=e^{tZ}|\Phi_{0}\rangle,\end{split}$
(3)

where $|\Phi_{0}\rangle$ is the reference configuration. The first derivative
of $|\Phi(\mathbf{Z},t)\rangle$ with respect to $t$ at $t=0$ is a linear
combination of singly-excited determinants,

$\begin{split}\left.\frac{\partial|\Phi(\mathbf{Z},t)\rangle}{\partial
t}\right|_{t=0}=Z|\Phi_{0}\rangle=\sum_{a=1}^{N_{\rm vir}}\sum_{i=1}^{N_{\rm
occ}}Z_{ai}|\Phi_{i}^{a}\rangle,\end{split}$ (4)

where $|\Phi_{i}^{a}\rangle=E_{i}^{a}|\Phi_{0}\rangle$. This derivative can be
written as a two-point central difference with finite step size $2\delta t$

$\begin{split}\left.\frac{\partial\left|\Phi(\mathbf{Z},t)\right\rangle}{\partial
t}\right|_{t=0}=\lim_{\delta t\rightarrow
0}\frac{\left|\Phi\left(\mathbf{Z},\delta
t\right)\right\rangle-\left|\Phi\left(\mathbf{Z},-\delta
t\right)\right\rangle}{2\delta t}.\end{split}$ (5)

Therefore, we can approximate a linear combination of
$N_{\mathrm{occ}}N_{\mathrm{vir}}$ singly excited configurations
$Z|\Phi_{0}\rangle$ with only two NOSD

$\begin{split}|\Psi_{\rm S}\rangle=Z|\Phi_{0}\rangle=\lim_{\delta t\rightarrow
0}\frac{\left|\Phi\left(\mathbf{Z},\delta
t\right)\right\rangle-\left|\Phi\left(\mathbf{Z},-\delta
t\right)\right\rangle}{2\delta t}.\end{split}$ (6)

For double excitations, we write the second derivative of
$|\Phi(\mathbf{Z},t)\rangle$ with respect to $t$

$\begin{split}{Z}^{2}|\Phi_{0}\rangle&=\left.\frac{\partial^{2}|\Phi(\mathbf{Z},t)\rangle}{\partial
t^{2}}\right|_{t=0}\\\ &=\lim_{\delta t\rightarrow
0}\frac{|\Phi(\mathbf{Z},2\delta t)\rangle+|\Phi(\mathbf{Z},-2\delta
t)\rangle-2|\Phi_{0}\rangle}{4\delta t^{2}},\end{split}$ (7)

where ${Z}^{2}|\Phi_{0}\rangle$ is approximated by three determinants: two
Thouless-parametrized determinants plus the reference.

A linear combination of double excitations is defined as

$\begin{split}|\Psi_{\text{D}}\rangle=\sum_{\begin{subarray}{c}a,b\in\mathrm{vir}\\\
i,j\in\mathrm{occ}\end{subarray}}w_{ji}^{ba}|\Phi_{ji}^{ba}\rangle,\end{split}$
(8)

where $w_{ji}^{ba}$ are the coefficients, and
$|\Phi_{ji}^{ba}\rangle=E_{j}^{b}E_{i}^{a}|\Phi_{0}\rangle$. Assigning
$p=(b-1)N_{\mathrm{occ}}+j$ and $q=(a-1)N_{\mathrm{occ}}+i$, we define the
matricization of $w_{ji}^{ba}$ as

$\begin{split}W_{pq}=w_{ji}^{ba},\end{split}$ (9)

where $\mathbf{W}$ is a Hermitian matrix of size $L\times L$, with
$L=N_{\mathrm{vir}}N_{\mathrm{occ}}$, and can be chosen real. Hence, we can
diagonalize $\mathbf{W}$ and write

$\begin{split}W_{pq}=\sum_{k=1}^{L}\lambda_{k}U_{kp}U_{kq}=\sum_{k=1}^{L}\lambda_{k}U_{k,bj}U_{k,ai},\end{split}$
(10)

where $\lambda_{k}$ is the $k$th eigenvalue, and rows of $\mathbf{U}$ are
corresponding eigenvectors. Eq. (8) is rewritten as

$\begin{split}|\Psi_{\text{D}}\rangle&=\sum_{k=1}^{L}\lambda_{k}\left(\sum_{j,b}U_{k,bj}{E}_{j}^{b}\right)\left(\sum_{i,a}U_{k,ai}{E}_{i}^{a}\right)|\Phi_{0}\rangle\\\
&=\sum_{k=1}^{L}\lambda_{k}{Z}_{k}^{2}|\Phi_{0}\rangle.\end{split}$ (11)

Combined with Eq. (7), in the limit of $\delta t\rightarrow 0$, we arrive at

$\begin{split}|\Psi_{\text{D}}\rangle=\sum_{k=1}^{L}\lambda_{k}\frac{|\Phi(\mathbf{Z}_{k},2\delta
t)\rangle+|\Phi(\mathbf{Z}_{k},-2\delta t)\rangle-2|\Phi_{0}\rangle}{4\delta
t^{2}}.\end{split}$ (12)

Therefore, we can represent $|\Psi_{\text{D}}\rangle$ with at most
$(2N_{\mathrm{vir}}N_{\mathrm{occ}}+1)$ NOSD, including the reference state
$|\Phi_{0}\rangle$. The compression of double excitations brings a square-root
reduction of the determinant count. In practical applications, introducing a
threshold $\lambda_{\min}$ for $|\lambda_{k}|$ allows exclusion of double
excitations with negligible contributions. This means that the actual number
of $\mathbf{Z_{k}}$ matrices can be significantly less than the theoretical
maximum of $N_{\mathrm{vir}}N_{\mathrm{occ}}$.

A CISD wavefunction is defined as

$\begin{split}\left|\Psi_{\rm
CISD}\right\rangle=w_{0}\left|\Phi_{0}\right\rangle+\sum_{\begin{subarray}{c}a\in\mathrm{vir}\\\
i\in\mathrm{occ}\end{subarray}}w_{i}^{a}\left|\Phi_{i}^{a}\right\rangle+\sum_{\begin{subarray}{c}a,b\in\mathrm{vir}\\\
i,j\in\mathrm{occ}\end{subarray}}w_{ji}^{ba}|\Phi_{ji}^{ba}\rangle,\end{split}$
(13)

where $w_{0}$, $\\{w_{i}^{a}\\}$, and $\\{w_{ji}^{ba}\\}$ are the expansion
coefficients for the reference determinant, single excitations, and double
excitations, respectively.

Following Eq. (6) and (12), we arrive at our non-orthogonal representation of
$\left|\Psi_{\rm CISD}\right\rangle$ in the limit of $\delta t\rightarrow 0$
as

$\begin{split}\left|\Psi_{\rm
CISD}\right\rangle&=\tilde{w}_{0}|\Phi_{0}\rangle\\\ &+\frac{1}{2\delta
t}\left(\left|\Phi(\mathbf{Z}^{(1)},\delta
t)\right\rangle+\left|\Phi(\mathbf{Z}^{(1)},-\delta t)\right\rangle\right)\\\
&+\sum_{k=1}^{L}\frac{\lambda_{k}}{4\delta
t^{2}}\left(|\Phi(\mathbf{Z}_{k}^{(2)},2\delta
t)\rangle+|\Phi(\mathbf{Z}_{k}^{(2)},-2\delta t)\rangle\right),\end{split}$
(14)

where $\tilde{w}_{0}=w_{0}-\sum_{k=1}^{L}\frac{\lambda_{k}}{2\delta t^{2}}$.
The Thouless parameters are given by

$\begin{split}Z_{ai}^{(1)}&=w_{i}^{a},\\\
Z_{k,ai}^{(2)}&=U_{k,ai}.\end{split}$ (15)

In practice, the two-point central difference must be approximated by a
nonzero $\delta t$ value, which determines the accuracy of Eq. (14). Fig. 1(a)
examines the compression accuracy in unrestricted CISD (UCISD) simulations for
the N2 molecule across various $\delta t$ values, employing STO-3G and 6-31G
basis sets, and adjusting the bond length to reflect weak to strong
correlation regimes. We observe that the finite difference error decays with
decreasing $\delta t$, hitting a point where numerical instability arises. An
acceptable $\delta t$ range, highlighted in grey shade, corresponds to a
compression error below $10^{-3}$ Hartree in the total energy. For most
$\delta t$ values between $0$ and $0.5$, the error can be well controlled.
Figure 1(b) illustrates how the number of linearly independent determinants
varies with $\delta t$. With an increase in $\delta t$, we observe a rise in
the number of determinants, eventually reaching a plateau for the total number
of compressed UCISD determinants. In the cases examined in this paper, we
found that $\delta t=0.05$ ensures numerical stability and accuracy.
Therefore, we have chosen to fix $\delta t=0.05$ throughout this work.

(a) Energy error

(b) Number of determinants

Figure 1: Analysis of UCISD compression as a function of $\delta t$ for N2
with STO-3G (solid lines) and 6-31G (dashed lines) at varying bond lengths
($R$). (a) Absolute value of energy error; (b) Number of linearly independent
determinants post-compression. All energies in Hartrees.

(a) STO-3G

(b) 6-31G

Figure 2: Numbers of determinants of UCISD (grey), after compression (dark
blue), after selection based on the metric test (medium blue) and after the
energy test (light blue) for N2 and the H10 chain using (a) STO-3G and (b)
6-31G basis sets at equilibrium ($R_{\rm eq}$) and extended $3R_{\rm eq}$ bond
lengths.

###  II.3 Selection of NOCISD determinants

To further refine the NOCISD expansion, we follow a procedure similar to that
outlined in work by Dutta et. al.: determinants meeting specific criteria are
selected and added to a given set, while others are discarded. Dutta  _et al._
(2021) We employ a metric test with threshold $m_{0}>0$ and an energy test
with threshold $h_{0}>0$. We use Latin letters
$\left|p\right\rangle,\left|q\right\rangle\cdots$ for determinants in the
given set $\mathcal{R}$, and Greek letters
$\left|\mu\right\rangle,\left|\nu\right\rangle\cdots$ for determinants to be
assessed against the criteria.

We start by defining the projector

$\begin{split}{Q}={I}-\sum_{pq}|p\rangle{X}_{pq}\langle q|,\end{split}$ (16)

where $\mathbf{X}=\mathbf{M}^{-1}$ and $M_{pq}=\langle p|q\rangle$. ${Q}$
projects a state onto the orthogonal space of $\mathcal{R}$. The metric test
requires that

$\begin{split}\frac{\left\|{Q}\left|\mu\right\rangle\right\|}{\left\|\left|\mu\right\rangle\right\|}\geq
m_{0}.\end{split}$ (17)

The energy test determines whether a new determinant $|\mu\rangle$ makes a
sufficient contribution to the total energy. Assuming $E_{0}$ is the NOCI
energy of configurations in $\mathcal{R}$, and $\varepsilon$ is the new NOCI
energy after adding $|\mu\rangle$ to $\mathcal{R}$, we require

$\begin{split}\frac{E_{0}-\varepsilon}{|E_{0}|}>h_{0}.\end{split}$ (18)

Details of evaluating $\varepsilon$ are provided in Appendix A.

Figure 3: Impact of varying $\lambda_{\min}$ on the ground state energy (in
Hartrees) and SNOCISD length for the N2 molecule at various bond lengths.
$E_{0}$ is the energy setting $\lambda_{\min}=10^{-7}$.

Figure 2 illustrates the reduction in the number of determinants following
UCISD compression and subsequent selection for two systems: the N2 molecule
and a H10 chain ($10$ H atoms equally distributed in the $z$ direction), using
STO-3G and 6-31G basis sets. We chose $\lambda_{\min}=10^{-7}$, $\delta
t=0.05$, $m_{0}=10^{-5}$ and $h_{0}=10^{-6}$. The equilibrium bond lengths
$R_{\rm eq}$ are computed by selected heat bath configuration interaction
(SHCI).Holmes, Tubman, and Umrigar (2016); Sharma  _et al._ (2017); Li  _et
al._ (2018) The reduction in the number of NOSD after compression and
selection quite is significant compared to the original UCISD expansion. The
reduction is pronounced at larger bond lengths ($3R_{\rm eq}$) regarding the
energy test, illustrating the diminishing role of CISD in strongly correlated
scenarios.

Practically, given the relatively high computational cost of the energy test,
bypassing it can expedite the process. Yet, the energy test is essential for
minimizing the NOCI expansion in strongly correlated systems, presenting a
dilemma: opt for a compact NOCI expansion by conducting the energy test, or
avoid it to reduce computational load at the expense of a larger NOCI
expansion. An alternative to removing less important configurations involves
stricter truncation of the eigenvalue $\lambda_{k}$ in Eq. (14), i.e., using a
larger $\lambda_{\min}$.

Figure 3 shows the impact of varying $\lambda_{\min}$ on the ground state
energy and SNOCISD determinant count for the N2 molecule at various bond
lengths within a 6-31G basis set. The reference NOCI uses two determinants
from a ResHF optimization, setting $\delta t=0.05$ and $m_{0}=10^{-5}$,
skipping the energy test. Energy differences relative to that with
$\lambda_{\min}=10^{-7}$ are presented in the top panel, while the length of
the SNOCISD wavefunction is shown in the bottom panel. Notably, the energy
starts to increase for $\lambda_{\min}>10^{-3}$, while the determinant count
starts to decrease beyond $\lambda_{\min}>10^{-5}$, with an even early drop at
$R=3.0$Å. This trend identifies an optimal region where the choice of
$\lambda_{\min}$ guarantees a compact SNOCISD without compromising the
accuracy. The shorter SNOCISD expansion at $R=3.0$Å corroborates with earlier
findings that fewer determinants pass the energy test at larger bond lengths.
In the following, we bypass the energy test. Instead, meticulous selection of
the $\lambda_{\min}$ value ensures a balance between accuracy and efficiency.

###  II.4 Computational details

A homemade Python package is used to carry out the SNOCISD simulations. The
optimization in FED and ResHF employed the adaptive moment estimation
algorithm (ADAM) Kingma and Ba (2017) implemented in Jax. DeepMind  _et al._
(2020) The optimization process begins with initial guesses that are near-
singly-excited determinants, achieved by setting one Thouless parameter much
larger than the others ($5:0.1$). PySCF Sun  _et al._ (2018, 2020) is used for
quantum chemistry tasks such as integral acquisition, HF, CISD, CCSD, and SHCI
simulations. Block2 Zhai  _et al._ (2023) is used for density matrix
renormalization group (DMRG) calculations.

Figure 4: Dissociation energy curve of N2 molecule. Inset: SNOCISD energy with
respect to the number of reference determinants at $R_{\mathrm{eq}}$. The
energy unit is in Hartree.

##  III Applications

###  III.1 Dissociation energy

Generating accurate dissociation profiles for molecules is a conventional test
for quantum chemistry methods, where systems transition from weak to strong
correlations with bond stretching. This section assesses SNOCISD by computing
the ground-state energy of the N2 molecule across bond lengths from $0.8$ to
$3.0$Å, using the 6-31G basis set, and compares results with SHCI, UCISD, and
UCCSD methods. SHCI serves as the benchmark reference, denoted as "Exact" in
Fig. 4. A half-spin projectionMayhall  _et al._ (2014) is applied to the
SNOCISD wavefunction with $4$ reference configurations generated by FED, shown
in Fig. 4. The compression and selection process yields around 350 NOSD for
each reference without applying the energy test. We fix
$\lambda_{\min}=10^{-5}$, $m_{0}=10^{-5}$, and $\delta t=0.05$.

In the inset plot of Fig. 4, we show the improvement of the SNOCISD energy
with respect to the number of references. The $x$-axis represents the number
of reference configurations and the $y$-axis the energy at the equilibrium
bond length $R_{\rm eq}$. From top to bottom, the horizontal lines correspond
to the UCISD, UCCSD and SHCI energies, respectively. We observe that even with
$3$ reference configurations, the energy is already comparable to the UCCSD
energy. SNOCISD exhibits notable improvements over UCISD and performs
comparably to UCCSD, yet it needs further improvement to achieve FCI-level
accuracy. One promising direction beyond this work is to include higher-order
excitations with careful adaptation of the compression and selection
strategies.

###  III.2 Two-Dimensional Hubbard model

The two-dimensional Hubbard model is a challenging system in electronic
structure studies. Its scientific significance is underscored by applications
in describing the behavior of electrons in real materials, such as Mott
insulators and high-temperature superconductors. Assaad and Imada (1996);
Maier  _et al._ (2005); Sordi, Haule, and Tremblay (2010); Gull, Parcollet,
and Millis (2013); Zheng  _et al._ (2017) The Hamiltonian is given by

$\begin{split}H=-t\sum_{\langle
i,j\rangle,\sigma}\left(a_{i,\sigma}^{\dagger}a_{j,\sigma}+\mathrm{h.c.}\right)+U\sum_{i}n_{i\uparrow}n_{i\downarrow},\end{split}$
(19)

where $t$ represents the hopping amplitude, and $U>0$ is the on-site Coulomb
interaction strength. $\langle i,j\rangle$ represents nearest-neighbor sites,
and $\sigma=\uparrow,\downarrow$ denotes the spin degree of freedom.
Conventionally, $t$ is set to $1$, and energies are expressed in units of $t$.
As $U$ increases, the system transitions from a weakly correlated to a
strongly correlated regime.

####  III.2.1 Ground state energy

Table 1: Ground state energy per site of the two-dimensional Hubbard model of size $4\times 4$ with periodic boundary conditions (PBC) in both directions.  $U/t$ |  | 2 | 4 | 8 | 12
---|---|---|---|---|---
$f=0.875$ | Exact111DMRG results with bond dimension equal to $2000$. | -1.1956 | -0.9805 | -0.7377 | -0.6257
UCISD | -1.1930 | -0.9558 | -0.6252 | -0.5202
UCCSD | -1.1952 | -0.9710 | -0.6717 | -0.5541
SNOCISD | -1.1945 | -0.9682 | -0.6896 | \- 0.5560
$f=0.75$ | Exact | -1.2678 | -1.1077 | -0.9314 | -0.8496
UCISD | -1.2643 | -1.0811 | -0.7927 | -0.7133
UCCSD | -1.2663 | -1.0983 | -0.8673 | -0.7644
SNOCISD | -1.2659 | -1.0932 | -0.8676 | -0.7718

Accurately evaluating the ground state of the two-dimensional Hubbard model
away from half-filling is challenging due to the interplay of competing
phases. LeBlanc  _et al._ (2015) In our study, summarized in Table 1, we
compute the ground state energy per site for a $4\times 4$ Hubbard lattice
under periodic boundary conditions (PBC) in the underdoped regime, comparing
SNOCISD with various numerical methods. The number of spin-up and spin-down
electrons are set to be the same. The SNOCISD simulations used $14$ ResHF-
optimized references to address static correlations, and set $\delta t=0.05$
and $m_{0}=10^{-5}$. For different $U$ values, slightly different
$\lambda_{\min}$ values are picked. In particular, $\lambda_{\min}=0.001$ is
used for $U=2$ and $4$, while $\lambda_{\min}=0.005$ for $U=8$ and $12$. As
previously discussed in Section II.3, in the strongly correlated scenario,
CISD contributes little to the correlation energy. Therefore, it is plausible
to introduce a harsher $\lambda_{\min}$ for larger $U$ values. DMRG results,
with a bond dimension of $2000$, serve as our benchmark, denoted as the exact
value in Table 1. The SNOCISD method not only shows considerable improvement
over UCISD but also outperforms UCCSD at higher $U$ values, showing potential
in accurately capturing the ground state energies in complex correlated
systems.

####  III.2.2 Spin-spin correlation

The existence of broken spin symmetry suggests the presence of degenerate or
near-degenerate spin configurations, where starting from a UHF solution and
gradually adding excited configurations reveals a consistent antiferromagnetic
(AFM) pattern. However, achieving the FCI solution eliminates the AFM order,
as FCI encompasses all configurations. SNOCISD aims to include the low-lying
near-degenerate states, and thus may not precisely capture the long-range AFM
order. Nevertheless, employing SNOCISD for magnetic order studies clarifies
the contribution of each optimized reference state to the overall
wavefunction. The paramagnetic (PM) and AFM phases can be distinguished by the
spin-spin correlation function, $C_{s}(i,j)=\langle
n_{i\uparrow}n_{j\downarrow}\rangle-\langle n_{i\uparrow}\rangle\langle
n_{j\downarrow}\rangle$, which evaluates the spin-spin correlation (or
magnetic correlation) between site-$i$ and site-$j$.

Figure 5: Spin-spin correlation between site-$0$ and site-$i$
($i=0,\cdots,11$) for the two-dimensional Hubbard model of size $2\times 6$
with open boundary conditions (OBC). The numbers in the parentheses represent
the number of reference determinants.

Figure 5 presents the SNOCISD-computed spin-spin correlations for a half-
filled 8×2 Hubbard model with open boundary conditions (OBC). $C_{s}(0,i)$
between the lower-left site (denoted as site-$0$) and other sites
($i=1,\cdots,11$) is computed. The SNOCISD wavefunction employs parameters
$\lambda_{\min}=10^{-5}$, $\delta t=0.05$, and $m_{0}=10^{-5}$. Circle radii
are proportional to the absolute value of $C_{s}(0,i)$ with blue indicating
positive and red negative values. The decaying absolute value with respect to
the distance from site-$0$ and an alternating sign pattern confirms the AFM
order. At a lower interaction strength ($U=2$), a faint AFM order is
noticeable, diminishing with distance. However, at a moderate interaction
strength ($U=6$), the AFM pattern is observed only among nearest neighbors and
disappears at longer distances. This pattern diminishes further with an
increased number of reference determinants (lowest panel). At lower $U$
values, with less spin degeneracy, the optimization of reference determinants
for NOCI energy captures low excitations from the UHF solution, preserving a
weak AFM order. However, at higher $U$ values, the optimization begins to
encompass near-degenerate spin configurations, reducing the bias towards any
broken-symmetry configuration. This hypothesis is confirmed by the reduced
$C_{s}(0,i)$ values as more reference determinants are added.

##  IV Conclusions

In this work, we have presented a systematically improvable approach to
incorporate dynamic correlation into non-orthogonal configuration interaction
(NOCI) wavefunctions, achieved by adding a selected set of singly and doubly
excited configurations. Notably, we prove that any linear combination of
orthogonal excited configurations can be accurately represented by a much
smaller set of non-orthogonal configurations, leading to significant
efficiency gains. Further refinement is achieved through a selection process
based on metric and energy tests. This combined approach, termed selected non-
orthogonal configuration interaction with single- and double-electron
excitations (SNOCISD), offers an efficient method for simultaneously capturing
both dynamic and static correlations.

SNOCISD relies on two key components: compressing the orthogonal CISD
wavefunction and selecting configurations via metric and energy tests. For
compression, we used a two-point central difference strategy for approximating
orthogonal configurations with Thouless-parameterized non-orthogonal
expansions. The metric test assesses the overlap of new determinants with the
reference set, while the energy test uses a simplified energy evaluation
problem. We have analyzed the impact of key parameters, including two-point
difference step size $\delta t$, the eigenvalue cutoff of double excitations
$\lambda_{\min}$, metric threshold $m_{0}$ and the Hamiltonian threshold
$h_{0}$, on the accuracy and computational cost. SNOCISD’s capabilities are
demonstrated by replicating the N2 molecule’s energy dissociation curve and
studying the ground state energy and magnetic correlation of the two-
dimensional Hubbard model. Compared to UCISD, SNOCISD achieved significant
improvements and even surpasses UCCSD in the strong correlation domain.

In summary, SNOCISD provides a promising path towards a computationally
efficient and accurate multireference tool. Nonetheless, further development
is necessary to improve the method’s potential for large-scale simulations and
even higher accuracy. Incorporating higher-order excitations and further
refining the configuration selection processes are key directions for future
studies.

##  Appendix A Selecting determinants based on energy contribution

(a) Metric test.

(b) Energy test.

Figure A.1: Number of NOSDs selected from compressed UCISD (blue lines) and
NOCI energy (red lines) with respect to the values of (a) $m_{0}$ and (b)
$h_{0}$. The system simulated is N2 molecule at $R_{\rm eq}$ with 6-31G. The
energy is shifted by $-109.0$ Hartree.

We derive the expression of $\varepsilon$ in the energy test Eq. (18). We
define the solution in the given set $\mathcal{R}$ as

$\begin{split}|\Psi_{0}\rangle=\sum_{p}w_{p}|p\rangle,\quad
E_{0}=\frac{\langle\Psi_{0}|H|\Psi_{0}\rangle}{\langle\Psi_{0}|\Psi_{0}\rangle}.\end{split}$
(20)

Adding ${Q}|\mu\rangle$ to $|\Psi_{0}\rangle$ results in the wavefunction

$\begin{split}|\Psi\rangle=\eta_{0}|\Psi_{0}\rangle+\eta_{1}{Q}|\mu\rangle.\end{split}$
(21)

The new ground-state energy and corresponding coefficients
$[\eta_{0},\eta_{1}]$ are derived from solving a $2\times 2$ generalized
eigenvalue problem

$\begin{split}\mathbf{Hv}=\varepsilon\mathbf{Sv},\end{split}$ (22)

where

$\begin{split}\mathbf{S}&=\begin{pmatrix}\langle\Psi_{0}|\Psi_{0}\rangle&0\\\
0&\langle\mu|{Q}^{2}|\mu\rangle\end{pmatrix}=\begin{pmatrix}s_{0}&0\\\
0&s_{\mu}\end{pmatrix},\\\
\mathbf{H}&=\begin{pmatrix}\langle\Psi_{0}|{H}|\Psi_{0}\rangle&\langle\Psi_{0}|{H}{Q}|\mu\rangle\\\
\langle\mu|{Q}{H}|\Psi_{0}\rangle&\langle\mu|{Q}{H}{Q}|\mu\rangle\end{pmatrix}=\begin{pmatrix}H_{0}&T\\\
T^{*}&H_{\mu}\end{pmatrix}.\end{split}$ (23)

We require that the new energy after adding ${Q}|\mu\rangle$ passes the energy
test

$\begin{split}\frac{E_{0}-\varepsilon}{|E_{0}|}>h_{0}.\end{split}$ (24)

Here we emphasize that $\varepsilon$ is not the exact NOCI energy of
$|\Psi\rangle=\sum_{p}w_{p}|p\rangle+w_{\mu}|\mu\rangle$, where all
coefficients are relaxed and solved simultaneously. However, $\varepsilon$ is
a good enough approximation to guide the selection. Moreover, in each
selection step, we only need to solve a $2\times 2$ generalized eigenvalue
problem, which is much more economical than solving for the exact NOCI energy.

In Fig. A.1, we profile the truncation error with respect to the values of
$m_{0}$ and $h_{0}$. The set of NOSD to be selected is from compression of
UCISD for the N2 molecule at $R_{\rm eq}$ with 6-31G. We choose $\delta
t=0.05$ and $\lambda_{\min}=10^{-7}$. Fig. A.1(a) shows the number of SDs and
the energy with respect to the value of $m_{0}$, and it suggests that
$m_{0}<10^{-5}$ produces small truncation errors. On the other hand, when
$m_{0}$ is too small, the linear dependency is not fully removed, and
singularity can occur when calculating the NOCI energy. Fig. A.1(b) shows the
number of NOSDs and the energy with respect to the value of $h_{0}$, where we
first removed linear dependency with fixed $m_{0}=10^{-6}$. Understandably,
the energy is more sensitive to the value of $h_{0}$. The error comes from
three aspects: i) $\varepsilon$ being not the exact NOCI energy, ii) error
accumulation, and iii) correlation among discarded determinants. In practice,
one should choose the $h_{0}$ value according to the desired number of SDs and
computing resources.

###### Acknowledgements.

This work was supported by the U.S. Department of Energy, Office of Basic
Energy Sciences, Computational and Theoretical Chemistry Program under Award
DE-SC0001474. G.E.S. acknowledges support as a Welch Foundation Chair (Grant
No. C-0036). This work was also partly supported by the Big-Data Private-Cloud
Research Cyberinfrastructure MRI-award funded by NSF under grant CNS-1338099
and by Rice University’s Center for Research Computing (CRC). C.S. would like
to thank Austin Cheng, Carlos Jimenez-Hoyos, Ruiheng Song, Thomas Henderson,
Xing Zhang, and Zhihao Cui for discussions and suggestions.

## Data Availability Statement

The data that support the findings of this study are openly available on
GitHub with the link:

https://github.com/sunchong137/data_for_nocisd_paper

## References

  * Roos (1980) B. O. Roos, “The complete active space scf method in a fock-matrix-based super-ci formulation,” International Journal of Quantum Chemistry 18, 175–189 (1980).
  * Thom and Head-Gordon (2009) A. J. W. Thom and M. Head-Gordon, “Hartree–Fock solutions as a quasidiabatic basis for nonorthogonal configuration interaction,” The Journal of Chemical Physics 131, 124113 (2009).
  * Jiménez-Hoyos, Rodríguez-Guzmán, and Scuseria (2013) C. A. Jiménez-Hoyos, R. Rodríguez-Guzmán, and G. E. Scuseria, “Multi-component symmetry-projected approach for molecular ground state correlations,” The Journal of Chemical Physics 139, 204102 (2013).
  * Jimenez-Hoyos (2013) C. Jimenez-Hoyos, “Variational approaches to the molecular electronic structure problem based on symmetry-projected hartree–fock configurations,” Ph.D. Thesis, Rice University (2013).
  * Barca, Gilbert, and Gill (2014) G. M. J. Barca, A. T. B. Gilbert, and P. M. W. Gill, “Communication: Hartree-Fock description of excited states of H2,” The Journal of Chemical Physics 141, 111104 (2014).
  * Jensen  _et al._ (2018) K. T. Jensen, R. L. Benson, S. Cardamone, and A. J. W. Thom, “Modeling electron transfers using quasidiabatic hartree–fock states,” Journal of Chemical Theory and Computation 14, 4629–4639 (2018).
  * Mahler and Thompson (2021) A. D. Mahler and L. M. Thompson, “Orbital optimization in nonorthogonal multiconfigurational self-consistent field applied to the study of conical intersections and avoided crossings,” The Journal of Chemical Physics 154, 244101 (2021).
  * Smeyers and Doreste-Suarez (1973) Y. G. Smeyers and L. Doreste-Suarez, “Half-projected and projected hartree-fock calculations for singlet ground states. i. four-electron atomic systems,” International Journal of Quantum Chemistry 7, 687–698 (1973).
  * Smeyers and Delgado-Barrio (1974) Y. G. Smeyers and G. Delgado-Barrio, “Half-projected and projected hartree-fock calculations for singlet ground states. ii. lithium hydride,” International Journal of Quantum Chemistry 8, 733–743 (1974).
  * Jiménez-Hoyos  _et al._ (2012) C. A. Jiménez-Hoyos, T. M. Henderson, T. Tsuchimochi, and G. E. Scuseria, “Projected Hartree–Fock theory,” The Journal of Chemical Physics 136, 164109 (2012).
  * Garza, Jiménez-Hoyos, and Scuseria (2013) A. J. Garza, C. A. Jiménez-Hoyos, and G. E. Scuseria, “Capturing static and dynamic correlations by a combination of projected Hartree-Fock and density functional theories,” The Journal of Chemical Physics 138, 134102 (2013).
  * Tsuchimochi and Ten-no (2016) T. Tsuchimochi and S. Ten-no, “Communication: Configuration interaction combined with spin-projection for strongly correlated molecular electronic structures,” The Journal of Chemical Physics 144, 011101 (2016).
  * Thouless (1960) D. Thouless, “Stability conditions and nuclear rotations in the hartree-fock theory,” Nuclear Physics 21, 225–232 (1960).
  * Fukutome (1988) H. Fukutome, “Theory of Resonating Quantum Fluctuations in a Fermion System: Resonating Hartree-Fock Approximation,” Progress of Theoretical Physics 80, 417–432 (1988).
  * Schmid  _et al._ (1989) K. Schmid, Zheng Ren-Rong, F. Grümmer, and A. Faessler, “Beyond symmetry-projected quasi-particle mean fields: A new variational procedure for nuclear structure calculations,” Nuclear Physics A 499, 63–92 (1989).
  * Schmid (2004) K. Schmid, “On the use of general symmetry-projected hartree–fock–bogoliubov configurations in variational approaches to the nuclear many-body problem,” Progress in Particle and Nuclear Physics 52, 565–633 (2004).
  * Ten-no (1998) S. Ten-no, “Superposition of nonorthogonal slater determinants towards electron correlation problems,” Theoretical Chemistry Accounts 98, 182–191 (1998).
  * Burton and Thom (2020) H. G. A. Burton and A. J. W. Thom, “Reaching full correlation through nonorthogonal configuration interaction: A second-order perturbative approach,” Journal of Chemical Theory and Computation 16, 5586–5600 (2020).
  * Nite and Jiménez-Hoyos (2019) J. Nite and C. A. Jiménez-Hoyos, “Efficient multi-configurational wavefunction method with dynamical correlation using non-orthogonal configuration interaction singles and doubles (nocisd),”  (2019), chemriv:11369646.v1 .
  * Löwdin (1955) P.-O. Löwdin, “Quantum theory of many-particle systems. ii. study of the ordinary hartree-fock approximation,” Phys. Rev. 97, 1490–1508 (1955).
  * Ripka, Blaizot, and Ripka (1986) S. Ripka, J. Blaizot, and G. Ripka,  _Quantum Theory of Finite Systems_ (MIT Press, 1986).
  * Balian and Brezin (1969) R. Balian and E. Brezin, “Nonunitary bogoliubov transformations and extension of wick’s theorem,” Nuovo Cimento B 64, 37–55 (1969).
  * Chen and Scuseria (2023) G. P. Chen and G. E. Scuseria, “Robust formulation of Wick’s theorem for computing matrix elements between Hartree–Fock–Bogoliubov wavefunctions,” The Journal of Chemical Physics 158, 231102 (2023).
  * Dutta  _et al._ (2021) R. Dutta, G. P. Chen, T. M. Henderson, and G. E. Scuseria, “Construction of linearly independent non-orthogonal AGP states,” The Journal of Chemical Physics 154, 114112 (2021).
  * Holmes, Tubman, and Umrigar (2016) A. A. Holmes, N. M. Tubman, and C. J. Umrigar, “Heat-bath configuration interaction: An efficient selected configuration interaction algorithm inspired by heat-bath sampling,” Journal of Chemical Theory and Computation 12, 3674–3680 (2016).
  * 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,” Journal of Chemical Theory and Computation 13, 1595–1604 (2017).
  * Li  _et al._ (2018) J. Li, M. Otten, A. A. Holmes, S. Sharma, and C. J. Umrigar, “Fast semistochastic heat-bath configuration interaction,” The Journal of Chemical Physics 149, 214110 (2018).
  * Kingma and Ba (2017) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” (2017), arXiv:1412.6980 [cs.LG] .
  * DeepMind  _et al._ (2020) DeepMind, I. Babuschkin, K. Baumli, A. Bell, S. Bhupatiraju, J. Bruce, P. Buchlovsky, D. Budden, T. Cai, A. Clark, I. Danihelka, A. Dedieu, C. Fantacci, J. Godwin, C. Jones, R. Hemsley, T. Hennigan, M. Hessel, S. Hou, S. Kapturowski, T. Keck, I. Kemaev, M. King, M. Kunesch, L. Martens, H. Merzic, V. Mikulik, T. Norman, G. Papamakarios, J. Quan, R. Ring, F. Ruiz, A. Sanchez, L. Sartran, R. Schneider, E. Sezener, S. Spencer, S. Srinivasan, M. Stanojević, W. Stokowiec, L. Wang, G. Zhou, and F. Viola, “The DeepMind JAX Ecosystem,”  (2020).
  * Sun  _et al._ (2018) Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, “Pyscf: the python-based simulations of chemistry framework,” WIREs Computational Molecular Science 8, e1340 (2018).
  * Sun  _et al._ (2020) Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, “Recent developments in the PySCF program package,” The Journal of Chemical Physics 153, 024109 (2020).
  * Zhai  _et al._ (2023) H. Zhai, H. R. Larsson, S. Lee, Z.-H. Cui, T. Zhu, C. Sun, L. Peng, R. Peng, K. Liao, J. Tölle, J. Yang, S. Li, and G. K.-L. Chan, “Block2: A comprehensive open source framework to develop and apply state-of-the-art DMRG algorithms in electronic structure and beyond,” The Journal of Chemical Physics 159, 234801 (2023).
  * Mayhall  _et al._ (2014) N. J. Mayhall, P. R. Horn, E. J. Sundstrom, and M. Head-Gordon, “Spin–flip non-orthogonal configuration interaction: a variational and almost black-box method for describing strongly correlated molecules,” Phys. Chem. Chem. Phys. 16, 22694–22705 (2014).
  * Assaad and Imada (1996) F. F. Assaad and M. Imada, “Insulator-metal transition in the one- and two-dimensional hubbard models,” Phys. Rev. Lett. 76, 3176–3179 (1996).
  * Maier  _et al._ (2005) T. A. Maier, M. Jarrell, T. C. Schulthess, P. R. C. Kent, and J. B. White, “Systematic study of $d$-wave superconductivity in the 2d repulsive hubbard model,” Phys. Rev. Lett. 95, 237001 (2005).
  * Sordi, Haule, and Tremblay (2010) G. Sordi, K. Haule, and A.-M. S. Tremblay, “Finite doping signatures of the mott transition in the two-dimensional hubbard model,” Phys. Rev. Lett. 104, 226402 (2010).
  * Gull, Parcollet, and Millis (2013) E. Gull, O. Parcollet, and A. J. Millis, “Superconductivity and the pseudogap in the two-dimensional hubbard model,” Phys. Rev. Lett. 110, 216405 (2013).
  * Zheng  _et al._ (2017) B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, “Stripe order in the underdoped region of the two-dimensional hubbard model,” Science 358, 1155–1160 (2017).
  * LeBlanc  _et al._ (2015) J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), “Solutions of the two-dimensional hubbard model: Benchmarks and results from a wide range of numerical algorithms,” Phys. Rev. X 5, 041041 (2015).