Semistochastic Heat-Bath CI (SHCI)

• SHCI is a selected CI+PT method that constructs compact variational spaces using a heat-bath criterion and semistochastic PT2 to approach FCI accuracy.
• Its semistochastic PT2 step efficiently overcomes memory bottlenecks by sampling the perturbative space and capturing both static and dynamic correlations.
• SHCI is applied to a range of systems from strongly correlated molecules to transition metals, serving as a benchmark for high-accuracy electronic structure calculations.

Semistochastic Heat-Bath Configuration Interaction (SHCI) is a selected configuration interaction plus perturbation theory method that combines a heat-bath criterion for building a compact variational determinant space with a semistochastic evaluation of the multireference Epstein–Nesbet second-order correction. It is the semistochastic extension of Heat-Bath Configuration Interaction (HCI), systematically approaches the Full Configuration Interaction (FCI) limit as its thresholds are tightened, and is widely used as a high-accuracy classical reference for strongly correlated electronic structure, often alongside the Density Matrix Renormalization Group (DMRG) (Holmes et al., 2016, Sharma et al., 2016, Yoo et al., 11 Jan 2026).

1. Origins and formal placement

SHCI belongs to the family of selected configuration interaction plus perturbation theory (SCI+PT) methods. In this family, one first constructs a compact variational expansion over a selected subset of determinants and then adds a multireference second-order perturbative correction from the enormous external space. HCI introduced the heat-bath selection rule as an efficient deterministic analogue of heat-bath sampling; SHCI retained that variational stage and replaced the fully deterministic perturbative step with a semistochastic algorithm designed to remove the severe memory bottleneck of the original method (Holmes et al., 2016, Sharma et al., 2016).

Within the broader SCI+PT landscape, HCI was reported to achieve sub-millihartree accuracy in variational spaces containing up to $10^{20}$ determinants, while SHCI was described as allowing millihartree accuracy for active spaces of over $100$ orbitals (Yoo et al., 11 Jan 2026). In later applications, SHCI was used for problems with effective FCI spaces as large as $10^{38}$ determinants, while still delivering FCI-quality energies in the chosen one-electron basis (Chien et al., 2018, Jerzyk et al., 19 Nov 2025). These results explain why SHCI is repeatedly treated as a near-exact reference method for multireference molecules, transition-metal systems, heavy atoms, and benchmark Hamiltonian suites (Yao et al., 2020, Bellonzi et al., 14 Aug 2025).

A central point in current practice is that SHCI is not merely an approximate solver of convenience. In several recent studies it is explicitly positioned as part of the classical “gold standard” for strongly correlated electrons, and as the baseline against which emerging quantum algorithms are judged (Yoo et al., 11 Jan 2026, Bellonzi et al., 14 Aug 2025).

2. Variational selected-CI stage and the heat-bath criterion

The variational SHCI wavefunction is a CI expansion over a selected determinant space $V$,

$|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$

with the variational energy obtained by diagonalizing the Hamiltonian in that truncated space (Yao et al., 2021). The practical question is therefore how to enlarge $V$ efficiently.

HCI and SHCI use a heat-bath selection rule. A determinant $D_a$ connected to the current variational space is added if there exists at least one $D_i \in V$ such that

$|H_{ai} c_i| \ge \epsilon_1,$

where $H_{ai}=\langle D_a|\hat H|D_i\rangle$ and $100$0 is the variational selection threshold (Holmes et al., 2016, Yao et al., 2021). In another equivalent formulation, the importance function is

$100$1

and the determinant is selected if $100$2 (Holmes et al., 2016).

This criterion differs from CIPSI-style selection rules that use an approximate first-order coefficient or second-order energy estimate. The HCI criterion discards the denominator and replaces the full sum by the largest product $100$3. The data describe this as a major efficiency gain: because most nonzero matrix elements arise from double excitations and depend only on the four orbitals involved, SHCI can precompute and sort these couplings and then generate only those excitations that can pass the cutoff (Holmes et al., 2016, Smith et al., 2017). In consequence, the algorithm does not spend time generating obviously unimportant determinants.

The determinant search is therefore integral-driven rather than determinant-driven. For each pair of occupied spin-orbitals, SHCI walks through a pre-sorted list of double excitations and stops once the matrix element is too small to satisfy the threshold. Single excitations are treated similarly, though they are less dominant in cost (Holmes et al., 2016). As $100$4, the selected space approaches the full determinant space and the method reduces to FCI in the chosen basis (Jerzyk et al., 19 Nov 2025).

For excited states and state-averaged calculations, the same logic is generalized to multiple roots. One form used in excited-state SHCI is

$100$5

so that a common variational space contains determinants important for at least one targeted state (Chien et al., 2018). Relativistic and multistate variants also use state-averaged forms based on norms over several states rather than a single root, which is especially important when fine splittings or near-degeneracies are involved (Wang et al., 2023, Alaal et al., 2021).

3. Semistochastic Epstein–Nesbet perturbation theory

After diagonalization in the selected space, SHCI adds a second-order Epstein–Nesbet correction from determinants outside $100$6. With

$100$7

and variational energy $100$8, the standard SHCI correction is

$100$9

and the total energy is

$10^{38}$0

This expression appears throughout the SHCI literature and in later application papers as the canonical perturbative correction for the method (Sharma et al., 2016, Yoo et al., 11 Jan 2026, Jerzyk et al., 19 Nov 2025).

A first practical simplification is to screen the inner sum with a perturbative threshold $10^{38}$1, retaining only terms with $10^{38}$2 (Holmes et al., 2016, Yao et al., 2021). The remaining challenge is that even the screened perturbative space can be too large to store explicitly. The original deterministic HCI PT2 step therefore became memory-limited for realistic $10^{38}$3 values, which directly motivated SHCI (Sharma et al., 2016).

SHCI resolves this by sampling determinants from the variational wavefunction with probability

$10^{38}$4

using the Alias method rather than a Metropolis–Hastings random walk (Sharma et al., 2016, Jerzyk et al., 19 Nov 2025). The perturbative energy is then evaluated semistochastically: a deterministic contribution is computed with a looser threshold, and a stochastic correction estimates the remaining tail. In its simplest form,

$10^{38}$5

where $10^{38}$6, the superscript $10^{38}$7 denotes the deterministic part, and the two stochastic terms are evaluated with the same samples so that their fluctuations largely cancel (Sharma et al., 2016).

Several practical properties follow directly from this construction. The perturbative calculation is embarrassingly parallel; there is no sign problem; independent Alias-method sampling avoids autocorrelation issues associated with Markov-chain approaches; and memory can be traded against computer time by varying $10^{38}$8, $10^{38}$9, and the sampling effort (Sharma et al., 2016). For many systems, if a stochastic error of $V$0 mHa is acceptable, semistochastic PT2 is faster than the deterministic variant (Sharma et al., 2016).

Modern large-scale SHCI workflows refine this further with deterministic, pseudo-stochastic, and fully stochastic PT2 stages controlled by $V$1, $V$2, and $V$3 (Yao et al., 2020, Yao et al., 2021). This suggests a mature implementation pattern: the largest contributions are handled exactly, intermediate contributions are estimated from sampled perturbative batches, and the smallest contributions are sampled over both variational and perturbative spaces.

4. Orbital optimization, extrapolation, and implementation practice

At finite $V$4, SHCI is not invariant under unitary orbital rotations. This has a major practical consequence: better orbitals can make the selected expansion dramatically more compact, reduce the PT2 correction, and improve extrapolation to the FCI limit (Smith et al., 2017, Yao et al., 2021).

Orbital optimization in SHCI resembles CASSCF in that both CI coefficients and orbitals are optimized, but it differs in a fundamental way: in SHCI there is no predefined inactive/active/virtual partition, so effectively all orbitals are active and most rotations are nonredundant (Yao et al., 2021). The paper on orbital optimization in selected CI methods therefore distinguishes uncoupled, fully coupled, and quasi-fully coupled optimization strategies, and finds that taking CI–orbital coupling into account is crucial for fast convergence (Yao et al., 2021).

Two quasi-fully coupled methods are specifically recommended for SHCI applications: accelerated diagonal Newton and BFGS (Yao et al., 2021). Starting from natural orbitals, these methods were shown to converge much faster than uncoupled orbital updates while avoiding the cost and nonconvexity problems of fully coupled Newton steps (Yao et al., 2021). In parallel, HCISCF work established that active-space orbitals obtained from HCI-based optimization can be relatively insensitive to the accuracy of the underlying HCI solve, making cheap orbital optimization followed by tighter SHCI energetics a practical workflow (Smith et al., 2017).

Extrapolation is another standard part of SHCI methodology. In many applications, one computes $V$5 at several $V$6 values and fits the total energy against the residual perturbative correction. Weighted quadratic extrapolation in $V$7 is used in thermochemistry and benchmark studies (Yao et al., 2020), while linear or quadratic extrapolations in related variables are used in excited-state, atomic, and multicomponent settings (Chien et al., 2018, Jerzyk et al., 19 Nov 2025, Alaal et al., 2021). The common idea is that $V$8 acts as a proxy for the remaining correlation missing from the variational space.

Implementation-wise, SHCI has been used through the Arrow-SHCI code and closely related selected-CI infrastructures. In one HCISCF example, an Fe-porphyrin model complex with an active space of $V$9 required 412 seconds per iteration on a single node containing 28 cores, of which 185 seconds were spent in the HCI calculation and 227 seconds mainly in integral transformation (Smith et al., 2017). This suggests that, once the perturbative bottleneck is controlled, integral transformation and sparse Hamiltonian handling can dominate the wall time.

5. Benchmarking performance and scientific applications

SHCI has been applied across strongly correlated molecules, excited states, thermochemistry, transition-metal chemistry, heavy atoms, and bioinorganic benchmark systems. The following representative cases illustrate the method’s range.

Domain	Representative result	Paper
Large active-space correlation	Better than $|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$0 mHa accuracy for F$|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$1 $|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$2, Mn–Salen $|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$3, and Cr$|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$4 $|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$5 with wall times of 55 s, 37 s, and 56 min	(Sharma et al., 2016)
Excited states	Hexatriene in ANO-L-pVDZ, with Hilbert space $|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$6, gave $|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$7 and $|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$8 eV for $|\Psi_V\rangle = \sum_{D_i \in V} c_i |D_i\rangle,$9 and $V$0; ozone ring minimum lies $V$1 eV above the open-ring minimum and is separated by a $V$2 eV barrier	(Chien et al., 2018)
Thermochemistry	CBS-extrapolated atomization energies for the G2 set gave MADs of $V$3 kcal/mol without and $V$4 kcal/mol with the basis-set correction	(Yao et al., 2020)
Transition metals and monoxides	SHCI plus density-based basis-set correction was reported to converge total, ionization, and dissociation energies to the CBS limit within chemical accuracy	(Yao et al., 2021)
Heavy atoms	First ionization potentials of Cr, Mo, and W were reported as $V$5, $V$6, and $V$7 eV	(Jerzyk et al., 19 Nov 2025)
Quantum-algorithm benchmarks	For $V$8 and $V$9, HCI/SHCI serve as the classical reference, with HCI selected spaces reaching tens to hundreds of millions of determinants	(Yoo et al., 11 Jan 2026)

In strongly correlated bond breaking and transition-metal clusters, SHCI is used because it recovers both static and dynamic correlation in a systematically improvable way. In the $D_a$0 dissociation benchmark with an all-electron $D_a$1 active space and a CAS size of $D_a$2 determinants, HCI with $D_a$3 was treated as essentially FCI-quality, while in the $D_a$4 cluster with a $D_a$5 active space and a CAS size of $D_a$6 determinants, tighter HCI thresholds brought the energy within $D_a$7 mHa of CASCI (Yoo et al., 11 Jan 2026).

In thermochemistry, the G2-set study demonstrates the distinction between “exact within basis” and the CBS limit. SHCI makes the former essentially available, and then basis extrapolation or basis-set correction becomes the dominant remaining issue (Yao et al., 2020). In transition-metal atoms, ions, and monoxides, this division of labor is even more explicit: SHCI provides near-FCI energies in each basis, while density-based basis corrections accelerate convergence to the CBS limit (Yao et al., 2021).

The method is also noteworthy as a benchmark generator. In the QB ground-state energy benchmark, “optimized SHCI” was the only solver reported to attain a solvability score of $D_a$8, with 148 tasks solved out of 226 attempted, and SHCI was used directly or indirectly to define FCI-quality references across the dataset (Bellonzi et al., 14 Aug 2025).

6. Extensions, benchmarking role, and limits

SHCI has been extended well beyond nonrelativistic ground-state molecular CI. A relativistic formulation for arbitrary two-component and four-component Hamiltonians was presented for systems such as $D_a$9 and $D_i \in V$0, correlating more than $D_i \in V$1 spinors in both cases (Wang et al., 2023). Earlier one-step SOC treatments had already shown that HCI and its semistochastic extension could treat spin–orbit coupling and correlation on an equal footing in large active spaces; for the Au atom, converged excitation energies were obtained with just over $D_i \in V$2 determinants in an active space containing $D_i \in V$3, whose full determinant count exceeds $D_i \in V$4 (Mussard et al., 2017).

Atomic-structure applications likewise broadened the scope of SHCI. In calculations of first ionization potentials of Cr, Mo, and W, SHCI was combined with orbital optimization, effective core potentials, $D_i \in V$5-based extrapolation to the FCI-in-basis limit, and basis-set extrapolation to the CBS limit (Jerzyk et al., 19 Nov 2025). The same paper emphasizes a practical parameterization in which

$D_i \in V$6

linking perturbative thresholds to the variational threshold (Jerzyk et al., 19 Nov 2025).

Methodological generalizations show that the underlying SHCI logic is portable. Multicomponent HCI for protonic excited states introduced deterministic Epstein–Nesbet PT2 in a product space of electronic and protonic determinants and explicitly identified a semistochastic extension as natural but not yet implemented there (Alaal et al., 2021). Vibrational heat-bath CI later imported the same semistochastic PT2 machinery into vibrational structure theory, reporting stochastic errors controllable to less than $D_i \in V$7 cm$D_i \in V$8 (Tran et al., 2023).

The method’s present limits are also clear. The original SHCI paper explicitly states that, after removing the PT2 memory bottleneck, the variational Hamiltonian becomes the dominant memory object (Sharma et al., 2016). Application papers and benchmarks likewise emphasize memory constraints and runtime blow-ups for very large active spaces, very strong correlation, and large perturbative spaces (Bellonzi et al., 14 Aug 2025, Jerzyk et al., 19 Nov 2025). In one comparative discussion, CI diagonalization cost is stated to scale roughly as $D_i \in V$9 in the number of determinants $|H_{ai} c_i| \ge \epsilon_1,$0, so reductions in selected-space size are directly consequential (Yoo et al., 11 Jan 2026).

A further limitation is methodological rather than algorithmic: benchmark universality can be misleading. The QB benchmark reports near-universal solvability for optimized SHCI on its current dataset, but also states that many Hamiltonians originate from prior SHCI-oriented studies, introducing a bias that favors classical selected-CI solvers (Bellonzi et al., 14 Aug 2025). This means that SHCI’s apparent universality on such datasets should not be read as universal dominance across all chemically relevant or physically difficult Hamiltonians. Systems with extremely large active spaces, metallic behavior, long-range entanglement, or especially unfavorable determinant connectivity remain plausible hard cases (Bellonzi et al., 14 Aug 2025).

Taken together, these developments define SHCI as a mature, systematically improvable SCI+PT framework: fast determinant selection by the heat-bath criterion, semistochastic multireference Epstein–Nesbet PT2, strong synergy with orbital optimization and extrapolation, and a demonstrated ability to deliver near-FCI reference data across electronic, relativistic, atomic, multicomponent, and vibrational settings (Sharma et al., 2016, Yao et al., 2021, Wang et al., 2023).