Semistochastic Second-Order Perturbative Correction

• Semistochastic second-order perturbative correction is a computational approach that reformulates energy calculations into deterministic and stochastic parts.
• It employs importance and stratified sampling with dynamic promotion of determinants to reduce variance and memory usage while achieving rapid convergence.
• This method underpins advanced techniques like selected CI, MRPT, and QMC, enabling high-accuracy studies of large quantum many-body systems.

A semistochastic second-order perturbative correction refers to a class of algorithms for evaluating second-order perturbation energy corrections—typically in quantum many-body and quantum chemistry context—where the correction is reformulated as a sum or integral, the dominant contributions are computed deterministically, and the remainder is estimated by stochastic sampling. This methodology enables evaluation of exceedingly large perturbative spaces with controlled sampling error and reduced memory requirements, while delivering rapid convergence and, when desired, zero statistical error at costs matching deterministic schemes. The approach is central to modern selected configuration interaction+perturbation theory schemes, multireference perturbation theory, and related methods.

1. Theoretical Basis: Reformulation of Second-Order Energy Corrections

In most quantum chemistry applications, the second-order perturbative energy correction is given by an Epstein–Nesbet (EN) partitioning. For a variational reference wavefunction $|\Psi\rangle = \sum_{I\in D} c_I |I\rangle$ with energy $E^{(0)}$, and its orthogonal complement $A$, the second-order correction takes the form:

$$E^{(2)} = \sum_{\alpha\in A} \frac{|\langle\alpha|H|\Psi\rangle|^2}{E^{(0)} - H_{\alpha\alpha}}$$

where $A$ is typically the space of determinants (or configurations) external to $D$ but connected via $H$ (Garniron et al., 2017). This sum is often intractable due to the exponential size of $A$.

To optimize evaluation, one can partition $A$ into disjoint subsets $A_I$, each associated with a reference determinant $|I\rangle$. This allows a decomposition $E^{(2)} = \sum_{I=1}^{N_{\text{ref}}} e_I$ with

$$e_I \equiv \sum_{\alpha\in A_I} \frac{|\langle\alpha|H|\Psi_I\rangle|^2}{E^{(0)} - H_{\alpha\alpha}}$$

where $\Psi_I$ is a truncated reference wavefunction. This sum of elementary contributions $e_I$ provides the foundation for efficient deterministic and stochastic sampling.

2. Semistochastic Algorithmic Structure

The semistochastic workflow segregates the perturbative correction into deterministic and stochastic parts. At iteration $m$:

$$E^{(2)} = E_{\text{det}}^{(m)} + E_{\text{sto}}^{(m)}$$

where

• $$E_{\text{det}}^{(m)} = \sum_{I \in D_{\text{det}}^{(m)}} e_I$$ includes already-converged contributions,
• $$E_{\text{sto}}^{(m)} = \sum_{I \notin D_{\text{det}}^{(m)}} e_I$$ is sampled stochastically.

Reference determinants are dynamically promoted from the stochastic to the deterministic pool as their contributions converge or are fully sampled (Garniron et al., 2017, Sharma et al., 2016). The key ingredients include:

• Importance sampling: Determinants are sampled with probabilities proportional to $c_I^2$ or $|c_I|$.
• Stratified/comb sampling: Faster convergence is achieved by stratifying the sample space and sampling multiple points per batch.
• Lazy evaluation: Each $e_I$ is computed at most once and cached, reducing redundant work.

In SHCI and related methods, further structure is introduced:

• Multiple thresholds ($$\epsilon_2^{\text{det}} > \epsilon_2^{\text{psto}} > \epsilon_2$$) define deterministic, pseudo-stochastic, and stochastic regions, with separate evaluation strategies for each (Yao et al., 2020, Li et al., 2018).

3. Error Scaling, Memory and Computational Efficiency

A hallmark advantage is the error scaling:

• In a fully stochastic Monte Carlo estimator, the error decays as $\sim t^{-1/2}$ (with $t$ the computational time).
• In semistochastic algorithms, due to progressive deterministic accumulation, empirical error decay is polynomial: $\sim t^{-n}$, with $n \approx 3-4$ observed in practical systems (Garniron et al., 2017, Li et al., 2018).
• Dynamical promotion of high-variance contributions to the deterministic sum dramatically reduces sampling variance.

Memory consumption is minimal: only the deterministic pool and the currently sampled subset of stochastic contributions are stored. This scaling (often $O(\max(N_{\text{det}},N_{\text{sample}}))$) enables treatment of spaces otherwise requiring unfeasible storage.

The computational cost per reference determinant $e_I$ is $O(|A_I|)$, and as $|c_I|$ or $c_I^2$ decreases, $|A_I|$ typically shrinks rapidly in selected CI or MRPT contexts; selection is therefore biased towards high-contribution determinants for maximal efficiency (Garniron et al., 2017, Sharma et al., 2016).

4. Algorithmic Implementation and Parallelization

Efficient implementation depends on the following topics:

• Sampling and Data Structures: High-performance hash tables index reference determinants, accumulate partial sums, or store cached energies (Li et al., 2018). Sampling is based either on sorted lists, CDF arrays (for Alias method or bisection), or direct heat-bath criterion (Tran et al., 2023).
• Batching and Hash-based Partitioning: Batches (disjoint segments of the external space) are assigned via determinant hashes to optimize memory usage and support parallel sample collection (Yao et al., 2020, Li et al., 2018).
• Parallelization: The process is naturally parallel:
  • Each thread or MPI rank draws its own random samples or works on disjoint batches;
  • Minimal synchronization is necessary except for global reduction of energies and error estimates;
  • Demonstrated strong scaling efficiencies of $\sim92\%$ (4096-way parallel) have been reported for MC-based stochastic perturbation theory (Cruz et al., 2022).
• Error estimation and dynamic control: The standard error is monitored—computed from fluctuations among MC batches or block averages; sampling proceeds until the desired tolerance is reached (Garniron et al., 2017, Tran et al., 2023).

5. Applicability to Selected CI, MRPT, QMC, and Vibrational Problems

Semistochastic second-order perturbative corrections underpin a variety of methodologies:

• Selected CI+PT (SHCI, HCI, VHCI): High-accuracy selected CI methods now routinely use semistochastic schemes for the PT2 correction, allowing for active spaces in excess of $10^9$ determinants with sub-$10\,\mu$Ha error (Li et al., 2018, Yao et al., 2020, Tran et al., 2023).
• Multireference PT2 (MRPT2): Decomposition of $E^{(2)}$ into reference-determinant sums enables fast, memory-efficient evaluation in multireference correlated treatments (Garniron et al., 2017).
• Stochastic Many-Body PT (MC-MP2): Real-space and imaginary-time MC integration, as in stochastic MP2, achieves quadratic to cubic scaling for relativistic systems, bypassing $O(N^5)$ bottlenecks and enabling large-core, heavy-element calculations (Cruz et al., 2022).
• Quantum Monte Carlo (i-FCIQMC+PT2): In i-FCIQMC, perturbative corrections are computed from discarded walker spawnings and stochastically accumulated with replica sampling, correcting for the truncation bias at minimal cost (Blunt, 2018).
• Vibrational Structure (VHCI+SPT2/SSPT2): The semistochastic paradigm extends directly to vibrational Hamiltonians, where exact evaluation is replaced by selection/sampling over product states in a sum-of-products basis (Tran et al., 2023).

6. Numerical Performance and Benchmarks

Quantitative results illustrate the efficiency and scalability:

Application	# Reference Determinants	System	Wall Time (PT2)	Statistical Error	Reference
F$_2$ (14e,108o,cc-pVQZ)	$1.05 \times 10^6$	F$_2$	$22\,\text{s}$	$0.9\,\text{mHa}$	(Sharma et al., 2016)
Cr$_2$ (12e,190o,cc-pVQZ)	$9.5 \times 10^6$	Cr$_2$	$56\,\text{min}$	$0.7\,\text{mHa}$	(Sharma et al., 2016)
Cr$_2$ (28e,176o,cc-pVQZ)	$2\times 10^7$	Cr$_2$	$2.5$–$18.5$ h	$\lesssim 0.1\,\text{mHa}$	(Garniron et al., 2017)
Au$_2$ (MP2)	—	Au$_2$	$\lesssim 1$h	$<10^{-3}$ rel. error	(Cruz et al., 2022)
Copper (SHCI)	$19 \times 10^6$ (V)	Cu	$3.7$–$4.2$ h	$0.9 \times 10^{-5}$ Ha	(Li et al., 2018)

In i-FCIQMC+PT2, the correction removes 80–99% of truncation error, often reducing statistical uncertainty below $1\,\text{mHa}$ at essentially no additional spawning cost (Blunt, 2018). For vibrational SPT2/SSPT2 in VHCI, statistical errors $<1$ cm$^{-1}$ are maintained for systems with $10^5$–$10^6$ configurations at 10$\times$ reduced memory relative to deterministic PT2 (Tran et al., 2023).

7. Practical Variations, Extension, and Best Practices

Key recommendations and best practices have emerged:

• Promotion strategy: Dynamically assign configurations to deterministic or stochastic pools based on contribution thresholds or convergence of MC samples (Garniron et al., 2017, Li et al., 2018).
• Batch and segment size: Optimal batch/segment count for stratified sampling balances variance reduction with bookkeeping cost (typical: $M=100-200$) (Garniron et al., 2017).
• Screening: Tight screening in both selection and PT2 stages maximizes efficiency; stochastic sampling fills in neglected small elements (Sharma et al., 2016, Yao et al., 2020).
• Efficient data structures: Use of distributed hash tables, sorted lists, and compressed bitword storage are standard (Li et al., 2018).
• Parallelization: Embarrassingly parallel algorithms maximize multicore and multinode efficiency, often requiring only a global reduction step (Li et al., 2018, Cruz et al., 2022, Tran et al., 2023).
• Error control: User-specified stochastic error tolerances (energy, or properties) direct runtime adaptivity (Li et al., 2018, Sharma et al., 2016).

A plausible implication is that future extensions may incorporate semistochastic frameworks into coupled cluster, density matrix renormalization group (DMRG) perturbation, and explicitly correlated methods, as the principle is agnostic to the particular reference space or Hamiltonian structure.

• "Hybrid stochastic-deterministic calculation of the second-order perturbative contribution of multireference perturbation theory" (Garniron et al., 2017)
• "Fast Semistochastic Heat-Bath Configuration Interaction" (Li et al., 2018)
• "Semistochastic Heat-bath Configuration Interaction method: selected configuration interaction with semistochastic perturbation theory" (Sharma et al., 2016)
• "Almost exact energies for the Gaussian-2 set with the semistochastic heat-bath configuration interaction method" (Yao et al., 2020)
• "An efficient and accurate perturbative correction to initiator full configuration interaction quantum Monte Carlo" (Blunt, 2018)
• "Stochastic evaluation of four-component relativistic second-order many-body perturbation energies: A potentially quadratic-scaling correlation method" (Cruz et al., 2022)
• "Vibrational heat-bath configuration interaction with semistochastic perturbation theory using harmonic oscillator or VSCF modals" (Tran et al., 2023)