# Further Development of iCIPT2 for Strongly Correlated Electrons

Ning Zhang  Beijing National Laboratory for Molecular Sciences, Institute of
Theoretical and Computational Chemistry, College of Chemistry and Molecular
Engineering, Peking University, Beijing 100871, China  Wenjian Liu
<EMAIL_ADDRESS>Qingdao Institute for Theoretical and Computational Sciences,
Shandong University, Qingdao, Shandong 266237, China  Mark R. Hoffmann
Chemistry Department, University of North Dakota, Grand Forks, ND 58202-9024,
U.S.A.

###### Abstract

The efficiency of the recently proposed iCIPT2 [iterative configuration
interaction (iCI) with selection and second-order perturbation theory (PT2);
J. Chem. Theory Comput. 16, 2296 (2020)] for strongly correlated electrons is
further enhanced (by up to $20\times$) by using (1) a new ranking criterion
for configuration selection, (2) a new particle-hole algorithm for Hamiltonian
construction over randomly selected configuration state functions (CSF), and
(3) a new data structure for the quick sorting of the joint variational and
first-order interaction spaces. Meanwhile, the memory requirement is also
reduced greatly. As a result, this improved implementation of iCIPT2 can
handle one order of magnitude more CSFs than the previous version, as revealed
by taking the chromium dimer and an iron-sulfur cluster, [Fe2S2(SCH3)]42-, as
examples.

##  1 Introduction

As pointed out recently1, 2, among the three components of electronic
structure theory (i.e., relativity, correlation, and QED), it is correlation
that is most challenging, especially for strongly correlated systems. Roughly
speaking, a strongly correlated system features multiple open-shell orbitals
or nearly degenerate electronic states thanks to the existence of a dense set
of energetically adjacent frontier orbitals, such that there does not exist a
single, leading component in the wave functions. Because of this, single-
reference methods are bound to fail, especially for those low-spin states.
Instead, multi-reference methods3, 4, 5, 6, 7, 8 are needed even for a
qualitative description of such states. A unanimous assumption underlying
standard multi-reference methods is that the full many-electron Hilbert space
can be decomposed into a model space $P$ and the remaining complement $Q=1-P$,
which are responsible for the static and dynamic components of the overall
correlation, respectively. Then, according to when the static and dynamic
correlations are handled, such methods can be classified into three families,
viz. “static-then-dynamic”, “dynamic-then-static”, and “static-dynamic-static”
(SDS)9. Given their great success in investigating low-lying electronic states
of many chemical systems, it must be realized that such multi-reference
methods still have strong limitations. For instances, it is by no means
trivial to maintain the same model space, so as to produce smooth potential
energy surfaces. Even if this can be achieved, there is no guarantee that the
chosen model space is equally good for all geometries and all target states.
Even more seriously, those states that have little projections in the chosen
model space cannot well be described. Taking all these facts into account, one
can classify chemical systems into weakly correlated (which can well be
described by single-reference methods), moderately correlated (which can well
be described by standard multi-reference methods with a moderate model space,
typically much smaller than the complete active space (CAS) of 18 electrons in
18 orbitals), and strongly correlated (which require a large model space and
are therefore beyond the capability of standard multi-reference methods). As a
matter of fact, the concept of model space looses its original meaning for the
so-defined strongly correlated systems because the states in a large model
space $P$ and those formally in the orthogonal complement $Q$ are heavily
intersected in both energy and composition. In other words, the static and
dynamic correlations are in such cases strongly entangled and interchangeable.
Moreover, solving the large-model-space problem itself is already a heavy
task, needless to say the subsequent treatment of dynamic correlation (NB: the
larger the model space, the larger the first-order interacting space (FOIS)).
So the question is how to design an adaptive method that can adapt to the
variable static correlation automatically and can meanwhile handle the
residual dynamic correlation efficiently. It appears that the most promising
way to this end is to introduce some selection procedure, so as to build up
progressively a compact yet accurate variational space $P$, in full accordance
with the nature of the system under concern. In this spirit, a number of near-
exact approaches have been designed in the past, including density matrix
renormalization group (DMRG)10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, full
configuration interaction (FCI) quantum Monte Carlo (FCIQMC)21, 22, 23, 24,
25, 26, 27, 28, 29, 30, cluster-analysis-driven FCIQMC31, adaptive coupled-
cluster (CC)32, full CC reduction33, 34, deterministically projected FCI35,
many-body expanded FCI36, 37, 38, 39, incremental FCI40, 41, rank-reduction
FCI42, fast randomized interaction-based FCI43, intrinsic scaling-based
correlation expansion44, 45, as well as selected CI (sCI) 46, 47, 48, 49, 50,
51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 9, 65, 66, 67, 68, 69,
70. Among these, the sCI type of approaches, which are largely a revival of
the very old ideas71, 72, 73, 74, 75 in one way or another, are the simplest
in structure and rival those more sophisticated approaches as shown in a
recent blind test76. As the newest member of the sCI family, the iCIPT2
approach66 proposed by the present authors, i.e., the combination of iterative
CI (iCI)9 with selection and second-order perturbation theory (PT2) has a
number of salient features:

  1. (1)

Born from the restricted SDS framework64 for strongly correlated electrons,
iCI9 is a parameter-free, exact solver of the FCI eigenvalue problem. It
constructs and diagonalizes a very small Hamiltonian matrix $\bar{\mathbf{H}}$
(see Sec. 2) at each macro/micro-iteration but can converge quickly to the
exact solutions. The major bottleneck of iCI lies therefore only in the
construction of the Hamiltonian matrix $\mathbf{H}$ in the basis of Slater
determinants (SD) or configuration state functions (CSF) before contracted to
$\bar{\mathbf{H}}$ for diagonalization. Since $\mathbf{H}$ is extremely
sparse, the combination of iCI with selection is a natural choice, leading to
“selected iCI” (SiCI) as an near-exact solver of FCI. To be more practical,
the (one-step) SiCI can be performed in two steps, with the first step
accounting for the static correlation whereas the second step for a second-
order perturbative treatment of the dynamic correlation. This can be achieved
by a single parameter $C_{min}$, which controls the size of the variational
space and hence the final accuracy. This scheme has been dubbed as iCIPT266.

  2. (2)

The selection is carried out over the entire Hilbert space or a very large
active space (which can in this context be viewed as a coarse-grained
selection step). Those unimportant doubly excited orbital configurations (CFG)
are first screened out based on the upper bounds of the elements of
$\mathbf{H}$. After this integral-driven screening, the individual CSFs of a
selected CFG are further selected using their approximate first-order
coefficients. In short, although the selection is performed on individual
CSFs, it is CFGs that are used as the organizing units. Given a coefficient
pruning-threshold $C_{min}$, the selection of important CFGs/CSFs is performed
iteratively until convergence.

  3. (3)

At each iteration for the growth of wave function, the FOIS is decomposed into
disjoint subspaces, so as to reduce memory requirement on one hand and
facilitate parallelization on the other.

  4. (4)

A new technique called TUGA (tabulated unitary group approach) is introduced
to compute and reuse the basic coupling coefficients (BCC), i.e., the matrix
elements of $U(n)$ generators over CSF pairs.

  5. (5)

Because of the use of CSFs as the many-electron basis, full spin symmetry is
always maintained, which is of vital importance for describing low-spin states
of general open-shell systems as well as subsequent perturbative treatment of
spin-orbit couplings.

  6. (6)

Upon termination of the selection, dynamic correlation is estimated by using
the state-specific Epstein-Nesbet type of second-order perturbation theory
(ENPT2).

  7. (7)

The linear relationship between the total versus second-order correlation
energies allows for an accurate extrapolation.

The present contribution amounts to further enhancing the efficiency of iCIPT2
by using (a) a new ranking criterion for the selection of CFGs/CSFs (see Sec.
4), (b) a new particle-hole algorithm for the construction and update of
$\mathbf{H}$ in the variational space (see Sec. 5), and (c) the Timsort
algorithm77, 78 for sorting the joint variational and first-order interacting
spaces (see Sec. 6). Before these, we ought to outline the SDS family of
methods in Sec. 2 and the essentials of TUGA in Sec. 3. A close comparison of
TUGA with Table-CI79, 80, 81, 82 is also given in the latter, since they are
both configuration-driven alogithms. The chromium dimer and an iron-sulfur
cluster, [Fe2S2(SCH3)]42-, are finally taken as examples to reveal the
efficacy of the improved implementation of iCIPT2 (see Sec. 7). The
presentation is closed with concluding remarks in Sec. 8.

##  2 The SDS family of methods: SDSCI, SDSPT2, iCI, iCIPT2

The _restricted_ “static-dynamic-static” (SDS) framework64 for strongly
correlated electrons gives rise first to the SDSCI approach, which is a
minimal MRCI with the following form for the wave function $|\Psi_{i}\rangle$
of state $i\in[1,N_{P}]$,

$\displaystyle|\Psi_{i}\rangle$
$\displaystyle=\sum_{k=1}^{mN_{P}}|\tilde{\Phi}_{k}\rangle C_{k,i},\quad
m=3\mbox{ or }4,$ (1) $\displaystyle|\tilde{\Phi}_{k}\rangle$
$\displaystyle=|\Psi_{k}^{(0)}\rangle=\sum_{|J\nu\rangle\in
P}|J\nu\rangle\bar{C}_{\nu,k}^{J(0)},\quad k\in[1,N_{P}],$ (2)
$\displaystyle|\tilde{\Phi}_{k+N_{P}}\rangle$
$\displaystyle=|\Psi_{k}^{(1)}\rangle=Q\frac{1}{E_{k}^{(0)}-H_{0}}QH|\Psi_{k}^{(0)}\rangle=\sum_{|J\nu\rangle\in
Q}|J\nu\rangle\bar{C}_{\nu,k}^{J(1)},\quad k\in[1,N_{P}],$ (3)
$\displaystyle|\tilde{\Phi}_{k+2N_{P}}\rangle$
$\displaystyle=P_{\mathrm{s}}|\Theta_{k}\rangle,\quad k\in[1,N_{P}],$ (4)
$\displaystyle|\tilde{\Phi}_{k+3N_{P}}\rangle$
$\displaystyle=P_{\mathrm{s}}|\Theta_{k}^{\prime}\rangle,\quad
k\in[1,N_{P}]\quad(\mbox{ if }m>3),$ (5) $\displaystyle P$
$\displaystyle=\sum_{J\nu}^{d_{R}}|J\nu\rangle\langle
J\nu|=P_{\mathrm{m}}+P_{\mathrm{s}},\quad
P_{\mathrm{m}}=\sum_{k=1}^{N_{P}}|\Psi_{k}^{(0)}\rangle\langle\Psi_{k}^{(0)}|,\quad
Q=1-P.$ (6)

Here, the first two sets of $N_{P}$ functions $\\{|\tilde{\Phi}_{k}\rangle\\}$
are nothing but the zeroth-order (primary) and first-order (external)
functions, respectively. To account for changes in the static correlation
(described by $\\{|\Psi_{k}^{(0)}\rangle\\}_{k=1}^{N_{P}}$) due to the
presence of dynamic correlation (described by
$\\{|\Psi_{k}^{(1)}\rangle\\}_{k=1}^{N_{P}}$), $N_{P}$ secondary (buffer)
functions $\\{|\Theta_{k}\rangle\\}_{k=1}^{N_{P}}$, Eq. (4), are further
introduced, which can be defined in a number ways. For instance, the not-
energy-biased Lanczos-type functions (denoted as Buf(1)),
$\\{H|\Psi_{k}^{(1)}\rangle\\}_{k=1}^{N_{P}}$, turn out to be very
effective64. An alternative choice might be the second set of $N_{P}$ zeroth-
order functions $\\{|\Psi_{k}^{(0)}\rangle\\}_{k=N_{P}+1}^{2N_{P}}$ (denoted
as Buf(0)) in the spirit of intermediate Hamiltonian theory, the first set of
$N_{P}$ second-order functions $\\{|\Psi_{k}^{(2)}\rangle\\}_{k=1}^{N_{P}}$)
(denoted as Buf(2)) in the spirit of perturbation theory or the $N_{P}$
approximate eigenvectors of the preceding iteration (denoted as Buf(3)) in the
spirit of conjugate gradient theory83. While Buf(2) is very expensive, both
Buf(0) and Buf(3) are very easy to handle because the Hamiltonian matrix
elements can readily be evaluated83. As a matter of fact, Buf(0) and Buf(3)
can be applied simultaneously without increasing the computational overhead.
In this case, we will have $m=4$ in Eq. (1). Otherwise, we have $m=3$ in Eq.
(1) with Buf(1) alone. The yet unknown expansion coefficients are to be
determined by the generalized secular equation

$\displaystyle\bar{\mathbf{H}}\mathbf{C}=\mathbf{S}\mathbf{C}\mathbf{E}.$ (7)

Since all the functions $\\{|\tilde{\Phi}_{k}\rangle\\}_{k=N_{P}+1}^{mN_{P}}$
are fully contracted (i.e., linear combinations of the CSFs
$\\{|\nu\rangle\\}$ generated from CFGs $\\{|J\rangle\\}$) and are specific to
the primary state $|\Psi_{k}^{(0)}\rangle$, the dimension of Eq. (7) is just 3
(or 4) times the number ($N_{P}$) of target states, no matter how many
electrons and how many orbitals are correlated. Although restricted as such,
it has been shown84 that SDSCI is a very effective variational method for both
ground and excited states. Two extensions of SDSCI have been considered so
far, SDSPT264, 65 and iCI9, both with Buf(1) alone. The former amounts to
replacing the $QHQ$ block of the Hamiltonian matrix in Eq. (7) with $QH_{0}Q$.
Different from most variants of multi-reference second-order perturbation
theory (MRPT2), the CI-like SDSPT2 treats single and multiple states in the
same way and is particularly advantageous when a number of states are nearly
degenerate, thanks to the efficacy of the secondary states in revising the
coefficients of the primary states. In contrast, iCI takes the solutions of
Eq. (7) as new primary states $\\{|\Psi_{k}^{(0)}\rangle\\}_{k=1}^{N_{P}}$ and
repeats the SDS procedure (1) until convergence. Clearly, each iteration
(defined as macro-iteration) accesses a space that is higher by two ranks than
that of the preceding iteration. Up to $2i$-tuple excitations (relative to the
initial primary space) can be accessed if $i$ macro-iterations are carried
out. A few micro-iterations $j$ can be invoked at each macro-iteration so as
to relax the contraction coefficients. As such, every iCI($i,j$) corresponds
to a physically meaningful model. It has been shown both theoretically and
numerically that iCI can converge quickly from above to FCI even when starting
with a very poor initial guess. More generally, iCI can be viewed as a
particular sequential, exact partial diagonalization of a huge matrix, by
getting first the roots of one portion of the matrix and then those of an
enlarged portion, until the full matrix, whereas its micro-iterations can be
generalized to an iterative vector interaction (iVI) approach83 for the roots
of a given matrix treated as a whole. Again, each iteration of iVI constructs
and diagonalizes a $mN_{P}\times mN_{P}$ matrix ($m=3\mbox{ or }4$) in the
form of Eq. (7). In particular, by combining with the energy-vector following
technique, iVI can directly access interior roots belonging to a predefined
window, without knowing the number and characters of the roots85.

As stated before, iCIPT2 amounts to splitting iCI into two steps, the first of
which accounts for the static correlation by selecting important CFGs/CSFs in
an iterative manner, whereas the second of which accounts for the residual
dynamic correlation via a second-order perturbation theory. It is the
improvement of iCIPT2 that is under concern here.

##  3 TUGA

The very first issue pertinent to iCIPT2 is how to evaluate the Hamiltonian
matrix elements over randomly selected CSFs $\\{|J\nu\rangle\\}$. In this
context, the configuration-driven Table-CI approach79, 80, 81, 82 may be a
possible choice. However, it involves line-up permutations of open-shell
orbitals, which become ineffective for CFGs of many open-shell orbitals (say,
$>10$). To avoid such line-up permutations, we have introduced a “tabulated
unitary group approach” (TUGA)66 (originally called “tabulated orbital
configuration-based unitary group approach”). Since the UGA86, 87, 88 is
itself rather involved, we here recapitulate the essentials of TUGA. The first
step is to break the spin-free, second-quantized Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{i,j}h_{ij}E_{ij}+\frac{1}{2}\sum_{i,j,k,l}(ij|kl)e_{ij,kl},$
(8) $\displaystyle E_{ij}$ $\displaystyle=$
$\displaystyle\sum_{\sigma}a_{i\sigma}^{\dagger}a_{j\sigma}=E_{ji}^{\dagger},$
(9) $\displaystyle e_{ij,kl}$ $\displaystyle=$
$\displaystyle\sum_{\sigma,\tau}a_{i\sigma}^{\dagger}a_{k\tau}^{\dagger}a_{l\tau}a_{j\sigma}=E_{ij}E_{kl}-\delta_{jk}E_{il}=\\{E_{ij}E_{kl}\\}$
(10) $\displaystyle=$
$\displaystyle\\{E_{kl}E_{ij}\\}=e_{kl,ij}=e^{\dagger}_{ji,lk}$

into a form that is consistent with the diagrams employed in the UGA87 for the
BCC between CSFs. Specifically,

$\displaystyle H$ $\displaystyle=$ $\displaystyle
H_{1}^{0}+H_{2}^{0}+H_{1}^{1}+H_{2}^{1}+H_{2}^{2},$ (11) $\displaystyle
H_{1}^{0}$ $\displaystyle=$ $\displaystyle\sum_{i}h_{ii}E_{ii},$ (12)
$\displaystyle H_{1}^{1}$ $\displaystyle=$
$\displaystyle\sum_{i<j}h_{ij}E_{ij}+\sum_{i>j}h_{ij}E_{ij},$ (13)
$\displaystyle H_{2}^{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i}(ii|ii)E_{ii}(E_{ii}-1)+\sum_{i<j}\left[(ii|jj)E_{ii}E_{jj}+(ij|ji)e_{ij,ji}\right],$
(14) $\displaystyle H_{2}^{1}$ $\displaystyle=$ $\displaystyle\sum_{i\neq
j}\left[(ii|ij)(E_{ii}-1)E_{ij}+(ij|jj)E_{ij}(E_{jj}-1)\right]$ (15)
$\displaystyle+$ $\displaystyle\sum_{i\neq j\neq
k}\left[(ij|kk)E_{ij}E_{kk}+(ik|kj)e_{ik,kj}\right],$ $\displaystyle
H_{2}^{2}$ $\displaystyle=$ $\displaystyle\sum_{i\leq k}\sum_{j\leq
l}^{\prime}\left[2^{-\delta_{ik}\delta_{jl}}(ij|kl)e_{ij,kl}+(1-\delta_{ik})(1-\delta_{jl})(il|kj)e_{il,kj}\right],$
(16)

where the superscripts 0, 1, and 2 in $H_{i}$ ($i=1,2$) indicate that the
terms contribute to the Hamiltonian matrix elements over two CFGs that are
related by zero, single, and double excitations, respectively, whereas the
prime in Eq. (16) indicates that $\\{j,l\\}\cap\\{i,k\\}=\emptyset$. For two
identical CFGs $|I\rangle$, the Hamiltonian matrix elements can be calculated
most efficiently as66

$\displaystyle\langle I\mu|H_{1}^{0}+H_{2}^{0}|I\nu\rangle$
$\displaystyle=\delta_{\mu\nu}\\{\frac{1}{2}\sum_{i}n_{i}^{R}[h_{ii}+\epsilon_{i}+(n_{i}^{R}-2)g_{ii}]+\sum_{i}\Delta_{i}^{I}[\epsilon_{i}+(n_{i}^{R}-1)g_{ii}]$
$\displaystyle+\sum_{i\leq
j}\Delta_{i}^{I}g_{ij}\Delta_{j}^{I}\\}+\sum_{i<j}(ij|ji)\langle
I\mu|e_{ij,ji}^{1}|I\nu\rangle\delta(n_{i}^{I},1)\delta(n_{j}^{I},1),$ (17)
$\displaystyle f_{ij}$
$\displaystyle=h_{ij}+\sum_{k}n_{k}^{R}[(ij|kk)-\frac{1}{2}(ik|kj)],$ (18)
$\displaystyle\epsilon_{i}$
$\displaystyle=f_{ii}=h_{ii}+\sum_{k}n_{k}^{R}g_{ik},\quad
g_{ik}=(ii|kk)-\frac{1}{2}(ik|ki)=g_{ki},$ (19) $\displaystyle\Delta_{i}^{I}$
$\displaystyle=n_{i}^{I}-n_{i}^{R},$ (20)

where $\\{n_{i}^{R}\\}$ are the occupation numbers of the spatial orbitals of
a common reference CFG $|R\rangle$ (either Hartree-Fock (HF) or restricted
open-shell HF). The superscript 0 or 1 in $e_{ij,kl}$ (see Eq. (17)) refers to
the intermediate angular momentum $X$ due to the coupling of two
spin-$\frac{1}{2}$ functions. If $|I\rangle=E_{ij}|J\rangle$, the
corresponding Hamiltonian matrix elements read (under the Yamanouchi-Kotani
(YK) phase87)

$\begin{split}\langle I\mu|H_{1}^{1}+H_{2}^{1}|J\nu\rangle&=\left\langle
I\mu|E_{ij}|J\nu\right\rangle\left[f_{ij}+\sum_{k}\Delta_{k}^{J}[(ij|kk)-\frac{1}{2}(ik|kj)]\right.\\\
&+\left.\frac{1}{2}n_{i}^{J}(ii|ij)+(\frac{1}{2}n_{j}^{J}-1)(ij|jj)\right]\\\
&+\sum_{k\in\mbox{ exterior open}}(ik|kj)\langle
I\mu|e_{ik,kj}^{1}|J\nu\rangle\\\ &+\sum_{k\in\mbox{ interior
open}}(ik|kj)\left[\frac{1}{2}\left\langle
I\mu|E_{ij}|J\nu\right\rangle+\langle
I\mu|E_{kj}E_{ik}|J\nu\rangle\right],\end{split}$ (21)

where ‘exterior open’ and ‘interior open’ emphasize that level $k$ belongs to
the non-overlapping ($k<\min(i,j)$ or $k>\max(i,j)$) and overlapping
($\min(i,j)<k<\max(i,j)$) cases, respectively, and is singly occupied. On the
other hand, if $|I\rangle=e_{ij,kl}|J\rangle$ subject to
$\\{j,l\\}\cap\\{i,k\\}=\emptyset$, the Hamiltonian matrix elements read

$\displaystyle\langle
I\mu|H_{2}^{2}|J\nu\rangle=\left[2^{-\delta_{ik}\delta_{jl}}(ij|kl)\langle
I\mu|e_{ij,kl}|J\nu\rangle+(1-\delta_{ik})(1-\delta_{jl})(il|kj)\langle
I\mu|e_{il,kj}|J\nu\rangle\right].$ (22)

At this stage two important points can be observed:

  1. (1)

Thanks to the conjugacy (bra-ket inversion) relations $\langle
I\mu|E_{ij}|J\nu\rangle=\langle I\mu|E_{ij}|J\nu\rangle^{*}=\langle
J\nu|E_{ji}|I\mu\rangle$ and $\langle I\mu|e_{ij,kl}|J\nu\rangle=\langle
I\mu|e_{ij,kl}|J\nu\rangle^{*}=\langle J\nu|e_{ji,lk}|I\mu\rangle$ in the
absence of spin-orbit couplings, only generators $E_{ij}$ with $i>j$ and
$e_{ij,kl}$ with $k\geq i$ and $k>l\geq j$ need to be considered explicitly.
They correspond to the s2, c$x$ and d$x$ ($x\in[1,7]$) types of diagrams shown
in Figs. 2 and 3 of Ref. 66. More specifically, the s2, c4, c6 and d2 diagrams
are required for the evaluation of Eq. (21), while the other c$x$ and d$x$
types of diagrams are needed for the evaluation of Eq. (22). Once such BCCs
are available, those conjugate ones can be obtained simply by matrix
transpose.

  2. (2)

The BCCs $\langle I\mu|E_{ij}|J\nu\rangle$ and $\langle
I\mu|e_{ij,kl}|J\nu\rangle$ depend only on the relative occupations of the CFG
pair $(I,J)$ but not on the individual orbitals $\\{i,j,k,l\\}$. Moreover, the
common zero or doubly occupied orbitals need not be considered explicitly
because their segment values are just one under the YK phase. Therefore, they
can be rewritten as $\langle I\mu|E_{\bar{i}\bar{j}}|J\nu\rangle$ and $\langle
I\mu|e_{\bar{i}\bar{j},\bar{k}\bar{l}}|J\nu\rangle$, respectively, in terms of
the reduced orbital indices (ROI) $\\{\bar{i},\bar{j},\bar{k},\bar{l}\\}$
after deleting the common zero or doubly occupied orbitals. Note that the ROIs
have one-to-one correspondence with the original orbital indices (see the
example in Table 1). This way, the same BCCs can be used for very many
integrals sharing the same ROIs.

After deleting the common zero or doubly occupied orbitals in the bra and ket
CFGs, we are left with seven occupation patterns, each of which can be
assigned a code number as shown in Table 2. Every pair of CFGs can therefore
be characterized by a reduced occupation table (ROT) consisting of two
sequences, ROT_Orb and ROT_Code. The former records the orbital indices to
fetch molecular integrals, whereas the latter records the corresponding code
sequence (e.g., (012430) for the example shown in Table 1) to determine the
generators, with the following rules: (1) if CFG $|I\rangle$ arises from
$|J\rangle$ by exciting one electron from orbital $j$ to orbital $i$,
$\bar{j}$ would then correspond to code 2 or 4, while $\bar{i}$ to code 1 or 3
because of the relations $n_{j}^{I}=n_{j}^{J}-1$ and $n_{i}^{I}=n_{i}^{J}+1$.
(2) If CFG $|I\rangle$ arises from $|J\rangle$ by exciting two electrons from
orbitals $j$ and $l$ $(\geq j)$ to orbitals $i$ and $k$ $(\geq i)$ [NB:
$\\{i,k\\}\cap\\{j,l\\}=\emptyset$], $\bar{j}$ and $\bar{l}$ $(\geq\bar{j})$
would correspond to code 2, 4 or 6, while $\bar{i}$ and $\bar{k}$
$(\geq\bar{i})$ to code 1, 3 or 5 because of the relations
$n_{j}^{I}=n_{j}^{J}-1$, $n_{l}^{I}=n_{l}^{J}-1$, $n_{i}^{I}=n_{i}^{J}+1$ and
$n_{k}^{I}=n_{k}^{J}+1$. Codes 6 and 5 further imply $\bar{j}=\bar{l}$ and
$\bar{i}=\bar{k}$, respectively. (3) If $\bar{i}<\bar{j}$ in
$E_{\bar{i}\bar{j}}$ or $\bar{k}<\bar{l}$ in
$e_{\bar{i}\bar{j},\bar{k}\bar{l}}$, a bra-ket inversion should be invoked
when calculating the BCCs in terms of the diagrams documented in Table 3. To
search and sort the code sequences efficiently, each ROT_Code will be
converted to an array of 64-bit integers (3 bits per code). It is clear that,
to reutilize the BCCs efficiently, the CFG pairs with the same ROT_Code must
be grouped together. Instead of the red-black tree type of bilinear search
employed before66, we adopt here an array-based sorting algorithm: (a)
generate all connected CFG pairs and store them in one array; (b) loop over
each CFG pair and identify its ROT_Orb and ROT_Code; (c) sort the array based
on ROT_Code. This way, those CFG pairs sharing the same ROT_Code (and hence
the same BCCs) are adjacent to each other and form a segment in the array. To
expedite the sorting process, the code sequences can further be characterized
by their lengths and category numbers (CN) defined in Table 3. The latter just
number different ranges of the generators ($\mathrm{CN}=0$ for singles and
$\mathrm{CN}\in[1,7]$ for doubles) and are hence directly related to the
diagrams required for the evaluation of the BCCs. This new sorting algorithm
is significantly faster than the red-black trees employed before, although the
latter is not really expensive for this purpose.

Table 1: Illustration on the reduced orbital indices (ROI)

OrbIndx |  0 |  1 |  2 |  3 |  4 |  5 |  6 |  7
---|---|---|---|---|---|---|---|---
bra occ |  1 |  0 |  1 |  0 |  2 |  1 |  2 |  1
ket occ |  1 |  0 |  0 |  1 |  2 |  2 |  1 |  1
ROT_Orb |  0 |  $\times$ |  2 |  3 |  $\times$ |  5 |  6 |  7
ROI |  $\bar{0}$ |  $\times$ |  $\bar{1}$ |  $\bar{2}$ |  $\times$ |  $\bar{3}$ |  $\bar{4}$ |  $\bar{5}$
ROT_Code* |  0 |  $\times$ |  1 |  2 |  $\times$ |  4 |  3 |  0

  * *

See Table 2.

Table 2: Code numbers for the relative occupation patterns of CFG pairs bra occ |  1 |  1 |  0 |  2 |  1 |  2 |  0 |  0 |  2
---|---|---|---|---|---|---|---|---|---
ket occ |  1 |  0 |  1 |  1 |  2 |  0 |  2 |  0 |  2
code |  0 |  1 |  2 |  3 |  4 |  5 |  6 |  $\times$ |  $\times$
Table 3: Correspondence between UGA diagrams and category numbers (CN) of code
sequences recorded in ROT_code

CN | Generator | Range | Diagrama |  $\tilde{H}^{IJ}$ b
---|---|---|---|---
0 | $E_{ij}$ | $i>j$ | s2 |
$e_{ik,kj}^{1}$ | $i>j>k$ | c4 |
$k>i>j$ | c6 |
$i>k>j$ | d2 |
1 | $e_{ij,kl}$ | $k>l>j>i$ | c1 |  $\max(|(ij|kl)+(il|kl)|$, $\sqrt{3}|(ij|kl)-(il|kj)|)$
$e_{il,kj}$ | c3
2 | $e_{ij,kl}$ | $k>l>i>j$ | d1
$e_{il,kj}$ | c5
3 | $e_{ij,kl}$ | $k>i>l>j$ | d3
$e_{il,kj}$ | d5
4 | $e_{ij,kj}$ | $k>j>i$ | c2 | 2$|(ij|kj)|$
5 | $e_{ij,kj}$ | $k>i>j$ | d4
6 | $e_{il,ij}$ | $i>l>j$ | d6 | 2$|(il|ij)|$
7 | $e_{ij,ij}$ | $i>j$ | d7 | $|(ij|ij)|$

  * a

See Figs. 1-3 in Ref. 66.

  * b

Estimates of $|\langle I\mu|H_{2}^{2}|J\nu\rangle|$, which are upper bounds66
for $\mathrm{CN}\in[4,7]$, exact for $\mathrm{CN}\in[1,3]$ if $|I\rangle$ or
$|J\rangle$ is a closed shell CFG, but not guaranteed to be upper bounds for
$\mathrm{CN}\in[1,3]$ if both $|I\rangle$ and $|J\rangle$ are open-shell CFGs
(NB: slight deviations from the upper bounds can be taken care of by adjusting
$C_{min}$).

To be compatible with the desired generator ranges, we have to represent and
store properly the CFGs and corresponding CSFs. A CFG can first be represented
by an array of occupation numbers, Occ[$i$]=$n_{i}$, $i\in[0,N_{orb}-1]$, with
$N_{orb}$ being the total number of spatial orbitals. The comparison of a CFG
pair can then be performed efficiently with Algorithm 3 presented
previously66. A CFG pair $(I,J)$ is said to be singly or doubly connected if
$|I\rangle$ can be obtained by exciting one or two electrons from $|J\rangle$.
CFG $|I\rangle$ is stored preceding CFG $|J\rangle$ (i.e., $I<J$) if and only
if $n_{p}^{I}>n_{p}^{J}$ for $p=\max\\{x|n_{x}^{I}\neq n_{x}^{J}\\}$. A more
compact representation of CFGs is to use two bits to store the occupancy
number of each spatial orbital. Specifically, $(00)_{2}$, $(01)_{2}$ and
$(11)_{2}$ are used for $n_{i}=0$, 1, and 2, respectively. A CFG can then be
represented by an array, OrbOccBinary, of 64-bit integers. This representation
was used for both bra and ket CFGs in our previous implementation. Now it is
discovered that the single occupancy in the bra CFG can be represented by
$(10)_{2}$ instead of $(01)_{2}$, so as to make the identification of ROT_Code
more efficiently. As can be seen from Table 4, the bitwise XOR operations on
the so-defined binary occupation codes (BOC) always give $(00)_{2}$ for common
zero or doubly occupied orbitals, such that the identification of ROT_Code can
be achieved based fully on bit operations. Such double representation requires
$M_{CFG}=2\times 8\left(\lfloor\frac{N_{orb}-1}{32}\rfloor+1\right)$ bytes of
memory for storing a CFG.

Table 4: Bitwise XOR operations on binary occupation codes (BOC) Bra Occ | BOC | Ket Occ | BOC | XOR
---|---|---|---|---
0 | $(00)_{2}$ | 0 | $(00)_{2}$ | $(00)_{2}$
2 | $(11)_{2}$ | 2 | $(11)_{2}$ | $(00)_{2}$
1 | $(10)_{2}$ | 1 | $(01)_{2}$ | $(11)_{2}$
0 | $(00)_{2}$ | 1 | $(01)_{2}$ | $(01)_{2}$
0 | $(00)_{2}$ | 2 | $(11)_{2}$ | $(11)_{2}$
1 | $(10)_{2}$ | 0 | $(00)_{2}$ | $(10)_{2}$
1 | $(10)_{2}$ | 2 | $(11)_{2}$ | $(01)_{2}$
2 | $(11)_{2}$ | 0 | $(00)_{2}$ | $(11)_{2}$
2 | $(11)_{2}$ | 1 | $(01)_{2}$ | $(10)_{2}$

A CFG $|I\rangle$ can generate $N^{I}_{S,S}$, genealogically coupled CSFs
$\\{|I\mu\rangle\\}_{\mu=0}^{N_{S,S}^{I}-1}$ of total spin $S$, with
$N^{I}_{S,S}$ being

$N_{S,M=S}^{I}=\frac{2S+1}{\frac{1}{2}N_{o}^{I}+S+1}C_{N_{o}^{I}}^{\frac{1}{2}N_{o}^{I}-S},$
(23)

where $N_{o}^{I}$ is the (seniority) number of open-shell orbitals in CFG
$|I\rangle$. Such CSFs are characterized uniquely by the Shavitt step number89
sequences $\\{(d_{0}^{\mu}d_{1}^{\mu}\cdots)\\}|_{\mu=0}^{N_{S,S}^{I-1}}$,
which can be arranged in a lexical order (i.e., $|I\mu\rangle$ precedes
$|I\nu\rangle$ if and only if $d_{p}^{\mu}<d_{p}^{\nu}$ for
$p=\min\\{x|d_{x}^{\mu}\neq d_{x}^{\nu}\\}$), so as to fix the relative
ordering $\mu$ uniquely. A 32-bit integer is required to store the relative
ordering of CSFs $\\{|I\mu\rangle\\}$. Since both its BOC and lexical order
have to be stored, the memory requirement for storing a CSF is $M_{CFG}+8$
bytes due to memory alignment. The odering of CSFs is defined as follows: if
CSF $|I\mu\rangle$ and $|J\nu\rangle$ stem from the same CFG, then
$|I\mu\rangle<|J\nu\rangle$ if $\mu<\nu$; otherwise, let $Occ_{Y}$ ($Y=I,J$)
be the corresponding occupation array of CFG $|Y\rangle$ and
$i=\arg\max_{j}\\{Occ_{I}[j]\neq Occ_{J}[j]\\}$, then
$|I\mu\rangle<|J\nu\rangle$ if $Occ_{I}[i]>Occ_{J}[i]$. In view of Eq. (23), a
large number of CSFs can be generated from a CFG of high seniority $N_{o}^{I}$
but only a few of them may contribute discernibly. Therefore, it is absolutely
necessary to do individual selections of CSFs.

After the above detailed description of TUGA, a gross comparison with the
Table-CI approach79, 80, 81, 82 is in order. Although both are configuration-
driven algorithms for calculating the Hamiltonian matrix elements over
randomly selected CSFs, they actually differ significantly. First of all,
Table-CI is closely related to the symmetric group approach (SGA)90 in view of
its explicit tabulation of line-up permutations, whereas TUGA is plainly an
implementation of the UGA 86, 87, 88, in terms of the tabulated ROTs rather
than the famous graphical representation89, 91 (which is efficient only for
well-structured wave functions). Secondly, TUGA and Table-CI differ in the
handling of open-shell orbitals. While Table-CI always defines CFGs with
doubly occupied orbitals sitting together and open-shell orbitals sitting
together (which involves reordering of the orbitals for each CFG), TUGA adopts
a universal orbital ordering defined in the very beginning. The common open-
shell orbitals are eliminated in Table-CI by virtue of line-up permutations,
leaving a relatively small number of unique interaction patterns between SDs
resulting from non-identical open-shell orbitals. However, the BCCs between
CSFs are dependent on the relative positions of common open-shell orbitals.
Therefore, they are not permuted in TUGA but stay where they are after
deleting the common zero or doubly occupied orbitals. This increases merely
the lengths of code sequences but not the number of ROTs (NB: every CFG pair
corresponds to a ROT). In essence, the ROTs in TUGA play the same role in
classifying unique interacting CFG pairs at the integral level as the Tables
defined in Table-CI. The former are calculated on an as-needed basis, whereas
the latter, characterized by a number of parameters (i.e., $\Delta K$, $P$,
$Q$, and $R$ determined by the relative occupations of CFG pairs; see Ref. 82
for more details), are created and stored in advance. Thirdly, TUGA evaluates
the Hamiltonian matrix elements directly between CSFs, whereas Table-CI has to
assemble them via transformations of those between SDs (which can be viewed as
an effective means to avoid the overly complicated matrix elements of line-up
permutations $Q$ over CSFs, as required by the genuine SGA90).

In short, TUGA has been based on the single, simple fact that the segment
values for common zero or doubly occupied orbitals are always one under the YK
phase87, such that the BCCs are determined solely by those common open-shell
orbitals or orbitals with different occupancy numbers (which are encoded in
the code sequences) and can be calculated directly through products of their
segment values.

##  4 Selection of CSFs

The aim of selection is to find a better variational space $P$ by feeding in a
guess space $P_{0}$. This can be achieved in two steps, ranking and pruning.
In the ranking step, proper rank values for the CSFs in the FOIS $Q$
$(=1-P_{0})$ will be evaluated. Those CSFs with rank values larger than the
given ranking-threshold are put into $P_{0}$, leading to an expanded space
$P_{1}$. In the pruning step, the Hamiltonian matrix in $P_{1}$ is first
constructed and diagonalized. Those CSFs with coefficients smaller in absolute
value than the given pruning-threshold are then discarded, so as to reduce
$P_{1}$ to $P$. The procedure is repeated until the energies obtained by
diagonalizing $P$ and $P_{0}$ are sufficiently close. Such a two-step
selection scheme has been adopted by most sCIs. Yet, a number of refinements
has been introduced in our own implementation66:

  1. (1)

A combined integral- and coefficient-driven algorithm is used for the
selection, viz. prior to the selection of individual CSFs with their
(approximate) first-order coefficients as rank values, a integral-driven
screening of unimportant, doubly connected CFG pairs is first carried out
based on upper bounds $\tilde{H}^{IJ}$ of $\langle I\mu|H_{2}^{2}|J\nu\rangle$
(which are independent of CSFs, see Table 3).

  2. (2)

To minimize the usage of memory, an increment access of the $Q$ space is
invoked, with a dynamically adjusted integral-threshold $\varepsilon_{k}$
($=\frac{1}{2}\varepsilon_{k-1})$. For a given $\varepsilon_{k}$ of iteration
$k$ (i.e., $Q[\varepsilon_{k}]$), the selection is repeated until
$P[\varepsilon_{k}]$ and $P_{0}[\varepsilon_{k}]$ are sufficiently similar in
composition (i.e., $\frac{|P_{0}[\varepsilon_{k}]\bigcap
P[\varepsilon_{k}]|}{|P_{0}[\varepsilon_{k}]\bigcup P[\varepsilon_{k}]|}\geq
0.95$).

  3. (3)

Natural orbitals (NO) are also generated dynamically (i.e., for each converged
$P[\varepsilon_{k}]$; for more details, see Algorithm 4 in Ref. 66).

  4. (4)

A residue-based algorithm92, 62, 63 is used for establishing the connections
between randomly selected CFGs and CSFs.

  5. (5)

The diagonalization of Hamiltonian matrix is performed with the iVI-Buf(0+3)
approach83, 85, which has the potential to access directly the roots of a
given energy window.

  6. (6)

The size of the resulting variational space $P[\varepsilon_{k}=0]$ is
controlled by a single parameter $C_{min}$ (coefficient pruning-threshold),
for other parameters have been fixed to conservative values.

It has been demonstrated that the above selection can generate a very compact
yet accurate variational space. However, its efficiency is not yet
satisfactory because of the following reasons related to the first four points
above: (a) the adopted ranking criterion is too good (and hence expensive) for
the purpose of selection. (b) Although the incremental access of the $Q$ space
is very effective in reducing memory requirement and does not affect the final
result, it does result in additional iterations and hence slow down the
selection procedure. (c) The generation of NOs and hence transformation of
integrals are performed too many times, which is unnecessary. (d) The residue-
based algorithm is not optimal for identifying the connections between the
CFGs in $P$. In particular, it does not allow for an easy construction of the
Hamiltonian matrix in CSR format. To resolve these issues, we here introduce a
simpler ranking criterion, remove the incremental access of $Q$, simplify the
generation of NO, and replace the residue-based algorithm with a particle-hole
algorithm for Hamiltonian construction (see Sec. 5).

In view of many-body perturbation theory (MBPT), the most rigorous ranking of
a CSF $|I\mu\rangle\in Q$ is its first-order coefficient in absolute value

$f(|I\mu\rangle,C_{min})=\left|\frac{\sum_{|J\nu\rangle\in
P_{0}}H^{IJ}_{\mu\nu}C^{J}_{\nu}}{E_{0}-H_{\mu\mu}^{II}}\right|\geq
C_{min},\quad H^{IJ}_{\mu\nu}=\langle I\mu|H|J\nu\rangle,$ (24)

which is to be called ‘CIPSI criterion’, with CIPSI standing for
‘configuration interaction with perturbative selection made iteratively’73. If
those $|J\nu\rangle\in P_{0}$ with small coefficients are excluded from the
sum in the numerator of Eq. (24), it may be termed ‘pruned CIPSI criterion’ or
simply ‘ASCI criterion’, as adopted by adaptive sampling CI60. Noticing that
the computational expense of Eq. (24) results from both the sum in the
numerator and the evaluation of the denominator, a dramatically simplified
criterion is considered in heat-bath CI (HBCI)50, viz.,

$f(|I{\mu}\rangle,|J\rangle,\epsilon_{Q})=\max_{\nu}|H^{IJ}_{\mu\nu}C^{J}_{\nu}|\stackrel{{\scriptstyle\mathrm{CSF}\rightarrow\mathrm{SD}}}{{\longrightarrow}}\max_{\nu}|H_{\mu\nu}C_{\nu}|\geq\epsilon_{Q},$
(25)

in the basis of SDs though. This way, those determinants $\\{|\mu\rangle\\}$
that are doubly excited from $|\nu\rangle\in P_{0}$ are never accessed if
$|H_{\mu\nu}|<\epsilon_{Q}/|C_{\nu}|$93. However, such integral-driven
selection leads usually to a variational space that is much less compact than
that by the coefficient-driven selection (24), particularly when a large basis
set is used. A possible remedy to this is to introduce94 an approximate
denominator to condition (25), without sacrificing the efficiency. Our
previous ‘iCI criterion’66 amounts to combining _dynamically_ these integral-
and coefficient-driven algorithms for selecting doubly excited CSFs, viz.,

$f(|I\mu\rangle,\varepsilon_{k},C_{min})=\left|\frac{\sum_{J\in
P_{0}^{(k)}}^{(\varepsilon_{k})}(\sum_{\nu,|J\nu\rangle\in
P_{0}^{(k)}}H^{IJ}_{\mu\nu}C^{J(k)}_{\nu})}{E^{(k)}-H_{\mu\mu}^{II}}\right|\geq
C_{min},$ (26)

where $\varepsilon_{k}$ is a dynamically adjusted integral-threshold and the
summation over CFG $|J\rangle\in P_{0}^{(k)}$ is subject to the following
condition

$\max_{\nu}|\tilde{H}^{IJ}C_{\nu}^{J(k)}|\geq\varepsilon_{k}=\frac{1}{2}\varepsilon_{k-1},\quad|J\rangle\in
P_{0}^{(k)}.$ (27)

The original motivation for this particular criterion is to reduce memory
usage as much as possible and meanwhile to end up with a variational space
that is as compact as possible. It appears that both are not really necessary:
the memory requirement in the PT2 step is much larger and a slight sacrifice
of the compactness can be overcompensated by the dramatic gain in efficiency.
Therefore, we introduce here a new ‘iCI criterion’, that is, for a given
$|I\mu\rangle\in Q=1-P_{0}$, if there exists a $|J\rangle\in P_{0}$ for which
$f(|I\mu\rangle,|J\rangle,C_{min})$ is true, then $|I\mu\rangle$ is selected.
The boolean function $f(|I\mu\rangle,|J\rangle,C_{min})$ is defined as

  1. (A)

If $|I\rangle$ is identical with or singly excited from $|J\rangle$, then

$f(|I\mu\rangle,|J\rangle,C_{min})=\left(\max_{\nu}(|H^{IJ}_{\mu\nu}C_{J\nu}|)\geq
C_{min}\right)\text{ and
}\left(\max_{\nu}\left(\left|\frac{H^{IJ}_{\mu\nu}C_{J\nu}}{E_{0}-H^{II}_{\mu\mu}}\right|\right)\geq
C_{min}\right);$ (28)

  2. (B)

If $|I\rangle$ is doubly excited from $|J\rangle$, then

$\displaystyle f(|I\mu\rangle,|J\rangle,C_{min})$ $\displaystyle=$
$\displaystyle\left(\max_{\nu}(|\tilde{H}^{IJ}C_{J\nu}|)\geq
C_{min}\right)\text{ and }\left(\max_{\nu}(|H^{IJ}_{\mu\nu}C_{J\nu}|)\geq
C_{min}\right)$ (29) $\displaystyle\text{ and
}\left(\max_{\nu}\left(\left|\frac{H^{IJ}_{\mu\nu}C_{J\nu}}{E_{0}-H^{II}_{\mu\mu}}\right|\right)\geq
C_{min}\right).$

Literally, for case (A), loop over all $|I\mu\rangle$ in $Q$ and evaluate
$H^{IJ}_{\mu\nu}$ for all CSFs $|J\nu\rangle\in P_{0}$. If
$\max_{\nu}(|H^{IJ}_{\mu\nu}C_{J\nu}|)$ is larger than $C_{min}$ then evaluate
$H^{II}_{\mu\mu}$; otherwise discard $|I\mu\rangle$. If
$\max_{\nu}(|\frac{H^{IJ}_{\mu\nu}C_{J\nu}}{E_{0}-H^{II}_{\mu\mu}}|)$ is
larger than $C_{min}$ then $|I\mu\rangle$ is selected. As for case (B), only
those doubly excited CFGs $|I\rangle$ with $\tilde{H}^{IJ}$ larger than
$C_{min}/\max_{\nu}|C_{J\nu}|$ need to be generated (i.e., those unimportant
ones are never touched, as illustrated in Fig. 1). For such $\\{|I\rangle\\}$,
the remaining step is the same as case (A).

Figure 1: Screening of doubly excited configurations based on upper bounds.

It should be clear that, for the same $C_{min}$, the revised ‘iCI criterion’
(28)/(29) will yield a smaller number of CSFs than the original one (26)/(27),
for the latter may accumulate those CSFs of small coefficients until the
condition is fulfilled (cf. the second sum in Eq. (26)), whereas the former
just picks up the CSF with the largest coefficient. Note in passing that, to
be compared directly with this ‘iCI criterion’, the ‘HBCI criterion’ (25)
should be rewritten as

  1. (a)

If $|I\rangle$ is identical with or singly excited from $|J\rangle$, then

$f(|I\mu\rangle,|J\rangle,C_{min})=\max_{\nu}(|H^{IJ}_{\mu\nu}C_{J\nu}|)\geq
C_{min};$ (30)

  2. (b)

If $|I\rangle$ is doubly excited from $|J\rangle$, then

$f(|I\mu\rangle,|J\rangle,C_{min})=\left(\max_{\nu}(|\tilde{H}^{IJ}C_{J\nu}|)\geq
C_{min}\right)\text{ and }\left(\max_{\nu}(|H^{IJ}_{\mu\nu}C_{J\nu}|)\geq
C_{min}\right),$ (31)

which amounts to skipping the final step in Eqs. (28)/(29).

The overall selection procedure is illustrated in Fig. 2. Four additional
points still deserve to be mentioned:

  1. (1)

For multiple states, the selection will be carried out in a state-collective
way. That is, the boolean function $f_{k}$ defined in Eq. (24), (28)/(29) or
(30)/(31) is examined for each state $k$ and $|I\mu\rangle$ is selected if
$\exists k$ and $|J\nu\rangle$ s.t. $f_{k}=\mathrm{true}$.

  2. (2)

Duplicates in the selected CSFs $\\{|I\mu\rangle\\}\in Q$ are removed as
follows. The for-loop over CFGs $|J\rangle$ in $P_{0}$ is distributed among
working threads. For each thread, once the memory usage has reached a given
size (e.g., 2 GB), the selected CSFs $\\{|I\mu\rangle\\}$ are sorted
immediately, with duplicates removed by a single loop. To expedite this, the
CSFs are classified according to their seniorities and highest occupied
orbitals. The surviving CSFs $\\{|I\mu\rangle\\}$ are then compared with those
in $P_{0}$ to further remove the duplicates. Once the loop over $|J\rangle$ is
done, the surviving CSFs $\\{|I\mu\rangle\\}$ from each thread are merged to
the main thread and a final removal of duplicates is then executed.

  3. (3)

The selection is terminated once the measure of composition similarity,
$\frac{|P_{0}\bigcap P|}{|P_{0}\bigcup P|}$, exceeds the threshold $S_{P}$
(e.g., 0.95).

  4. (4)

If desired, NOs can be generated in a partially dynamic manner, e.g., going
from $C^{\prime}_{min}=4C_{min}$ to $2C_{min}$ and finally to $C_{min}$. For
each $C_{min}^{\prime}$, the above selection is performed to find the
corresponding variation space $P$, from which the NOs can be generated. After
integral transformation, the Hamiltonian matrix is reconstructed and
diagonalized, with unimportant CSFs pruned. This extra step is important for
the expansion coefficients of the wave function are affected by orbital
rotations. The NOs generated this way usually give rise to a more compact
variational space than those generated directly at the target $C_{min}$.

Figure 2: Flowchart for iCIPT2.

##  5 Connections between selected CSFs

Having identified the most important CSFs spanning space $P$, their
connections must be established for constructing and updating the Hamiltonian
matrix. In the previous implementation of iCIPT266, we adopted the residue
array-based algorithm92, 62, 63. The $n$th-order residues of an
$N_{e}$-electron CFG $|I\rangle$ are those $(N_{e}-n)$-electron CFGs that can
be generated by removing $n$ electrons from $|I\rangle$ in all possible ways.
A nice point of residues lies in that they provide connection information on
the parent CFGs. For instance, different CFG pairs sharing the same first-
order residue (FOR) are identical or singly connected, while those sharing the
same second-order residue (SOR) are identical, singly or doubly connected.
However, for the purpose of constructing and updating the Hamiltonian matrix
in the $P$ space, the residues have some drawbacks: (1) the number of SORs
scale linearly with respect to the number $N_{cfg}$ of CFGs in $P$ but scale
quadratically with respect to the number $N_{e}$ of correlated electrons, such
that both the consumption of memory and the repeated update of the SORs (which
is necessary during the iterative selection) are unfavorable when $N_{e}$ is
large. (2) There may exist numerous redundant residues that do not provide any
connection information, especially in the presence of many doubly occupied
orbitals. (3) More seriously, the CFGs connected to a given CFG may be
scattered widely in location, such that it is not easy to construct the
Hamiltonian matrix in CSR format for efficient matrix-vector multiplications
(unless a sorting step is taken but which itself is not cheap). It turns out
that all these issues can be resolved by going to a particle-hole
representation of the CFGs in $P$. The holes and particles are just the doubly
occupied and the remaining singly or zero occupied orbitals in the common
reference CFG $|R\rangle$ (which is already adopted in Eq. (17)),
respectively. Then, the array OrbOccBinary representing a CFG can be
decomposed into a hole (HoleStr) and a particle part (PartStr), such that
connections between CFG pairs can be identified from those between HoleStr’s
and between PartStr’s. To this end, we first define single transitions within
the hole or particle orbitals as single inner excitations (SIE), whereas those
between the two sets as single outer excitations (SOE) or de-excitations
(SDE). While SIEs do not change the lengths of HoleStr and PartStr, a SOE
(SDE) will increase (decrease) their lengths by one. Two HoleStr’s (PartStr’s)
are connected only if their lengths differ less than three. A given
HoleStr/PartStr StrA with length len can be connected with a HoleStr/PartStr
StrB in eight possible ways,

  1. 1.

StrB with length len+2 and two SOEs;

  2. 2.

StrB with length len+1 and one SOE;

  3. 3.

StrB with length len+1, one SIE and SOE;

  4. 4.

StrB with length len and one SIE;

  5. 5.

StrB with length len and two SIE;

  6. 6.

StrB with length len-1 and one SDE;

  7. 7.

StrB with length len-1, one SIE and one SDE;

  8. 8.

StrB with length len-2 and two SDEs,

in addition to the case (denoted as 0) when StrB is identical with StrA.

The data structure of HoleStr is very simple. It consists of three parts: (1)
an array of 64-bit integers recording the corresponding HoleStr; (2) eight
hole connection arrays $\mathcal{C}_{ph}$ recording the above connections to
other HoleStr’s; (3) one CFG array $\mathcal{G}_{ph}$ recording the indices of
CFGs sharing the same HoleStr. The data structure of PartStr is exactly the
same. The connections between CFG pairs can then be identified in a simple
way, since singly connected CFG pairs can only be one of the three single-
connection types (SCT) $\\{(2,2),(4,0),(6,6)\\}$, whereas doubly connected CFG
pairs can only be one of the eight double-connection types (DCT)
$\\{(1,1),(2,3),(3,2),(4,4),(5,0),(6,7),(7,6),(8,8)\\}$. For instance, for a
given CFG $|I\rangle$ with HoleStr HoleI and PartStr PartI, one can first find
all HoleJ’s and PartJ’s according to the SCTs. Then, for each pair of HoleJ
and PartJ, the intersection of their $\mathcal{G}_{ph}$’s contains all the
CFGs that are singly connected to CFG $|I\rangle$. The doubly connected CFG
pairs can be identified in the same way.

The Hamiltonian matrix can be constructed as follows: (A) loop over all CFGs
$|I\rangle\in P$. (B) For each $|I\rangle$, find all CFGs $|J\rangle\in P$
connected with $|I\rangle$ but with indices smaller than that of $|I\rangle$.
(C) Loop over $|J\rangle$ and calculate the Hamiltonian matrix elements
between $|I\mu\rangle$ and $|J\nu\rangle$, which can readily be compressed in
CSR format because the indices are now contiguous. As for the update of the
Hamiltonian matrix, suppose that the old and new CSF spaces are $P_{1}\subset
P_{2}$ and $P_{2}$, respectively, and the CSFs in $P_{2}-P_{1}$ have been
grouped according to their CFGs. After updating the connections between
HoleStr’s and between PartStr’s, just loop over all newly added CFGs
$|I\rangle$ in $P_{2}-P_{1}$ and do (B) and (C) above.

In summary, as far as the handling of the $P$ space is concerned, the present
particle-hole algorithm is advantageous over the residue-based algorithm92,
62, 63 in several aspects, especially when $N_{e}$ is large. First of all, the
metadata of the HoleStr’s and PartStr’s is much smaller in size than that of
the FORs and SORs, which is reflected by the much reduced memory consumption.
For instance, for a variational space with 1164710 unique CFGs or 3059395 CSFs
in the all-electron calculation of Cr2 with the Ahlrichs SV basis set95, the
memory cost of the HoleStr’s and PartStr’s is only 0.55 GB, whereas those of
the FORs and SORs are 2.48 GB and 29.1 GB, respectively. Secondly, the simple
structure of the HoleStr’s and PartStr’s allows for an easy construction of
the Hamiltonian matrix in CSR format. For the case of Cr2, the Hamiltonian
construction is speeded up by $2.5\times$. Thirdly, the connections between
HoleStr’s and between PartStr’s can readily be reutilized. For example, when
updating the Hamiltonian matrix from $P_{1}$ to $P_{2}$ mentioned previously,
if some CSFs of a CFG in $P_{2}-P_{1}$ are present in $P_{1}$, the connections
of this CFG to other CFGs in $P_{1}$ have already been encoded in the
connections between HoleStrs and between PartStrs of $P_{1}$. More
interestingly, such connections can also be shared by CSF spaces of different
spatial and/or spin symmetries, which facilitates the simultaneous calculation
of several states of different spatial and/or spin symmetries with a common
set of orthonormal orbitals.

As a final note, it deserves to be mentioned that, because of the random
nature of the $P$ space, the number of CFG pairs therein sharing the same
Rot_Code is usually very small, such that the reutilization of BCCs is
ineffective: it cannot compensate the overhead necessary for sorting the CFG
pairs (which requires both substantial memory and synchronous operations).
Instead, recalculating the BCCs whenever needed turns out to be more
efficient. However, the situation is different for the PT2 step, where the
reutilization of BCCs by TUGA is essential.

##  6 Constraint-Based ENPT2

As long as the selectively determined variational space $P$ and hence the
zeroth-order eigenpairs $\\{E_{k}^{(0)},|\Psi_{k}^{(0)}\rangle\\}$ are good
enough, the remaining dynamic correlation can be accounted for accurately at
the lowest level of theory. Here, we adopt the state-specific ENPT2

$\displaystyle E_{c,k}^{(2)}$ $\displaystyle=$
$\displaystyle\sum_{|I\mu\rangle\in Q}\frac{|\langle
I\mu|H|\Psi^{(0)}_{k}\rangle|^{2}}{E^{(0)}_{k}-H^{II}_{\mu\mu}}$ (32)
$\displaystyle=$ $\displaystyle\sum_{|I\mu\rangle\in
Q}\frac{\left|\sum_{|J\nu\rangle\in
P}H^{IJ}_{\mu\nu}C^{J}_{\nu,k}\right|^{2}}{E_{k}^{(0)}-H_{\mu\mu}^{II}}.$ (33)

Before evaluating this energy expression, several important points should be
observed. Firstly, what the interaction really does is to excite one or two
electrons from space $P$ to space $Q$ in all possible ways, leaving $N_{e}-1$
or $N_{e}-2$ CFGs in $P$, which are nothing but the respective FORs and SORs
mentioned in the previous section. Reversely, the excited CFGs can be
generated by inserting one and two electrons into the FORs and SORs,
respectively. Therefore, the FORs and SORs provide natural connections between
the two spaces and are hence the proper organizing units. Secondly, the number
of CFG pairs sharing the same Rot_Code can be very large, such that is
essential to reutilize the BCCs by TUGA. Thirdly, the summation over
$|I\mu\rangle\in Q$ can be extended to include also the CSFs in $P$, by
precomputing and finally subtracting their PT2-like energy
$\bar{E}_{c,k}^{(2)}$63, viz.,

$\displaystyle E_{c,k}^{(2)}$
$\displaystyle=\tilde{E}_{c,k}^{(2)}-\bar{E}_{c,k}^{(2)},$ (34)
$\displaystyle\tilde{E}_{c,k}^{(2)}$ $\displaystyle=\sum_{|I\mu\rangle\in
W}\frac{\left|\sum_{|J\nu\rangle\in
P,|J\nu\rangle\neq|I\mu\rangle}H_{\mu\nu}^{IJ}C_{\nu,k}^{J}\right|^{2}}{E_{k}^{(0)}-H_{\mu\mu}^{II}},\quad
W=P\cup Q,$ (35) $\displaystyle\bar{E}_{c,k}^{(2)}$
$\displaystyle=\sum_{|I\mu\rangle\in P}\frac{\left|\langle
I\mu|H|\Psi_{k}^{(0)}\rangle-
C_{\mu,k}^{I}H^{II}_{\mu\mu}\right|^{2}}{E_{k}^{(0)}-H_{\mu\mu}^{II}}$ (36)
$\displaystyle=\sum_{|I\mu\rangle\in
P}(C_{\mu,i}^{I})^{2}(E_{i}^{(0)}-H^{II}_{\mu\mu}),$ (37)

where use of the relation $\langle
I\mu|H|\Psi_{k}^{(0)}\rangle=C_{\mu,k}^{I}E_{k}^{(0)}$ has been made when
going from Eq. (36) to Eq. (37). The negative term in the numerator of Eq.
(36) arises from the fact that the diagonal terms have been excluded in Eq.
(35). This way, there is no need to double check whether the excited CSFs
belong to $Q$ or $P$, which is every expensive. Fourthly, since there is no
coupling between the $|I\mu\rangle$ functions, the whole $W$ space can be
decomposed into a series of disjoint subspaces $\\{W_{i}\\}_{i}^{N_{s}}$, so
as to reduce memory requirement and meanwhile facilitate the parallelization.

Following the idea of constraint PT263, a CFG subspace $W_{i}$ of $W$ can be
defined by a constraint consisting of $L_{c}$ highest occupied orbitals (say,
$p_{1},p_{2},\cdots,p_{L_{c}}$ in ascending order) as well as their occupation
numbers (say, $n_{p_{1}},n_{p_{2}},\cdots,n_{p_{L_{c}}}$). The number $N_{s}$
of such subspaces is bounded by $C_{N_{orb}}^{L_{c}}\times 2^{L_{c}}$. To
achieve this, the FORs and SORs of $P$ are first generated and sorted, with
the unique ones recorded in arrays $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$,
respectively. Here, each unique residue is associated with an array of records
{indx,orbj,orbl} to trace how it arises: indx records the CFG in $P$, whereas
orbj and orbl record the orbitals from which electrons are removed (i.e., the
$j$ and $l$ indices of $e_{ij,kl}$). Note that orbj can be equal to orbl for a
SOR, whereas only orbsj is needed for a FOR. Then, those FORs and SORs that
can generate the CFGs belonging to $W_{i}$ must be identified. Since the
occupation pattern of a CFG in $W_{i}$ must be
$(0,\cdots,0,n_{p_{L_{c}}},0,\cdots,0,n_{p_{L_{c}-1}},0,\cdots,0,n_{p_{1}})$
for orbitals of indices not smaller than $p_{1}$, a valid residue with
occupation numbers $\\{m_{n}\\}_{n=0}^{N_{orb-1}}$ must be subject to the
restrictions $m_{p}\leq n_{p}$ for $p\geq p_{1}$ and
$R_{t}=\sum_{p=p_{1}}^{N_{orb}-1}(n_{p}-m_{p})\in[0,2]$. For the case of FORs,
$R_{t}$ can be either 0 or 1, which are labeled c1 and c2, respectively. For
the former, $n_{p}=m_{p}$ for all $p\geq p_{1}$. For the latter, there exists
one orbital $q$ such that $n_{q}=m_{q}+1$ for $q\in[p_{1},p_{L_{c}}]$ and
$n_{p}=m_{p}$ for $p\geq p_{1}$ other than $q$. For the case of SORs, $R_{t}$
can also be 2, which leads to two additional conditions: (a) c3: if there
exists $q\geq p_{1}$ with $n_{q}=2$ and $m_{q}=0$, then $n_{p}=m_{p}$ for
$p\geq p_{1}$ other than $q$. (b) c4: if there exist $q$ and $r$ with
$n_{q}-m_{q}=1$ and $n_{r}-m_{r}=1$, then $n_{p}=m_{p}$ for $p\geq p_{1}$
other than $q$ and $r$. The valid FORs and SORs can be classified according to
their occupation patterns for orbitals $p\geq p_{1}$, such that those of the
same occupation pattern are located contiguously in
$\mathcal{R}_{1}$/$\mathcal{R}_{2}$, thereby forming a segment. The numbers of
such patterns/segments are $1+C_{L_{c}}^{1}$ for FORs (1 for c1 and
$C_{L_{c}}^{1}$ for c2) and $1+C_{L_{c}}^{1}+C_{L_{c}}^{1}+C_{L_{c}}^{2}$ for
SORs ($C_{L_{c}}^{1}$ for c3 and $C_{L_{c}}^{2}$ for c4). The head $R_{min}$
(tail $R_{max}$) of a segment is determined by distributing the remaining
$N_{e}-\sum_{p\geq p_{1}}^{N_{orb}-1}m_{p}-N_{R}$ electrons ($N_{R}=1$ for
FORs and 2 for SORs) to orbitals as close to (far away from) $p_{1}$ as
possible. They can be used to locate the relevant residues by binary search of
the segments. For a given residue, each way of adding electrons will generate
a CFG $|I\rangle\in W_{i}$, whose connections with the CFGs $\\{|J\rangle\\}$
in $P$ can be established by looping over the records {indx,orbj,orbl}
associated with the residue (NB: each indx corresponds to a CFG $|J\rangle$).
Such pair connection records (PCR) are stored in array $\mathcal{C}_{IJ}$,
which is to be sorted to group those records sharing the same $|I\rangle$
together. Note in passing that in principle only the SORs are needed to
generate such PCRs. However, very many duplicate singly connected CFG pairs
can be generated from the SORs and the removal of them can be very costly.
Instead, it is more favorable to take the (non-duplicate) singly connected CFG
pairs from the FORs and simply dump the ones from the SORs. The price to pay
is just some extra memory. For a given PCR, an interaction record
$\\{ROT\\_Code,I,J,integrals\\}$ is further needed to calculate
$\sum_{\nu}\langle I\mu|H|J\nu\rangle C_{\nu,k}^{J}$ for all CSFs of
$|I\rangle\in W_{i}$. It can be constructed by first identifying the ROT for
$(I,J)$ and then fetching the integrals based on ROT_Orb. All the interaction
records are stored in array $\mathcal{V}_{IJ}$, which is to be sorted
according to ROT_Code to reuse the BCCs by TUGA.

At this moment two key factors should be observed: (1) the sizes of different
subspaces $W_{i}$ may differ by orders of magnitude, which renders the load in
parallelization extremely imbalanced on one hand and the reutilization of BCCs
of very small subspaces virtually impossible. (2) Both the generation and
sorting of the connection and interac- tion data can be very expensive. The
former can be resolved by precomputing the sizes of subspaces, so as to merge
small subspaces together or further split very large subspaces into smaller
ones. Although not very cheap, the expense of this null step is
overcompensated by the gain in parallelization efficiency. As for the latter,
the red-black tree-based algorithm adopted in our previous implementation66
turns out to be inefficient: too much time is wasted on the searching and
insertion, due to the fact that the nodes in red-black trees are scattered in
memory, thereby resulting in a very low cache-hit rate. Since array appending
and sorting are more efficient than tree insertion and searching, we here
employ arrays instead of trees for the data. What is essential is then how to
impose a particular structure to the array of raw data, according to which a
possibly optimal sorting algorithm can be found. For instance, if the array
occurs automatically as sorted subarrays with an equal length, the Timsort
algorithm77, 78 would be the best choice (NB: Timsort appears to work well
even if some subarrays have different lengths). To this end, the above
connection array $\mathcal{C}_{IJ}$ is to be generated as follows: (1) for the
c2 segment of FORs or the c3 and c4 segments of SORs, the electron-accepting
orbitals are already fixed and are the same for all the residues within the
segment. Since the FORs and SORs are already sorted, the PCRs generated from
such segment are automatically sorted. (2) For the c1 segment of FORs or the
c2 segment of SORs, one electron should be added to orbital $s\in[0,p_{1}-1]$.
If the loop over $s$ proceeds in descending order, the generated CFGs will be
in ascending order, leading to ordered PCRs for each residue. Therefore, the
PCRs for such segment consist of several ordered subarrays, thereby well
suited for the Timsort algorithm. (3) For the c1 segment of SORs, two
electrons should be added to orbitals $i,k\in[0,p_{1}-1]$ with $k\geq i$ (NB:
$i,k$ are the indices of $e_{ij,kl}$). If the outer loop over $k$ and the
inner loop over $i$ for each $k$ both proceed in descending order, the PCRs
will be of the same structure as case (2) and hence fit the Timsort algorithm
as well. After looping over all residue segments, the whole $\mathcal{C}_{IJ}$
array will consist of several sorted subarrays and can finally be sorted by
Timsort. In contrast, the interaction array $\mathcal{V}_{IJ}$ cannot be
prepared in any particular structure and is hence sorted simply by std::sort()
in C++.

The efficacy of the above PT2 algorithm can quickly be revealed by again
taking Cr2/SV as an example. It turns out that the smaller the $C_{min}$
(i.e., the larger the variational space $P$), the larger the speedup of PT2
over the previous, red-black tree-based implementation66, amounting to
21$\times$ for $P$ with 1.1$M$ CSFs. If desired, two additional cutoffs can be
employed for the PT2 correction: (1) those doubly excited CFGs $|I\rangle$
with $\tilde{H}^{IJ}\max_{\nu}|C^{J}_{\nu,k}|$ smaller than a threshold (e.g.,
$10^{-3}C_{min}$) can be neglected. (2) Those CSFs $|I\mu\rangle$ with
$|\langle I\mu|H|\Psi_{k}^{(0)}\rangle|$ smaller than a threshold (e.g.,
$5\times 10^{-9}$) can further be neglected.

##  7 Results and discussion

All the calculations were performed with the BDF program package96, 97, 98,
99, 100 on one node with 4 Intel(R) Xeon(R) Gold 6240 CPUs (in total 72
physical cores) and 768 GB memory.

###  7.1 Comparison of different ranking criteria

The very first point to be checked is the compactness of the final variational
space $P$ determined by different ranking criteria discussed in Sec. 4. To
this end, all-electron nonrelativistic calculations were performed on the
ground state of Cr2 with the Ahlrichs SV basis set95 at two interatomic
distances, $R_{eq}$ (=1.68 Å) and $2.0R_{eq}$. It can be seen from Fig. 3
that, as the most rigorous ranking criterion, CIPSI (24) does lead to the
smallest number of CSFs and hence the most compact wave function for the same
variational energy (NB: although not documented here, the original iCI
criterion (26)/(27) is virtually identical with CIPSI (24)). It appears that
HBCI (30)/(31) follows closely CIPSI. However, as can be seen from Fig. 4,
this arises with a high price: the variational space $P_{1}$ determined by
HBCI (30)/(31) is actually very large and is only reduced (by a factor of more
than 10!) by the pruning step. That is, the HBCI ranking criterion is too
loose, thereby bringing in too many unimportant CSFs, more than 90% of which
are pruned away. This is clearly a waste of time for Hamiltonian construction
and diagonalization. Note in passing that the HBCI ranking used here is
different from the original one50, which does not invoke any pruning. It is
also clear that the iCI criterion (28)/(29) leads to the least compact
variational space $P$ among the three considered ranking criteria. However,
the situation gets improved steadily by reducing $C_{min}$. More importantly,
the iCI selection is much more efficient than the other two (cf. Fig. 5): the
cost of the iCI selection is only 20% of that of PT2, whereas the other two
are approximately two times more expensive than their PT2’s. Therefore, the
slight loss of compactness is overcompensated by the gain in efficiency. An
additional point deserves to be mentioned here: while typically less than 1%
of the CSFs in the CIPSI and iCI variational spaces have coefficients somewhat
smaller in absolute value than $C_{min}$, such CSFs can be up to 10% in the
pruned HBCI variational space.

|
---|---
(a) $E_{var}$ at $R_{eq}$ | (b) $E_{var}$ at $2.0R_{eq}$
Figure 3: The all-electron variational energy $E_{var}$ of Cr2/SV as function
of the number $N_{csf}$ of CSFs selected by the CIPSI (24), iCI (28)/(29) and
HBCI (30)/(31) ranking criteria. Figure 4: The size ratio
$N_{csf}^{\prime}/N_{csf}$ between the unpruned and pruned variational spaces
determined by the CIPSI (24), iCI (28)/(29) and HBCI (30)/(31) ranking
criteria.  Figure 5: The time ratio $T_{var}/T_{pt}$ between selection and PT2
by the CIPSI (24), iCI (28)/(29) and HBCI (30)/(31) ranking criteria.

###  7.2 Cr2

Given the improved efficiency of iCIPT2, we are now ready to perform more
accurate calculations on the ground state of Cr2, a classic strongly
correlated system for testing various methods101, 51, 94, 102, 62, 63, 18,
103. To this end, the spin-free exact two-component (sf-X2C) relativistic
Hamiltonian104, 105 and the cc-pVDZ-DK basis106 were adopted, along with an
active space of (28e, 76o) that is composed of $4.8\times 10^{27}$ CSFs. While
the Ne-core orbitals are just HF orbitals, the NOs were used for correlation,
which were generated in two ways, (A) ‘partially dynamic’ (see point (4) in
the end of Sec. 4) and (B) ‘static’. In the latter, the NOs were generated
with $C_{min}=1.0\times 10^{-4}$ and then used for all smaller $C_{min}$’s. In
this case, the variational space determined at a $C_{min}$ can directly be
used as initial guess for a smaller $C_{min}$. The results with the two sets
of NOs are documented in Tables 5 and 6, respectively. It is of interest to
see that the extrapolated value (-2099.9223(4) $E_{h}$) by scheme B is very
close to the HBCI value (-2099.9224(6) $E_{h}$)102.However, it is believed
that the extrapolated value (-2099.9240(2) $E_{h}$) by scheme A is more
accurate (by $-1.7$ $mE_{h}$) than that by scheme (B) because the partially
dynamically generated NOs in the former are of better quality than the fixed
ones in the latter. It is for sure that the two sets of calculations would
agree with each other completely by further reducing $C_{min}$, so as to
sample even larger portions of the full Hilbert space and hence minimize
orbital rotation effects. For instance, in the all-electron nonrelativistic
calculations on Cr2 with the SV basis, the difference between the energies by
scheme A (-2086.44468(2) $E_{h}$) and scheme B (-2086.44466(2) $E_{h}$) is
only $-0.02$ $mE_{h}$. Here, the largest FOIS ($1.3\times 10^{11}$) sampled by
iCIPT2/SV is ca. $10^{-8}$% of the full space ($1.4\times 10^{21}$). In
contrast, the largest FOIS ($3.2\times 10^{12}$) sampled by sf-X2C-iCIPT2/cc-
pVDZ-DK is only $10^{-14}\%$ of the full space ($4.8\times 10^{27}$).

Table 5: Frozen-core sf-X2C-iCIPT2/cc-pVDZ-DK calculations of the ground state
of Cr2 at $R_{eq}=1.68$ Å a

C${}_{\text{min}}$b |  $N_{cfg}$c |  $N_{csf}$d |  $\tilde{N}_{csf}$e |  $\tilde{N}_{det}$f | $E_{var}(+2099)/E_{h}$ | $E_{tot}(+2099)/E_{h}$ | $T_{wall}/s$
---|---|---|---|---|---|---|---
$5.0\times 10^{-5}$ | 455404 | 3136853 | 896527 | 3666946 | -0.81095 | -0.90466 | 464
$4.0\times 10^{-5}$ | 615595 | 4512684 | 1250882 | 5197093 | -0.81929 | -0.90612 | 641
$3.0\times 10^{-5}$ | 903347 | 7115256 | 1904606 | 8059852 | -0.82909 | -0.90785 | 973
$2.0\times 10^{-5}$ | 1560133 | 13369034 | 3470331 | 15033913 | -0.84161 | -0.90100 | 1775
$1.5\times 10^{-5}$ | 2281773 | 20737948 | 5268475 | 23176915 | -0.84926 | -0.91128 | 2765
$1.0\times 10^{-5}$ | 3959566 | 38959372 | 9616690 | 43154698 | -0.85945 | -0.91300 | 5194
$9.0\times 10^{-6}$ | 4565597 | 45851900 | 11235054 | 50662572 | -0.86191 | -0.91341 | 6309
0.0g |  |  |  |  -2099.9240$\pm$0.0002 |

  * a

D2h symmetry; HF energy: -2098.536329 $E_{h}$; active space: (28e, 76o);
partially dynamically generated NOs (see point (4) in the end of Sec. 4).

  * b

Threshold for pruning CSF in the variational (var) space.

  * c

Number of orbital configurations in the variational space.

  * d

Number of CSFs corresponding to $N_{cfg}$.

  * e

Number of CSFs after selection, among which about 0.5% have coefficients
slightly smaller in absolute values than $C_{\mathrm{min}}$ (which is due to a
final diagonalization).

  * f

Estimated number of determinants according to the expression
$\sum_{I}\frac{\tilde{N}_{\mathrm{csf}}^{I}}{N_{\mathrm{csf}}^{I}}N_{\mathrm{det}}^{I}$,
with $N_{\mathrm{det}}^{I}$ being the numbers of determinants of CFG
$|I\rangle$.

  * g

Extrapolated value by linear fit of the $E_{\mathrm{total}}$ vs.
$|E_{c}^{(2)}|$ plot, with uncertainty being half the length of 95% confidence
interval.

Table 6: Frozen-core sf-X2C-iCIPT2/cc-pVDZ-DK calculations of the ground state
of Cr2 at $R_{eq}=1.68$ Å a

C${}_{\text{min}}$ | $N_{cfg}$ | $N_{csf}$ | $\tilde{N}_{csf}$ | $\tilde{N}_{det}$ | $E_{var}(+2099)/E_{h}$ | $E_{tot}(+2099)/E_{h}$ | $T_{wall}/s$
---|---|---|---|---|---|---|---
5.0$\times 10^{-5}$ | 425431 | 2892927 | 840568 | 3438022 | -0.809883 | -0.903447 | 434
4.0$\times 10^{-5}$ | 573307 | 4145556 | 1168810 | 4857157 | -0.818110 | -0.904792 | 599
3.0$\times 10^{-5}$ | 840891 | 6513720 | 1779100 | 7530355 | -0.827922 | -0.906404 | 916
2.0$\times 10^{-5}$ | 1449299 | 12182405 | 3233090 | 14006309 | -0.840470 | -0.908481 | 1651
1.5$\times 10^{-5}$ | 2153481 | 19286806 | 4989382 | 21922383 | -0.847846 | -0.909859 | 2645
1.0$\times 10^{-5}$ | 3716908 | 35893887 | 9034413 | 40531897 | -0.858790 | -0.911610 | 4978
9.0$\times 10^{-6}$ | 4284579 | 42309173 | 10542471 | 47528961 | -0.861208 | -0.912028 | 5608
0.0 |  |  |  |  -2099.9223$\pm$0.0004 |
0.0b |  |  |  |  -2099.9224$\pm$0.0006 |
0.0c |  |  |  |  -2099.9195$\pm$0.0027 |

  * a

Fixed NOs generated with $C_{min}=1.0\times 10^{-4}$. For additional
explanations see Table 5.

  * b

HBCI102.

  * c

DMRG ($M=16000$)103.

###  7.3 [2Fe-2S]

As a final example, we consider [Fe2S2(SCH3)]42- (abbreviated as [2Fe-2S]).
Albeit a simplest system for modeling iron-containing enzymes19, even the
minimal chemically meaningful active space (30e, 20o) already contains
$240374016$ CSFs ($52581816$, $99419400$, $64438500$, $20575100$, $3174444$,
and $184756$ CSFs for spins from 0 to 5, respectively). The 20 active orbitals
here include Fe $3d$ and S $3p$ of the core [2Fe-2S] as well as their four
$\sigma$-bonds with ligands. To make a direct comparison with the previous
work, the (nonrelativistic) molecular orbital integrals are taken simply from
Ref. 19. The extrapolated (valence) energy of the singlet ground state is
-116.60574(6) $E_{h}$ (cf. Table 7), very close to the DMRG value of
-116.60561 $E_{h}$ with $M=8000$107. It is more appealing to take a closer
look at the wave function (selected with $C_{min}=5.0\times 10^{-6}$): it is
composed of 41.13% Fe(II) – Fe(II), 22.69% Fe(III) – Fe(II), 22.01% Fe(II) –
Fe(I), 6.05% Fe(III) – Fe(I), 3.40% Fe(III) – Fe(III), and 2.86% Fe(I) –
Fe(I), but the leading CSF has a weight only of 3.36% and stems from Fe(III) –
Fe(III) instead of the leading structure Fe(II) – Fe(II). It is not clear
whether this peculiar picture will hold for a larger active space, especially
in conjunction with localized orbitals. Nevertheless, it is perfectly
legitimate to further calculate more spin states, just to reveal the efficacy
of iCIPT2.

Although the particle-hole algorithm presented in Sec. 5 allows for
simultaneous calculations of several states of different spins with a common
set of orbitals, the function is not yet at our disposal. Instead, states of
different spins are calculated here separately. One then has to face the issue
that the quality of states of different spins may be different for the same
$C_{min}$, so as to affect the relative energies. To circumvent this problem,
we decompose the vertical excitation energy $\Delta E_{S,i}$
($=E_{S,i}-E_{0,0}$) of state $i$ with spin $S$ into two terms, $\Delta
E_{S,i}-\Delta E_{S,0}$ and $\Delta E_{S,0}-\Delta E_{0,0}$. The former is
calculated with the same $C_{min}$, whereas the latter is set to the
extrapolated value, i.e.,

$\Delta
E_{S,i}(0.0)=\left[E_{S,i}(C_{min})-E_{S,0}(C_{min})\right]+\left[E_{S,0}(0.0)-E_{0,0}(0.0)\right].$
(38)

In essence, what is assumed here is that different states of the same spin
have the same extrapolation distance. The so-calculated excitation energies
for all states below 1.8 eV are documented in Table 8 and further plotted in
Fig. 6. Interestingly, the lowest 6 states originate just from the lowest
state of each of the 6 possible spins $S\in[0,5]$ that result from the
coupling of the two high-spin-$\frac{5}{2}$ ionic configurations Fe(III)
$3d^{5}$ and are energetically ordered in increasing spin. However, the
remaining 30 states of different spins are intersected within an interval of
0.91 eV. That is, they are separated in average by only 0.03 eV, just like the
gap between the lowest singlet and triplet states. It is clear that such
complicated spin structure cannot be handled by determinant-based methods due
to severe spin contaminations.

Table 7: The singlet ground state of [Fe2S2(SCH3)]42- (for additional explanations see Table 5) C${}_{\text{min}}$ | $N_{cfg}$ | $N_{csf}$ | $\tilde{N}_{csf}$ | $\tilde{N}_{det}$ | $E_{var}/E_{h}$ | $E_{tot}/mE_{h}$
---|---|---|---|---|---|---
$2.0\times 10^{-5}$ | 313428 | 6242941 | 568382 | 3056319 | -116.602301 | -116.604325
$1.5\times 10^{-5}$ | 396163 | 7436740 | 754898 | 4031212 | -116.603153 | -116.604691
$1.0\times 10^{-5}$ | 540358 | 9332213 | 1108867 | 5866360 | -116.604010 | -116.605046
$9.0\times 10^{-6}$ | 586437 | 9908183 | 1228988 | 6485715 | -116.604196 | -116.605120
$7.0\times 10^{-6}$ | 744932 | 11779934 | 1673904 | 8765188 | -116.604554 | -116.605256
$5.0\times 10^{-6}$ | 938275 | 13485301 | 2236247 | 11619598 | -116.604908 | -116.605385
0.0 |  |  |  |  |  -116.60574$\pm$0.00006
Table 8: Vertical excitation energies ($\Delta E$ in eV) of low-lying excited states of [2Fe-2S] State | $2S+1$ | $\Delta E$ | State | $2S+1$ | $\Delta E$
---|---|---|---|---|---
0 | 1 | 0.000 | 18 | 5 | 1.502
1 | 3 | 0.034 | 19 | 7 | 1.578
2 | 5 | 0.112 | 20 | 9 | 1.605
3 | 7 | 0.246 | 21 | 3 | 1.635
4 | 9 | 0.454 | 22 | 9 | 1.645
5 | 11 | 0.767 | 23 | 7 | 1.652
6 | 3 | 0.889 | 24 | 5 | 1.664
7 | 5 | 1.044 | 25 | 3 | 1.701
8 | 5 | 1.159 | 26 | 5 | 1.706
9 | 7 | 1.225 | 27 | 1 | 1.718
10 | 7 | 1.269 | 28 | 3 | 1.727
11 | 9 | 1.340 | 29 | 7 | 1.738
12 | 3 | 1.376 | 30 | 1 | 1.774
13 | 7 | 1.391 | 31 | 9 | 1.774
14 | 5 | 1.409 | 32 | 5 | 1.776
15 | 9 | 1.434 | 33 | 9 | 1.786
16 | 7 | 1.480 | 34 | 3 | 1.787
17 | 9 | 1.482 | 35 | 5 | 1.803
Figure 6: Vertical excitation energies of low-lying excited states of [2Fe-2S]

##  8 Conclusions and Outlook

Sticking to the parlance of static and dynamic correlations, there could be
two paradigms for handling strongly correlated systems of electrons, viz.
“more static, less dynamic” and “less static, more dynamic”. It is the former
that is followed by the family of sCI+PT2 methods, among which iCIPT2 stands
out in several aspects as already highlighted in the Introduction. Its
efficiency has been improved herein by up to $20\times$ via three major
techniques: a new criterion for configuration selection, a new particle-hole
algorithm for Hamiltonian construction, and a new data structure for the quick
sorting of the FOIS, in addition to the workhorse TUGA for computing and
reutilizing the basic coupling coefficients between randomly selected CSFs.
All these can be applied to other types of sCI+PT2 as well. As revealed by
several examples, iCIPT2 can indeed be characterized as an near-exact
approach. However, it is still memory intensive, even though the FOIS has been
decomposed into disjoint subspaces. A possible way to resolve this issue is to
evaluate the PT2 correction in a stochastic manner. Other immediate extensions
of iCIPT2 include simultaneous treatment of states of different spatial and/or
spin symmetries (which is furnished by the particle-hole algorithm),
perturbative treatment of spin-orbit couplings (via the sf-X2C+soc-DKH1
Hamiltonian104, 105), and direct access of high-lying states of a given energy
window (which is furnished by the iVI eigensolver83, 85), etc. Work along
these directions are being undertaken at our laboratory.

The research of this work was supported by National Natural Science Foundation
of China (Grant Nos. 21833001 and 21973054) and the North Dakota University
System.

## References

  * Liu 2020 Liu, W. Essentials of relativistic quantum chemistry. _J. Chem. Phys._ 2020, _152_ , 180901
  * Liu 2020 Liu, W. Relativistic quantum chemistry: today and tomorrow. _Sci. Sin. Chim._ 2020, DOI: 10.1360/SSC–2020–0120
  * Chaudhuri et al. 2005 Chaudhuri, R. K.; Freed, K. F.; Hose, G.; Piecuch, P.; Kowalski, K.; Włoch, M.; Chattopadhyay, S.; Mukherjee, D.; Rolik, Z.; Szabados, Á.; Tóth, G.; Surján, P. R. Comparison of low-order multireference many-body perturbation theories. _J. Chem. Phys._ 2005, _122_ , 134105
  * Hoffmann et al. 2009 Hoffmann, M. R.; Datta, D.; Das, S.; Mukherjee, D.; Szabados, Á.; Rolik, Z.; Surján, P. R. Comparative study of multireference perturbative theories for ground and excited states. _J. Chem. Phys._ 2009, _131_ , 204104
  * Lyakh et al. 2012 Lyakh, D. I.; Musiaz, M.; Lotrich, V. F.; Bartlett, R. J. Multireference nature of chemistry: The coupled-cluster view. _Chem. Rev._ 2012, _112_ , 182–243
  * Sherrill and Schaefer III 1999 Sherrill, C. D.; Schaefer III, H. F. The configuration interaction method: Advances in highly correlated approaches. _Adv. Quantum Chem._ 1999, _34_ , 143–269
  * Szalay et al. 2012 Szalay, P. G.; Muller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfiguration self-consistent field and multireference configuration interaction methods and applications. _Chem. Rev._ 2012, _112_ , 108–181
  * Lischka et al. 2018 Lischka, H.; Nachtigallova, D.; Aquino, A. J.; Szalay, P. G.; Plasser, F.; Machado, F. B.; Barbatti, M. Multireference approaches for excited states of molecules. _Chem. Rev._ 2018, _118_ , 7293–7361
  * Liu and Hoffmann 2016 Liu, W.; Hoffmann, M. R. iCI: Iterative CI toward full CI. _J. Chem. Theory Comput._ 2016, _12_ , 1169–1178, (E) 2016, _12_ , 3000
  * White 1992 White, S. R. Density matrix formulation for quantum renormalization groups. _Phys. Rev. Lett._ 1992, _69_ , 2863
  * White 1993 White, S. R. Density-matrix algorithms for quantum renormalization groups. _Phys. Rev. B_ 1993, _48_ , 10345–10356
  * White and Martin 1999 White, S. R.; Martin, R. L. Ab initio quantum chemistry using the density matrix renormalization group. _J. Chem. Phys._ 1999, _110_ , 4127–4130
  * Mitrushenkov et al. 2001 Mitrushenkov, A. O.; Fano, G.; Ortolani, F.; Linguerri, R.; Palmieri, P. Quantum chemistry using the density matrix renormalization group. _J. Chem. Phys._ 2001, _115_ , 6815–6821
  * Chan and Head-Gordon 2002 Chan, G. K.-L.; Head-Gordon, M. Highly correlated calculations with a polynomial cost algorithm: A study of the density matrix renormalization group. _J. Chem. Phys._ 2002, _116_ , 4462–4476
  * Chan and Sharma 2011 Chan, G. K.-L.; Sharma, S. The density matrix renormalization group in quantum chemistry. _Annu. Rev. Phys. Chem._ 2011, _62_ , 465–481
  * Sharma and Chan 2012 Sharma, S.; Chan, G. K.-L. Spin-adapted density matrix renormalization group algorithms for quantum chemistry. _J. Chem. Phys._ 2012, _136_ , 124121
  * Wouters and Van Neck 2014 Wouters, S.; Van Neck, D. The density matrix renormalization group for ab initio quantum chemistry. _Eur. Phys. J. D_ 2014, _68_ , 272
  * Olivares-Amaya et al. 2015 Olivares-Amaya, R.; Hu, W.; Nakatani, N.; Sharma, S.; Yang, J.; Chan, G. K.-L. The ab-initio density matrix renormalization group in practice. _J. Chem. Phys._ 2015, _142_ , 034102
  * Li and Chan 2017 Li, Z.; Chan, G. K.-L. Spin-Projected Matrix Product States: Versatile Tool for Strongly Correlated Systems. _J. Chem. Theory Comput._ 2017, _13_ , 2681–2695
  * Baiardi and Reiher 2020 Baiardi, A.; Reiher, M. The density matrix renormalization group in chemistry and molecular physics: Recent developments and new challenges. _J. Chem. Phys._ 2020, _152_ , 040903
  * Booth et al. 2009 Booth, G. H.; Thom, A. J. W.; Alavi, A. Fermion Monte Carlo without fixed nodes: A game of life, death, and annihilation in Slater determinant space. _J. Chem. Phys._ 2009, _131_ , 054106
  * Cleland et al. 2010 Cleland, D.; Booth, G. H.; Alavi, A. Communications: Survival of the fittest: Accelerating convergence in full configuration-interaction quantum Monte Carlo. _J. Chem. Phys._ 2010, _132_ , 041103
  * Blunt et al. 2017 Blunt, N.; Booth, G. H.; Alavi, A. Density matrices in full configuration interaction quantum Monte Carlo: Excited states, transition dipole moments, and parallel distribution. _J. Chem. Phys._ 2017, _146_ , 244105
  * Petruzielo et al. 2012 Petruzielo, F. R.; Holmes, A. A.; Changlani, H. J.; Nightingale, M. P.; Umrigar, C. J. Semistochastic projector monte carlo method. _Phys. Rev. Lett._ 2012, _109_ , 230201
  * Holmes et al. 2016 Holmes, A. A.; Changlani, H. J.; Umrigar, C. J. Efficient heat-bath sampling in Fock space. _J. Chem. Theory Comput._ 2016, _12_ , 1561–1571
  * Ten-no 2013 Ten-no, S. Stochastic determination of effective Hamiltonian for the full configuration interaction solution of quasi-degenerate electronic states. _J. Chem. Phys._ 2013, _138_ , 164126
  * Ghanem et al. 2019 Ghanem, K.; Lozovoi, A. Y.; Alavi, A. Unbiasing the initiator approximation in full configuration interaction quantum Monte Carlo. _J. Chem. Phys._ 2019, _151_ , 224108
  * Blunt 2019 Blunt, N. S. A hybrid approach to extending selected configuration interaction and full configuration interaction quantum Monte Carlo. _J. Chem. Phys._ 2019, _151_ , 174103
  * Dobrautz et al. 2019 Dobrautz, W.; Smart, S. D.; Alavi, A. Efficient formulation of full configuration interaction quantum Monte Carlo in a spin eigenbasis via the graphical unitary group approach. _J. Chem. Phys._ 2019, _151_ , 094104
  * Neufeld and Thom 2020 Neufeld, V. A.; Thom, A. J. Accelerating Convergence in Fock Space Quantum Monte Carlo Methods. _J. Chem. Theory Comput._ 2020, _16_ , 1503–1510
  * Deustua et al. 2018 Deustua, J. E.; Magoulas, I.; Shen, J.; Piecuch, P. Communication: Approaching exact quantum chemistry by cluster analysis of full configuration interaction quantum Monte Carlo wave functions. _J. Chem. Phys._ 2018, _149_ , 151101
  * Lyakh and Bartlett 2010 Lyakh, D. I.; Bartlett, R. J. An adaptive coupled-cluster theory:@CC approach. _J. Chem. Phys._ 2010, _133_ , 244112
  * Xu et al. 2018 Xu, E.; Uejima, M.; Ten-no, S. L. Full Coupled-Cluster Reduction for Accurate Description of Strong Electron Correlation. _Phys. Rev. Lett._ 2018, _121_ , 113001
  * Xu et al. 2020 Xu, E.; Uejima, M.; Ten-No, S. L. Towards near-exact solutions of molecular electronic structure: Full coupled-cluster reduction with a second-order perturbative correction. _J. Phys. Chem. Lett._ 2020, _11_ , 9775–9780
  * Zhang and Evangelista 2016 Zhang, T.; Evangelista, F. A. A deterministic projector configuration interaction approach for the ground state of quantum many-body systems. _J. Chem. Theory Comput._ 2016, _12_ , 4326–4337
  * Eriksen et al. 2017 Eriksen, J. J.; Lipparini, F.; Gauss, J. Virtual orbital many-body expansions: A possible route towards the full configuration interaction limit. _J. Phys. Chem. Lett._ 2017, _8_ , 4633–4639
  * Eriksen and Gauss 2018 Eriksen, J. J.; Gauss, J. Many-Body Expanded Full Configuration Interaction. I. Weakly Correlated Regime. _J. Chem. Theory Comput._ 2018, _14_ , 5180–5191
  * Eriksen and Gauss 2019 Eriksen, J. J.; Gauss, J. Many-Body Expanded Full Configuration Interaction. II. Strongly Correlated Regime. _J. Chem. Theory Comput._ 2019, _15_ , 4873–4884
  * Eriksen and Gauss 2019 Eriksen, J. J.; Gauss, J. Generalized Many-Body Expanded Full Configuration Interaction Theory. _J. Phys. Chem. Lett._ 2019, _10_ , 7910–7915
  * Zimmerman 2017 Zimmerman, P. M. Incremental full configuration interaction. _J. Chem. Phys._ 2017, _146_ , 104102
  * Zimmerman 2017 Zimmerman, P. M. Strong correlation in incremental full configuration interaction. _J. Chem. Phys._ 2017, _146_ , 224104
  * Fales et al. 2018 Fales, B. S.; Seritan, S.; Settje, N. F.; Levine, B. G.; Koch, H.; Martínez, T. J. Large-Scale Electron Correlation Calculations: Rank-Reduced Full Configuration Interaction. _J. Chem. Theory Comput._ 2018, _14_ , 4139–4150
  * Greene et al. 2019 Greene, S. M.; Webber, R. J.; Weare, J.; Berkelbach, T. C. Beyond walkers in stochastic quantum chemistry: reducing error using fast randomized iteration. _J. Chem. Theory Comput._ 2019, _15_ , 4834–4850
  * Bytautas and Ruedenberg 2004 Bytautas, L.; Ruedenberg, K. Correlation energy extrapolation by intrinsic scaling. I. Method and application to the neon atom. _J. Chem. Phys._ 2004, _121_ , 10905–10918
  * Bytautas and Ruedenberg 2004 Bytautas, L.; Ruedenberg, K. Correlation energy extrapolation by intrinsic scaling. II. The water and the nitrogen molecule. _J. Chem. Phys._ 2004, _121_ , 10919–10934
  * Greer 1998 Greer, J. C. Monte Carlo configuration interaction. _J. Comput. Phys._ 1998, _146_ , 181–202
  * Coe and Paterson 2012 Coe, J. P.; Paterson, M. J. Development of Monte Carlo configuration interaction: Natural orbitals and second-order perturbation theory. _J. Chem. Phys._ 2012, _137_ , 204108
  * Ohtsuka and Hasegawa 2017 Ohtsuka, Y.; Hasegawa, J.-y. Selected configuration interaction method using sampled first-order corrections to wave functions. _J. Chem. Phys._ 2017, _147_ , 034102
  * Coe 2018 Coe, J. P. Machine Learning Configuration Interaction. _J. Chem. Theory Comput._ 2018, _14_ , 5739–5749
  * Holmes et al. 2016 Holmes, A. A.; Tubman, N. M.; Umrigar, C. J. Heat-bath configuration interaction: an efficient selected configuration interaction algorithm inspired by heat-bath sampling. _J. Chem. Theory Comput._ 2016, _12_ , 3674–3680
  * Garniron et al. 2017 Garniron, Y.; Scemama, A.; Loos, P.-F.; Caffarel, M. Hybrid stochastic-deterministic calculation of the second-order perturbative contribution of multireference perturbation theory. _J. Chem. Phys._ 2017, _147_ , 034101
  * Holmes et al. 2017 Holmes, A. A.; Umrigar, C. J.; Sharma, S. Excited states using semistochastic heat-bath configuration interaction. _J. Chem. Phys._ 2017, _147_ , 164111
  * Sharma et al. 2017 Sharma, S.; Holmes, A. A.; Jeanmairet, G.; Alavi, A.; Umrigar, C. J. Semistochastic Heat-Bath Configuration Interaction Method: Selected Configuration Interaction with Semistochastic Perturbation Theory. _J. Chem. Theory Comput._ 2017, _13_ , 1595–1604
  * Chien et al. 2018 Chien, A. D.; Holmes, A. A.; Otten, M.; Umrigar, C. J.; Sharma, S.; Zimmerman, P. M. Excited States of Methylene, Polyenes, and Ozone from Heat-Bath Configuration Interaction. _J. Phys. Chem. A_ 2018, _122_ , 2714–2722
  * Yao et al. 2020 Yao, Y.; Giner, E.; Li, J.; Toulouse, J.; Umrigar, C. Almost exact energies for the Gaussian-2 set with the semistochastic heat-bath configuration interaction method. _J. Chem. Phys._ 2020, _153_ , 124117
  * Evangelista 2014 Evangelista, F. A. Adaptive multiconfigurational wave functions. _J. Chem. Phys._ 2014, _140_ , 124114
  * Schriber and Evangelista 2016 Schriber, J. B.; Evangelista, F. A. Communication: An adaptive configuration interaction approach for strongly correlated electrons with tunable accuracy. _J. Chem. Phys._ 2016, _144_ , 161106
  * Schriber and Evangelista 2017 Schriber, J. B.; Evangelista, F. A. Adaptive configuration interaction for computing challenging electronic excited states with tunable accuracy. _J. Chem. Theory Comput._ 2017, _13_ , 5354–5366
  * Schriber et al. 2018 Schriber, J. B.; Hannon, K. P.; Li, C.; Evangelista, F. A. A Combined Selected Configuration Interaction and Many-Body Treatment of Static and Dynamical Correlation in Oligoacenes. _J. Chem. Theory Comput._ 2018, _14_ , 6295–6305
  * Tubman et al. 2016 Tubman, N. M.; Lee, J.; Takeshita, T. Y.; Head-Gordon, M.; Whaley, K. B. A deterministic alternative to the full configuration interaction quantum Monte Carlo method. _J. Chem. Phys._ 2016, _145_ , 044112
  * Lehtola et al. 2017 Lehtola, S.; Tubman, N. M.; Whaley, K. B.; Head-Gordon, M. Cluster decomposition of full configuration interaction wave functions: A tool for chemical interpretation of systems with strong correlation. _J. Chem. Phys._ 2017, _147_ , 154105
  * Tubman et al. 2020 Tubman, N. M.; Freeman, C. D.; Levine, D. S.; Hait, D.; Head-Gordon, M.; Whaley, K. B. Modern approaches to exact diagonalization and selected configuration interaction with the adaptive sampling CI method. _J. Chem. Theory Comput._ 2020, _16_ , 2139–2159
  * Tubman et al. 2018 Tubman, N. M.; Levine, D. S.; Hait, D.; Head-Gordon, M.; Whaley, K. B. An efficient deterministic perturbation theory for selected configuration interaction methods. 2018, arXiv preprint arXiv:1808.02049
  * Liu and Hoffmann 2014 Liu, W.; Hoffmann, M. R. SDS: the ‘static-dynamic-static’ framework for strongly correlated electrons. _Theor. Chem. Acc._ 2014, _133_ , 1481
  * Lei et al. 2017 Lei, Y.; Liu, W.; Hoffmann, M. R. Further development of SDSPT2 for strongly correlated electrons. _Mol. Phys._ 2017, _115_ , 2696–2707
  * Zhang et al. 2020 Zhang, N.; Liu, W.; Hoffmann, M. R. Iterative Configuration Interaction with Selection. _J. Chem. Theory Comput._ 2020, _16_ , 2296–2316
  * Garniron et al. 2018 Garniron, Y.; Scemama, A.; Giner, E.; Caffarel, M.; Loos, P.-F. Selected configuration interaction dressed by perturbation. _J. Chem. Phys._ 2018, _149_ , 064103
  * Wang et al. 2019 Wang, Z.; Li, Y.; Lu, J. Coordinate descent full configuration interaction. _J. Chem. Theory Comput._ 2019, _15_ , 3558–3569
  * Giner et al. 2016 Giner, E.; Assaraf, R.; Toulouse, J. Quantum Monte Carlo with reoptimised perturbatively selected configuration-interaction wave functions. _Mol. Phys._ 2016, _114_ , 910–920
  * Scemama et al. 2018 Scemama, A.; Benali, A.; Jacquemin, D.; Caffarel, M.; Loos, P.-F. Excitation energies from diffusion Monte Carlo using selected configuration interaction nodes. _J. Chem. Phys._ 2018, _149_ , 034108
  * Bender and Davidson 1969 Bender, C. F.; Davidson, E. R. Studies in configuration interaction: The first-row diatomic hydrides. _Phys. Rev._ 1969, _183_ , 23
  * Whitten and Hackmeyer 1969 Whitten, J. L.; Hackmeyer, M. Configuration interaction studies of ground and excited states of polyatomic molecules. I. The CI formulation and studies of formaldehyde. _J. Chem. Phys._ 1969, _51_ , 5584–5596
  * Huron et al. 1973 Huron, B.; Malrieu, J. P.; Rancurel, P. Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wave functions. _J. Chem. Phys._ 1973, _58_ , 5745–5759
  * Evangelisti et al. 1983 Evangelisti, S.; Daudey, J. P.; Malrieu, J. P. Convergence of an improved CIPSI algorithm. _Chem. Phys._ 1983, _75_ , 91–102
  * Buenker and Peyerimhoff 1974 Buenker, R. J.; Peyerimhoff, S. D. Individualized configuration selection in CI calculations with subsequent energy extrapolation. _Theor. Chem. Acta_ 1974, _35_ , 33–58
  * Eriksen et al. 2020 Eriksen, J. J.; Anderson, T. A.; Deustua, J. E.; Ghanem, K.; Hait, D.; Hoffmann, M. R.; Lee, S.; Levine, D. S.; Magoulas, I.; Shen, J.; Tubman, N. M.; Whaley, K. B.; Xu, E.; Yao, Y.; Zhang, N.; Alavi, A.; Chan, G. K.-L.; Head-Gordon, M.; Liu, W.; Piecuch, P.; Sharma, S.; Ten-no, S. L.; Umrigar, C. J.; Gauss, J. The ground state electronic energy of benzene. _J. Phys. Chem. Lett._ 2020, _11_ , 8922–8929
  * 77 Detail of timsort. https://en.wikipedia.org/wiki/Timsort
  * 78 https://github.com/timsort/cpp-TimSort
  * Buenker 1980 Buenker, R. J. In _Molecular Physics and Quantum Chemistry: into the 80’s_ ; Burton, P. G., Ed.; University of Wollongong Press: Wollongong, Australia, 1980; pp 1.5.1–1.5.37
  * Buenker 1986 Buenker, R. J. Combining perturbation theory techniques with variational CI calculations to study molecular excited states. _Int. J. Quantum Chem._ 1986, _29_ , 435–460
  * Krebs and Buenker 1995 Krebs, S.; Buenker, R. J. A new table-direct configuration interaction method for the evaluation of Hamiltonian matrix elements in a basis of linear combinations of spin-adapted functions. _J. Chem. Phys._ 1995, _103_ , 5613–5629
  * 82 Khait, Y. G.; Hoffmann, M. R. _Low-Lying Potential Energy Surfaces_ ; Chapter 8, pp 176–198
  * Huang et al. 2017 Huang, C.; Liu, W.; Xiao, Y.; Hoffmann, M. R. iVI: An iterative vector interaction method for large eigenvalue problems. _J. Comput. Chem._ 2017, _38_ , 2481–2499, (E) 2018, _39_ , 338
  * 84 Song, Y.; Guo, Y.; Lei, Y.; Liu, W. Benchmarking SDSPT2 and SDSCI for low-lying excited states of closed- and open-shell oragnic molecules. (unpublished).
  * Huang and Liu 2019 Huang, C.; Liu, W. iVI-TD-DFT: An iterative vector interaction method for exterior/interior roots of TD-DFT. _J. Comput. Chem._ 2019, _40_ , 1023–1037, (E) 2018, _39_ , 338
  * Paldus 1974 Paldus, J. Group theoretical approach to the configuration interaction and perturbation theory calculations for atomic and molecular systems. _J. Chem. Phys._ 1974, _61_ , 5321–5330
  * Paldus and Boyle 1980 Paldus, J.; Boyle, M. J. Unitary group approach to the many-electron correlation problem via graphical methods of spin algebras. _Phys. Script._ 1980, _21_ , 295–311
  * Shavitt 1981 Shavitt, I. In _The Unitary Group for the Evaluation of Electronic Energy Matrix Elements_ , 1st ed.; Hinze, J., Ed.; Lecture Notes in Chemistry 22; Springer-Verlag Berlin Heidelberg, 1981; pp 1–50
  * Shavitt 1977 Shavitt, I. Graph theoretical concepts for the unitary group approach to the many-electron correlation problem. _Int. J. Quantum Chem., Symp._ 1977, _12_ , 131–148
  * Duch and Karwowski 1985 Duch, W.; Karwowski, J. Symmetric group approach to configuration interaction methods. _Comp. Phys. Rep._ 1985, _2_ , 93–170
  * Shavitt 1978 Shavitt, I. Matrix element evaluation in the unitary group approach to the electron correlation problem. _Int. J. Quantum Chem._ 1978, _14_ , 5–32
  * Stampfuß and Wenzel 2005 Stampfuß, P.; Wenzel, W. Improved implementation and application of the individually selecting configuration interaction method. _J. Chem. Phys._ 2005, _122_ , 024110
  * Smith et al. 2017 Smith, J. E. T.; Mussard, B.; Holmes, A. A.; Sharma, S. Cheap and near exact CASSCF with large active spaces. _J. Chem. Theory Comput._ 2017, _13_ , 5468–5478
  * Li et al. 2020 Li, J.; Yao, Y.; Holmes, A. A.; Otten, M.; Sun, Q.; Sharma, S.; Umrigar, C. J. Accurate many-body electronic structure near the basis set limit: Application to the chromium dimer. _Phys. Rev. Research_ 2020, _2_ , 012015(R)
  * Schäfer et al. 1992 Schäfer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. _J. Chem. Phys._ 1992, _97_ , 2571–2577
  * Liu et al. 1997 Liu, W.; Hong, G.; Dai, D.; Li, L.; Dolg, M. The Beijing 4-component density functional theory program package (BDF) and its application to EuO, EuS, YbO and YbS. _Theor. Chem. Acc._ 1997, _96_ , 75–83
  * Liu et al. 2003 Liu, W.; Wang, F.; Li, L. _J. Theor. Comput. Chem._ 2003, _2_ , 257–272
  * Liu et al. 2004 Liu, W.; Wang, F.; Li, L. In _Recent Advances in Relativistic Molecular Theory_ ; Hirao, K., Ishikawa, Y., Eds.; World Scientific: Singapore, 2004; pp 257–282
  * Liu et al. 2004 Liu, W.; Wang, F.; Li, L. In _Encyclopedia of Computational Chemistry_ ; von Ragué Schleyer, P., Allinger, N. L., Clark, T., Gasteiger, J., Kollman, P. A., Schaefer III, H. F., Eds.; Wiley: Chichester, UK, 2004
  * Zhang et al. 2020 Zhang, Y.; Suo, B.; Wang, Z.; Zhang, N.; Li, Z.; Lei, Y.; Zou, W.; Gao, J.; Peng, D.; Pu, Z.; Xiao, Y.; Sun, Q.; Wang, F.; Ma, Y.; Wang, X.; Guo, Y.; Liu, W. BDF: A relativistic electronic structure program package. _J. Chem. Phys._ 2020, _152_ , 064113
  * Booth et al. 2014 Booth, G. H.; Smart, S. D.; Alavi, A. Linear-scaling and parallelisable algorithms for stochastic quantum chemistry. _Mol. Phys._ 2014, _112_ , 1855–1869
  * Li et al. 2018 Li, J.; Otten, M.; Holmes, A. A.; Sharma, S.; Umrigar, C. J. Fast semistochastic heat-bath configuration interaction. _J. Chem. Phys._ 2018, _149_ , 214110
  * Guo et al. 2018 Guo, S.; Li, Z.; Chan, G. K.-L. A Perturbative Density Matrix Renormalization Group Algorithm for Large Active Spaces. _J. Chem. Theory Comput._ 2018, _14_ , 4063–4071
  * Li et al. 2012 Li, Z.; Xiao, Y.; Liu, W. On the spin separation of algebraic two-component relativistic Hamiltonians. _J. Chem. Phys._ 2012, _137_ , 154114
  * Li et al. 2014 Li, Z.; Xiao, Y.; Liu, W. On the spin separation of algebraic two-component relativistic Hamiltonians: Molecular properties. _J. Chem. Phys._ 2014, _141_ , 054111
  * Balabanov and Peterson 2005 Balabanov, N. B.; Peterson, K. A. Systematically convergent basis sets for transition metals. I. All-electron correlation consistent basis sets for the 3d elements Sc–Zn. _J. Chem. Phys._ 2005, _123_ , 064107
  * 107 Li, Z. _(priviate communication)._