Source-linked AI summary
Semistochastic Heat-bath Configuration Interaction method: selected configuration interaction with semistochastic perturbation theory
Sandeep Sharma, Adam Holmes, Guillaume Jeanmairet, Ali Alavi, C. J. Umrigar
TL;DR
Large active-space multireference calculations are constrained by determinant growth and the memory required by deterministic HCI perturbative corrections. The paper introduces stochastic and semistochastic Epstein-Nesbet perturbation theory for HCI, using direct Alias sampling to remove this bottleneck. SHCI computes challenging systems with better than 1 mHa accuracy and reported wall times of 55 seconds, 37 seconds, and 56 minutes.
Problem
Selected configuration interaction perturbative calculations become memory-limited because the perturbative determinant space grows rapidly as the cutoff is reduced for greater accuracy.
Method
The paper introduces stochastic and semistochastic multireference Epstein-Nesbet perturbation theory for HCI, sampling determinants directly from the variational wavefunction with the Alias method.
Results
Better than 1 mHa accuracy was obtained for challenging large-active-space systems including Mn-Salen (28e, 22o) and Cr2 (12e, 190o).
Takeaways & Limitations
SHCI extends HCI to very large active spaces while completely removing the deterministic perturbative memory bottleneck.
Abstract
from arXiv · showhide
We extend the recently proposed heat-bath configuration interaction (HCI) method [Holmes, Tubman, Umrigar, J. Chem. Theory Comput. 12, 3674 (2016)], by introducing a semistochastic algorithm for performing multireference Epstein-Nesbet perturbation theory, in order to completely eliminate the severe memory bottleneck of the original method. The proposed algorithm has several attractive features. First, there is no sign problem that plagues several quantum Monte Carlo methods. Second, instead of using Metropolis-Hastings sampling, we use the Alias method to directly sample determinants from the reference wavefunction, thus avoiding correlations between consecutive samples. Third, in addition to removing the memory bottleneck, semistochastic HCI (SHCI) is faster than the deterministic variant for many systems if a stochastic error of 0.1 mHa is acceptable. Fourth, within the SHCI algorithm one can trade memory for a modest increase in computer time. Fifth, the perturbative calculation is embarrassingly parallel. The SHCI algorithm extends the range of applicability of the original algorithm, allowing us to calculate the correlation energy of very large active spaces. We demonstrate this by performing calculations on several first row dimers including F2 with an active space of (14e, 108o), Mn-Salen cluster with an active space of (28e, 22o), and Cr2 dimer with up to a quadruple-zeta basis set with an active space of (12e, 190o). For these systems we were able to obtain better than 1 mHa accuracy with a wall time of merely 55 seconds, 37 seconds, and 56 minutes on 1, 1, and 4 nodes, respectively.
I. INTRODUCTION
Multireference methods address strongly correlated systems but are limited by the exponential growth of active-space calculations and perturbative memory demands. SHCI improves selected configuration interaction efficiency, while its semistochastic extension removes HCI’s memory bottleneck.
- CCSD(T) and related single-reference methods can fail catastrophically for multireference systems such as reactions and transition-metal compounds.
- CAS-based reference spaces grow exponentially with active orbitals, limiting the active-space size available to multireference correlation methods.
- Selected configuration interaction reduces CAS size by retaining important determinants and can add multireference Epstein-Nesbet perturbative corrections.
- HCI accelerates both variational and perturbative stages but originally stores all perturbative determinants, creating a memory bottleneck.
- SHCI introduces semistochastic multireference Epstein-Nesbet perturbation theory to overcome HCI’s memory bottleneck.
- SHCI samples determinants directly with the Alias method, avoids correlations between samples, and can trade memory for modest additional computer time.
B. Perturbative Stage
The perturbative stage uses the variational wavefunction to define the zeroth-order Hamiltonian and perturbation, while HCI truncates small contributions to control computational cost. Deterministic evaluation can nevertheless require storing enormous partial sums, creating a severe memory bottleneck that stochastic or semistochastic perturbation theory is designed to eliminate.
- The variational wavefunction defines the zeroth-order Hamiltonian H0 and perturbation V.
- The second-order energy expression is expensive because it requires summing many small terms.
- HCI discards terms smaller in magnitude than ϵ2, retaining only contributions with |Haici| > ϵ2.The threshold ϵ2 is kept much smaller than ϵ1 because discarding small-amplitude determinants can significantly affect dynamical correlation.
- 10^12 determinants and over 10 terabytes of memory arise for n = 12, v = 50, and Nv = 10^7 in the perturbative space.The original HCI algorithm reduces this storage requirement by orders of magnitude, but deterministic perturbative evaluation still requires storing partial sums.
- Stochastic or semistochastic perturbation theory can completely eliminate the memory bottleneck without increasing ϵ2.
III. STOCHASTIC MULTIREFERENCE PERTURBATION THEORY
The paper next introduces stochastic computation of the perturbative correction before presenting the more efficient semistochastic method. A figure illustrates that the perturbative-space determinant count grows rapidly as ϵ2 decreases.
- The stochastic method for computing the perturbative correction is discussed before the semistochastic method.
- The number of perturbative-space determinants increases rapidly as ϵ2 is reduced.For the C2 dimer with a QZ basis, the calculations used ϵ1 = 2 × 10−4 Ha and 403071 variational determinants; the results section uses ϵ2 = 10−8 Ha.
A. Stochastic PT
The stochastic perturbative method estimates the second-order correction by sampling variational determinants directly, producing unbiased estimates without autocorrelation between batches. Its cost and precision depend on the sample size, with stochastic computation becoming advantageous for large variational spaces.
- A. Stochastic PT: The perturbative correction is a bilinear function of the zeroth-order wavefunction coefficients.
- A. Stochastic PT: For any Nd ≥2, averaging over Ns samples gives an unbiased estimate whose precision improves progressively with more samples.
- A. Stochastic PT: Independent batches contain independently chosen determinants, eliminating autocorrelation between consecutive batches.
- A. Stochastic PT: CPU time per sample increases nearly linearly with Nd, with an additional Nd log(Nd) contribution.
- A. Stochastic PT: Beyond about Nd = 200, increasing the sampled determinant count produces only a much shallower decrease in the time needed to reach 0.1 mHa standard deviation.
- A. Stochastic PT: For large Nv and small ϵ2, the stochastic method requires fewer samples than the deterministic calculation and can therefore be more efficient despite stochastic error.
B. Semistochastic PT
Semistochastic perturbation theory combines a deterministic calculation at a loose threshold with stochastic correction of the omitted contribution. This interpolates between deterministic and stochastic algorithms while allowing memory, time, and statistical error to be traded.
- B. Semistochastic PT: The semistochastic method splits the perturbative calculation into deterministic and stochastic steps to reduce the sampling time required for large variational spaces and basis sets.
- B. Semistochastic PT: The stochastic correction estimates the difference between tight- and loose-threshold second-order energies and adds it to the deterministic loose-threshold result.
- B. Semistochastic PT: Using the same sampled determinants for both stochastic calculations substantially cancels stochastic error with almost no additional memory or computer time.
- B. Semistochastic PT: The loose threshold ϵd2 changes statistical error at fixed computer time but does not change the expected energy.
- B. Semistochastic PT: The method provides a smooth transition from the fully deterministic to the fully stochastic algorithm.
- B. Semistochastic PT: Figure 2 reports near-linear CPU-time-per-batch scaling with Nd and a rapid initial reduction in time to reach 0.1 mHa standard deviation.
IV. IMPLEMENTATION
The implementation combines determinant selection, sparse Hamiltonian construction, diagonalization, Alias sampling, and hybrid OMP/MPI parallelization. Its computational design emphasizes efficient variational processing and memory-conscious stochastic perturbation.
- IV. IMPLEMENTATION: The variational stage identifies significant determinants, builds the Hamiltonian matrix, and diagonalizes it.
- IV. IMPLEMENTATION: Identifying important determinants costs O(kNv ln(Nv) + kNv ln(Np)), where k is the average number of qualifying Hamiltonian connections and Np = kNv.
- IV. IMPLEMENTATION: The current implementation stores nonzero Hamiltonian elements in memory using list-of-lists sparse storage.
- IV. IMPLEMENTATION: The stochastic perturbation step samples Nd determinants and evaluates connected determinants contributing to the perturbative correction.
- IV. IMPLEMENTATION: The Alias method has a one-time O(Nv) memory cost and an O(Nd) cost for each drawn sample.
- IV. IMPLEMENTATION: Hybrid OMP/MPI parallelization assigns one MPI process per node and multiple threads per node, replicating the variational wavefunction across nodes.
V. BENCHMARKS
The benchmarks cover first-row dimers across multiple basis sets, strongly correlated Cr2, and the open-shell Mn-Salen complex. These systems probe the method across weakly and strongly correlated molecular regimes.
- V. BENCHMARKS: The first-row dimer benchmarks include C2, N2, O2, NO, and F2 with cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets.
- V. BENCHMARKS: Cr2 is evaluated with cc-pVDZ, cc-pVTZ, and cc-pVQZ active spaces of (12e, 68o), (12e, 118o), and (12e, 190o), respectively.
- V. BENCHMARKS: The benchmark set includes Cr2 as a strongly correlated dimer for which most multireference methods use no more than the minimal active space.
- V. BENCHMARKS: Mn-Salen is tested as a strongly correlated open-shell inorganic molecule with nearly degenerate singlet and triplet ground states.
A. First row diatomics
SHCI calculations on first-row dimers address very large active spaces while maintaining sub-milliHartree accuracy, with stochastic perturbation reducing computational demands.
- A. First row diatomics: Initial variational iterations used a larger ϵ1 because determinant coefficients tend to be larger when the variational space contains few determinants.The threshold was lowered toward its final value over subsequent iterations.
- A. First row diatomics: The variational cutoff ϵ1 was reduced through successive values, with three iterations performed at each value for the cc-pVQZ example.The reported sequence was 10^-3, 5 × 10^-4, 3 × 10^-4, and 2 × 10^-4 Ha.
- A. First row diatomics: F2 used the largest active space, with 14 electrons in 108 orbitals and a Hilbert space exceeding 10^20 determinants.On one node, its energy converged to better than 1 mHa in less than 3 minutes.
- A. First row diatomics: The stochastic method required less memory than the deterministic algorithm and used less computer time where the deterministic calculation was feasible.This comparison concerns obtaining sub-milliHartree accuracy for the reported systems.
- A. First row diatomics: Using MP2 natural orbitals enabled approximately a factor-of-3 speedup for F2 across all three basis sets.The calculations could use a larger ϵ1 threshold with these orbitals.
B. Cr2 dimer
For challenging Cr2 calculations, SHCI combines stochastic and semistochastic perturbation with large active spaces, rapid energy convergence, and substantial computational savings.
- B. Cr2 dimer: Cr2 calculations included all virtual orbitals in active spaces using cc-pVDZ-DK, cc-pVTZ-DK, and cc-pVQZ-DK basis sets.The calculations used a 1.68 Å bond length, frozen cores, and second-order Douglas-Kroll-Hess relativistic effects.
- B. Cr2 dimer: 12th-order excitations occur in the variational wavefunction, so CI expansions truncated at doubles or quadruples are far from adequate.This conclusion was obtained with ϵ1 = 8 × 10^-5 Ha across the three basis sets.
- B. Cr2 dimer: Total energies including perturbative corrections converge rapidly as ϵ1 decreases, although the variational energies remain far from convergence.The table reports both variational and total energies for the three basis sets.
- B. Cr2 dimer: 33194 sec was required by semistochastic PT versus 82800 sec for stochastic PT at 0.1 mHa statistical error in the TZ example.The semistochastic time comprises 394 sec of deterministic work plus 32800 sec of stochastic work.
- B. Cr2 dimer: DZ and QZ energies at the smallest ϵ1 values were within 1 mHa of their extrapolated energies, whereas the TZ energies were not.TZ′ calculations using approximate natural orbitals from an SHCI wavefunction were much better converged.
- B. Cr2 dimer: For a given basis set, the CPU time needed to reach a fixed statistical error is relatively insensitive to variational-space size.The perturbative correction decreases as the variational space grows, while added determinants have relatively few connections meeting the ϵ2 threshold.
C. Mn-Salen
SHCI was applied to the strongly correlated Mn-Salen cluster using a 28-electron, 22-orbital active space and orbitals converged with DMRG-SCF.
- C. Mn-Salen: Mn-Salen is a strongly correlated molecule whose derivatives are used to catalyze enantioselective epoxidation of olefins.The mechanism of the catalysis reaction is stated to be unknown in the supplied passage.
- C. Mn-Salen: Table III compares DMRG and SHCI energies for the singlet and triplet Mn-Salen states.SHCI used ϵ1 = 2 × 10^-4 Ha, ϵ2 = 1 × 10^-8 Ha, and Nd = 200; its single-node wall time is also reported.
- C. Mn-Salen: The Mn-Salen calculations used a (28e, 22o) active space and converged orbitals obtained from prior DMRG-SCF calculations.The initial orbitals were taken from HOMO-13 through LUMO+7 Hartree-Fock orbitals before DMRG-SCF optimization.
VI. CONCLUSIONS
SHCI removes HCI’s perturbative memory bottleneck while retaining efficient correlation-energy calculations for very large active spaces. The method achieves sub-1 mHa accuracy for challenging multireference systems, though variational-space Hamiltonian storage remains the largest memory requirement.
- SHCI removes the perturbative memory bottleneck and is faster than the fully deterministic algorithm for most systems at 0.1 mHa stochastic noise.The semistochastic perturbative calculation avoids storing all contributing determinants.
- SHCI efficiently computes correlation energies for very large active spaces, including Mn-Salen (28e, 22o) and Cr2 (12e, 190o).
- Within the studied systems, correlation energies were accurate to within 1 mHa, with convergence independently checkable for Cr2 despite unavailable published values.
- After the perturbative bottleneck is removed, storing the Hamiltonian in the variational space becomes the largest memory requirement.
- A remaining development goal is obtaining the variational wavefunction without storing the Hamiltonian.