Numerically exact configuration interaction at quadrillion-determinant scale

Domain: physics
Published: 2025-12-10
Editorial outcome: major_revision

Abstract
The combinatorial growth of configuration interaction (CI) has long limited this formally exact quantum chemistry method to only the smallest molecules. Here, we report a numerically exact CI calculation exceeding one quadrillion (10^15) determinants, made possible by a lossless categorical compression strategy within the small-tensor-product distributed active space (STP-DAS) framework. This approach overcomes the traditional memory bottlenecks of CI by a numerically exact compression of the wavefunction representation and reformulating the most computationally demanding matrix–vector operations. Using this method, we performed a fully relativistic CI calculation of the ground state of HBrTe with over 10^15 complex-valued determinants in just 34.5 h on 1000 computing nodes—the largest CI calculation ever reported. We further achieved fast computation for systems with hundreds of billions of determinants on only a few compute nodes. Extensive benchmarks confirm that the method retains full numerical exactness while cutting memory and computational cost by orders of magnitude. Compared to previous state-of-the-art CI calculations, this work achieves a 1000 times increase in CI space, a 10^6-fold increase in floating-point operations performed, and a 10^6-fold improvement in computational speed.

Introduction
Full configuration interaction (FCI) offers the most complete and accurate description of a molecule’s electronic structure within a given basis set, providing the exact spectral solution to the non-relativistic electronic Schrödinger equation 1–6. Due to its variational nature, FCI is particularly well-suited for treating relativistic effects, such as spin–orbit and spin–spin couplings beyond perturbation theory 7–23, which are fundamentally rooted in the electronic Dirac equation.

At its core, solving the FCI problem reduces to diagonalizing a large many-electron Hamiltonian matrix. This matrix is Hermitian, sparse, and typically diagonally dominant, making it well-suited to iterative diagonalization techniques that can efficiently converge on a few eigenstates without requiring full storage or construction of the entire matrix. However, as the number of determinants grows factorially with system size, reflecting the combinatorial nature of Slater determinant enumeration within the full Hilbert space, even iterative methods become intractable beyond a certain threshold.

The CI wavefunction is often expressed as a linear combination of Slater determinants, typically generated by excitations from a mean-field self-consistent field (SCF) ground-state reference. In the relativistic regime, this framework must be reformulated using complex-valued 2- or 4-spinor wavefunctions, as required by the Dirac formalism 24,25.

Figure 1 illustrates the historical progression of CI implementations, highlighting major breakthroughs in the achievable scale of determinants. Prior to this work, over a span of 35 years (1990–2024), the field advanced from handling billions to trillions of determinants, driven largely by advances in computer hardware technologies. Although CI is amenable to large-scale parallel processing schemes 26–31, the explosive growth in memory requirements has historically restricted its applicability to only the smallest chemical systems. Relativistic CI is even more limited, due to the intrinsically larger spinor configuration space associated with complex-valued 2- or 4-component wavefunctions. Simply put, enabling CI for practical quantum chemistry applications demands alternative theoretical frameworks and data representations that can circumvent the brute-force enumeration of the CI space.

Many CI-based wavefunction methods aim to approximate the FCI solution. These methods use different types of approximations, which affect the accuracy of the resulting wavefunction. The two main approaches are complete active space CI (CASCI) and selected CI (SCI). While both CASCI and SCI methods effectively truncate the Hilbert space of the system to a subspace of significant determinants, this significance is determined differently and at different stages of the computation.

In the CASCI method 9,10,29,30,32–38, it is assumed that only a subspace of the full Hilbert space of the system contains meaningful correlation, and the FCI wavefunction is approximated as the CI wavefunction in the truncated space (the so-called active space). The truncation often leads to an underestimation of dynamic correlation. Applying more computationally demanding methods such as multiconfigurational self-consistent field (MCSCF) 6–8,12,13,16,39–48, multireference configuration interaction (MRCI) 6,19,48–54, and many-body perturbation theory (MRPT2, CASPT2, NEVPT2, MC-PDFT) 13,14,22,55–61 is typically required to achieve qualitative and quantitative agreement with experiment.

SCI-based methods estimate the importance of each configuration in the total wavefunction based on a predefined significance criterion, which depends on the chemical problem of interest 62–78. Once the most significant determinants are identified, the Hamiltonian is constructed and diagonalized within this reduced space to approximate the FCI wavefunction. As in CASCI, SCI approaches are often combined with perturbation theory to recover contributions from neglected configurations and improve quantitative accuracy 64,65,67–69.

Even with advances in dimensionality reduction, CI remains fundamentally constrained by memory limitations. State-of-the-art implementations still require explicit storage of either Hamiltonian matrix elements or excitation lists to support on-the-fly matrix-vector operations, commonly referred to as the \sigma-build 32,79–81. For large CI spaces, storing the full Hamiltonian matrix and performing direct diagonalization is clearly impractical. In on-the-fly CI algorithms, the one-electron excitation list, which encodes the allowed excitations between determinants for efficient Hamiltonian construction, scales as n_e \times (n_h + 1) \times N, where N is the number of determinants, and n_e and n_h denote the numbers of electrons and holes (unoccupied orbitals), respectively. Since N increases factorially with system size, the associated memory requirements grow rapidly, making conventional CI calculations infeasible for anything beyond the smallest systems.

For a CI problem involving N determinants, the size of each CI expansion vector scales linearly with N. For instance, a relativistic CI problem with one quadrillion (10^15, 100 orbitals, 88 electrons) determinants would require ~16 petabytes (PB) of memory just to store a single CI vector composed of complex-valued double-precision coefficients. While this memory footprint alone makes such problems challenging to tackle, the memory required to store the excitation list can easily scale to the exabyte (EB) regime. This poses a fundamental barrier to scalability, even before considering computational cost.

A recently introduced CI matrix-vector product algorithm 31 leverages the exact factorization of the active space into small tensor products of distributed active spaces—an approach known as the small-tensor product distributed active space (STP-DAS) framework, illustrated in Fig. 2A. The STP-DAS algorithm reformulates the large CI matrix-vector product into a sequence of small tensor products, each embedded within a distributed active space, computed on-the-fly using string-based methods. This advance exploits the mathematical condition governing the phase relationship between the global address and the local DAS address of any CI matrix element, enabling the use of only small local determinant address strings in the CI matrix-vector build and overcoming the memory bottleneck associated with storing the full excitation list. This formulation enables extensive reuse of Hamiltonian excitation lists, leading to a dramatic reduction in memory demands. For a CI problem involving one quadrillion (10^15) determinants, the STP-DAS framework reduces the excitation list memory requirement from 12 exabytes (EB) to 25 gigabytes (GB), an 8-orders-of-magnitude reduction! In addition, by evenly distributing the computation of small tensor products, the STP-DAS algorithm achieves excellent load balance with minimal node-to-node communication overhead, ensuring strong scalability across both single-node and large-scale parallel architectures.

Since the excitation list typically dominates the storage requirements in CI calculations, the STP-DAS framework overcomes a long-standing memory bottleneck and enables CI computations that were previously deemed intractable. By effectively eliminating the memory footprint of the excitation lists, the storage of CI vector coefficients now emerges as the bottleneck in many large-scale CI problems.

Revisiting the earlier CI example of 10^15 determinants: even after eliminating the memory bottleneck associated with storing excitation lists, the same calculation still demands 16 PB of memory to store the numerically exact coefficients of a single CI expansion vector in an iterative solver. Given that practical CI calculations typically require multiple subspace vectors for convergence, it becomes clear that storing these expansion vectors now represents the dominant memory bottleneck in large CI calculations. To address this challenge, we leverage the locally compressed nature of the STP-DAS framework to efficiently compute ultra-large-scale CI problems involving up to a quadrillion (10^15) determinants. This approach yields deterministic, numerically exact solutions and effectively shifts CI calculations from being memory-bound to compute-bound.

Results
Before presenting methodological details and performance benchmarks, we highlight the largest CASCI calculation to date for the ground state of HBrTe, enabled by the compression-compatible STP-DAS framework to be introduced herein. HBrTe is a substituted form of a hydrogen chalcogenide where one of the hydrogens was substituted with a bromine atom to decrease the symmetry to the molecule. We performed a relativistic exact two-component 21,22,82–93 CASCI (X2C-CASCI) calculation (100 2-spinor orbitals, 88 electrons, complex-valued 1.05 \times 10^15 2-spinor determinants) of the ground state of HBrTe. We also performed a calculation with the same number of determinants for the ground state of a magnesium atom (see Section S3 of the Supplementary Information). The calculation ran on the National Energy Research Scientific Computing Center’s Perlmutter high-performance supercomputer with a total of 1000 nodes (AMD EPYC 7763 Milan, 128,000 compute cores, 512 GB of RAM per node, 200 GB/s NIC, 2 message passing interface (MPI) processes per node, and 64 symmetric multiprocessing (SMP) threads per MPI process).

The excitation list is generated without assuming any symmetry of the target state. Consequently, the calculation is formally performed in a quadrillion-determinant space, with all determinants explicitly included. Because the CI coefficients are complex, the memory footprint is twice that of an analogous non-relativistic calculation. Moreover, the complex arithmetic makes the computational cost (in FLOP count) equivalent to that of a non-relativistic calculation with more than twice as many determinants.

Table 1 summarizes the results of the HBrTe calculation. The ground-state energy converged in 9 iterations to microhartree precision (<10^-6 a.u.) with a total runtime of 34.5 h. In each iteration, an additional CI expansion vector was introduced to accelerate convergence, while the compression algorithm dynamically adapted to the expanding vector space, leading to a gradual increase in computational cost. A total of 9 expansion vectors were involved in the \sigma-build. Despite the enormous configuration space of over one quadrillion (1.05 \times 10^15) complex-valued 2-spinor determinants, the average \sigma-build time per vector remained just 3.8 h.

The ground-state energy of the HBrTe molecule from this CI calculation is -9395.028004 a.u. Leveraging the gap theorem 94,95, we determine that our X2C-CASCI result lies within 10 \times 10^-6 a.u. of the true x2c-TZVPall96 ground state energy within that active space, well below any chemically meaningful threshold. A detailed analysis is provided in “Methods”.

This work represents a 3-orders-of-magnitude increase in CI space and a 6-orders-of-magnitude increase in FLOP count, which is estimated using O(N^2), compared to the previous state-of-the-art in CI calculations 30,31. Compared to previous state-of-the-art CI calculations, this work also achieves a 6-orders-of-magnitude speedup in time-to-solution. Such a large-scale CI calculation is enabled by the STP-DAS-based numerically exact categorical compression scheme, which reduces the memory required to store the 9 CI expansion vectors from 134 PB to less than 500 TB while maintaining a good load balance 31, making the computation feasible on most existing supercomputing infrastructures.

We now describe the algorithmic developments that enable such CI calculations to be performed on existing supercomputing resources. Detailed algorithms, parallel implementation strategies, and error bound analyses are provided in the Supplementary Information. The central concept of the STP-DAS framework is the systematic partitioning of the full CI orbital space into a collection of distributed active spaces. Within each active space, configurations are further classified into categorical subspaces, rigorously defined by distinct electron occupation patterns, as illustrated in Fig. 2A 31. The STP-DAS framework reformulates the CI \sigma-build as a sum of small-tensor products, each uniquely addressed via a global tensor looping structure. The major memory bottleneck associated with storing the excitation list is eliminated by allowing categorical subspaces to share compact, local excitation lists.

In the largest CASCI calculation (10^15 determinants) presented here, employing 13 distributed active spaces, the STP-DAS approach reduces the excitation list memory requirement from 12 \times 10^9 GB to just 25 GB. However, storing all nine CI expansion vectors for a system with 10^15 determinants would require ~134 PB of memory. On a high-performance computing system such as Perlmutter, this translates to more than 275,000 nodes, each equipped with 512 GB of memory. A straightforward element-wise sparsity treatment, however, does not meet the STP-DAS condition. To overcome this limitation, the following section introduces a categorical compression scheme that achieves the necessary memory reduction, making it possible to carry out large relativistic CI calculations on a medium-sized computing cluster.

With the STP-DAS framework, the memory bottleneck associated with storing CI expansion vectors is effectively eliminated by applying numerically exact categorical compression. The STP-DAS CI expansion vectors take the form

C = \sum_B C_B, (1)

where B is a category, defined by a unique electron occupation pattern within the distributed active spaces 31.

The compression scheme, categorical compression (see Fig. 2D), stores the C_B vectors as compressed sparse column (CSC) vectors. In this format, each categorical expansion vector C_B is represented by storing the values C_B \neq 0 and their corresponding local addresses. The tensor-loop structure in the STP-DAS \sigma-build algorithm is reformulated in terms of categorically compressed local addresses together with their corresponding global phase factors (see Fig. 2E).

In contrast to element-wise compression, categorical compression can eliminate all configurations within a category, i.e., skip an entire category at once. This approach is better able to preserve the vectorized structure compared to element-wise compression. The major difficulty of any sparse matrix-vector product algorithm is the lack of a priori knowledge of the location of the nonzero elements. This difficulty leads to a major bottleneck rooted in the nonuniform accesses to memory. The categorical compression localizes nonuniform memory accesses within a category, which is always orders-of-magnitude smaller than the size of the full CI space. This localized memory access pattern is the central strength of the categorical compression scheme.

Within each category, element-wise compression is still applicable to maximize sparsity. Most importantly, a category-based compression scheme naturally supports distributed small tensor products, taking advantage of the reduced memory footprint and improved parallel load balance of the STP-DAS framework.

In summary, the numerically exact categorical compression introduced here allows the STP-DAS \sigma-build algorithm to bypass entire categories of determinants while preserving both the vectorized structure and the local addressing scheme of STP-DAS. Because the compression is fully lossless, the omitted determinants have no impact on the resulting Ritz eigenvalue–eigenvector pair.

The final hurdle in reducing CI memory demands lies in the iterative solver. In CI, the Hamiltonian operator is diagonalized iteratively within the full Hilbert space of the system. As a result, the error in the computed Ritz value directly reflects the missing correlation energy in the associated approximate wavefunction. This enables the Ritz residual, defined as r \equiv \|HC - EC\|, to be computed and appended as an additional CI expansion vector. The norm of the residual provides rigorous bounds on the missing correlation energy relative to the true eigenpair (eigenvector and eigenvalue) of the Hamiltonian in the chosen basis 94,95,97–103. A widely used approach that leverages this principle is the Davidson iterative solver 104. Other related methods that use residual norm as the convergence criterion include the locally optimal block preconditioned conjugate gradient (LOBPCG) method 105, the Jacobi-Davidson 106 method and generalized preconditioned locally harmonic residual (GPLHR) method 107,108.

By applying numerical or convergence thresholds at various stages of the Davidson method 104, one can exploit the sparsity of newly generated CI expansion vectors. With sufficiently tight thresholds, these vectors can span the part of the Hilbert space required to accurately represent the desired wavefunction and drive the iterative diagonalization to any desired level of precision. The Davidson method 104 utilized the Davidson preconditioner, which generates the ith component of the next trial expansion vector according to

t_i = \begin{cases} 0, & \text{if } |\lambda - H_{ii}| \text{ is small} \\ \frac{r_i}{\lambda - H_{ii}}, & \text{else} \end{cases} (2)

Here, \lambda is the Ritz value of the current iteration, H_{ii} is the ith element of the diagonal of the Hamiltonian H, and r_i is the ith component of the current residual. Among the various preconditioners employed in the Davidson method 97,98,103,109–111, compression-compatible preconditioners, which discard terms t_i below a numerical threshold \epsilon, have been shown to achieve convergence to the exact same results as the traditional Davidson preconditioner 97–99,101,112.

In this work, we apply the compression-compatible categorical preconditioner

t_i = \begin{cases} \frac{r_i}{\lambda - H_{ii}}, & \text{if } |\lambda - H_{ii}| \geq 10^{-12} \text{ and } |s_i| \geq \epsilon \|s\| \\ 0, & \text{else} \end{cases} (3)

using the nonzero residual elements r_i. We do this on-the-fly to avoid explicitly storing the prohibitively large diagonal of H (a dense vector of size N_{dets}). Figure 2F illustrates the expansion of the CI vector space enabled by the compression-compatible categorical preconditioner used in the Davidson method implemented here. Eq. (3) closely resembles the preconditioner proposed in ref. 112, with the key distinction being the inclusion of the Davidson-preconditioned residual norm \|s\| in the dropping criterion. This facilitates dynamic threshold adjustment: as the iterations progress and \|s\| decreases, the criterion becomes more stringent. More importantly, the factor of \|s\| ensures that numerical thresholding is applied to the generated expansion vectors relative to the total norm of their exact (traditional Davidson) counterparts, rather than some fixed cutoff on the absolute values of their entries. This results in a very accurate CI space expansion scheme at the cost of computing and contracting a significant number of Hamiltonian matrix elements (to evaluate s exactly), which would be intractable without the STP-DAS framework.

Algorithms and pseudocodes of the compression-compatible STP-DAS method are presented in “Methods”, along with discussions on load balancing and parallel implementation. The convergence behavior of the STP-DAS framework equipped with the compression-compatible preconditioner defined in Eq. (3) was evaluated across three systems with varying degrees of electron correlation: the magnesium atom, diatomic nitrogen, and a model carbon nanotube. The results are provided in Section S1 of the Supplementary Information.

The analysis reveals that overly aggressive thresholding can cause the Davidson procedure to stagnate, thereby preventing convergence to the correct electronic wavefunction. When thresholds are too loose, newly generated trial vectors quickly become linearly dependent on the existing subspace vectors, signaling that the span of the modified subspace has saturated before achieving convergence. However, when the preconditioning threshold satisfies

\epsilon \lesssim \frac{1}{\sqrt{N_{dets}}}, (4)

the resulting energies agree with their exact values to better than 10^-7 E_h, and the residual norms become correspondingly small, demonstrating successful and reliable convergence.

The CI wavefunction of highly correlated systems is comprised of a large number of determinants with small CI coefficients. Because the Hilbert space of the problem is never truncated, and no determinants are discarded from the wavefunction itself. As a result, the compression-compatible preconditioner easily facilitates convergence to the exact wavefunction, provided that the preconditioning threshold \epsilon is sufficiently small for the true eigenvector to be accurately represented in the subspace spanned by the modified expansion vectors. With the capability to perform large CI calculations, energetic extrapolation to the correlation limit becomes feasible for many-electron systems. Figure 2B illustrates the correlation-consistent extrapolation for the Mg^2+ ion using double-zeta (DZ, 36 orbitals), triple-zeta (TZ, 68 orbitals), and quadruple-zeta (QZ, 118 orbitals) basis sets, involving 2.54 \times 10^8, 2.91 \times 10^11, and 9.75 \times 10^13 2-spinor determinants, respectively. The complete basis set limit of -199.16704295 a.u. was obtained using a mixed Gaussian extrapolation scheme tailored for correlation-consistent basis sets 113.

To demonstrate the strong scaling behavior of the compression-compatible STP-DAS \sigma-build algorithm, we performed relativistic X2C-CASCI calculations on the thallium hydride (TlH) molecule using active spaces of 40 and 41 2-spinor orbitals and 24 electrons, corresponding to 63 and 152 billion 2-spinor determinants, respectively. Five distributed active spaces (DASs) were employed. These calculations were executed on the University of Washington’s Hyak HPC system, a small-sized cluster where each node is equipped with two Intel Xeon 6230 Gold CPUs and a single 100 GB/s network interface card. As shown in Fig. 2C, even with just 5 compute nodes, relativistic X2C-CASCI calculations involving tens to hundreds of billions of determinants require only 4–6 min per \sigma-build on average. Increasing to 30 nodes further reduces the cost to just over 1 min. Past 30 nodes, the communication time dominates the runtime and the calculations no longer scale. This benchmark demonstrates that billion- and even trillion-determinant CI calculations are now feasible on a small-scale computing cluster.

Additionally, we performed X2C-CASCI calculations for the ground states of two highly correlated systems, showcasing the applicability of the compression-compatible STP-DAS framework to highly correlated systems. These systems were chosen from opposite ends of the correlation spectrum from strongly statically correlated to strongly dynamically correlated. Square Rb_4, a relativistic analog of H_4 114–126, displays strong static correlation and Xe_2 is a dynamically correlated noble gas dimer 127–135.

Tables 2 and 3 summarize the results of the Rb_4 and Xe_2 calculations. The Rb_4 calculation (50 2-spinor orbitals, 28 electrons, 8.9 \times 10^13 2-spinor determinants) ran on the National Energy Research Scientific Computing Center’s Perlmutter high-performance supercomputer with a total of 100 nodes (AMD EPYC 7763 Milan, 12,800 compute cores, 512 GB of RAM per node, 200 GB/s NIC, 1 MPI processes per node, and 128 SMP threads per MPI process) and took 6 iterations and 11.8 h to converge. The Xe_2 calculation (60 2-spinor orbitals, 12 electrons, 1.4 \times 10^12 2-spinor determinants) ran on the same platform with 256 nodes and took 7 iterations and 36.1 h to converge.

The ground-state energies of the Rb_4 and Xe_2 molecules from these calculations are -11916.152725 a.u. and -14889.646696 a.u., accordingly. The gap theorem 94,95 guarantees that these X2C-CASCI results lie within 0.53 \times 10^-6 a.u. and 17.56 \times 10^-6 a.u. of the true ground state energies within the corresponding active spaces and basis sets (see “Methods”).

In summary, by combining compression-compatible preconditioners with compression-compatible categorical CI vectors, the STP-DAS framework drastically reduces the memory footprint of both the excitation lists and the CI expansion vectors. In the largest relativistic CASCI calculation presented here, spanning 10^15 determinants across 13 distributed active spaces, the STP-DAS approach reduces the memory required for the excitation list from 12 \times 10^9 GB to just 25 GB, and for 9 CI expansion vectors from 134 PB to less than 500 TB. These reductions make quadrillion-determinant calculations tractable on current supercomputing architectures. While most of the community may not have access to the hundreds of compute nodes required for such runs, this work also demonstrates the practical feasibility of trillion-determinant calculations on just a few nodes and even on a laptop.

Methods
Lossless \sigma-build using the categorical compression of small tensor products

The categorical \sigma-build algorithm within STP-DAS 31 can be reformulated to exploit the categorical compression of the expansion vectors. The compact nature of the categorical representation enables fast and memory-efficient computation of \sigma-vectors. The categorical \sigma-build algorithm implements the evaluation of

\sigma^{LA} = ^1e\sigma^{LA} + ^2e\sigma^{LA}, (5)

^1e\sigma^{LA} = \sum_B \sum_{KB} \sum_{\mu \nu} P_{\mu \nu} \delta_{A} \delta_{B} h^0_{pq} h^{LA \mu \nu}_{KB \mu \nu} \langle LA \mu \nu | E^p_q | KB \mu \nu \rangle C_{KB}, (6)

^2e\sigma^{LA} = \sum_B \sum_{KB} \sum_{C} \sum_{JC} \delta_{A} \delta_{B} \delta_{C} \delta_{JC} P_{\mu \nu \kappa \lambda} g_{pqrs} h^{LA \mu \nu}_{KB \kappa \lambda} \langle LA \mu \nu | E^p_r E^q_s | JC \mu \nu \rangle C_{KB}, (7)

where p \in XA \mu, q \in XC \nu, r \in XC \kappa, s \in XB \lambda, \langle XA \mu \nu | E^p_q | XB \mu \nu \rangle are categorical one-electron excitation lists, P_{\mu \nu} are global phase factors, and h^0_{pq} (g_{pqrs}) are one (two) body Hamiltonian elements. See ref. 31 for algorithmic details.

Equations (5) to (7) define the categorical \sigma-vector in terms of local STP-DAS one-electron excitation lists and the categorical CI expansion vector. Notably, when the expansion coefficients C_{KB} are categorically compressed, the resulting \sigma-vector coefficients \sigma^{LA} are also categorically compressed. In such cases, the categorical \sigma-build reduces to a contraction between a categorically compressed expansion vector and categorically compressed STP-DAS one-electron excitation lists. Thus, Eqs. (5) to (7) yield a categorically compressed \sigma-vector, in a manner directly analogous to the compression-preserving behavior of sparse matrix-sparse vector products (SpMSpV). Importantly, this compression preservation is general and independent of the specific storage format used for the categorically compressed representations.

The categorically compressed representation can be implemented in various forms, with the choice of storage format guided primarily by computational efficiency. Since the categorical \sigma-build algorithm often involves reading numerous expansion coefficients with increasing local addresses during contraction, it is natural to adopt a compressed sparse column (CSC) format for storing the categorical expansion coefficients.

Eigenvalue bound analysis

We wish to apply the gap theorem 94,95,136 to bound the error in the computed X2C-CASCI ground state energy:

|\delta E| \leq \frac{\|r\|^2}{\gamma_0}, (8)

where r is the Ritz residual of the computed ground state and the gap \gamma_0 \equiv E_1 - E_0 is the difference between E_1, the (unknown) exact energy of the first excited state, and E_0, the computed Ritz value of the ground state. Because \gamma_0 is unknown, one can estimate its order-of-magnitude using approximate methods or use experimental values to compute a surrogate for the true gap. One can also obtain an exact lower bound on the gap by including the posterior error bound of the first excited state 136–138 in the Davidson calculation 100:

\gamma_0 = E_1 - E_0 \geq E_1 - \|r_1\| - E_0 - \|r_0\|, (9)

where r_1 is the residual associated with the Ritz value E_1.

Using the gap theorem 94,95,136, we can place an exact bound on the error in the computed X2C-CASCI ground-state energy. This requires an estimate of the energy gap between the ground and first excited states of the X2C-CASCI Hamiltonian. To obtain this, we performed an X2C-CISD calculation for the two lowest-lying states and determined a gap of ~0.095 a.u. for HBrTe. Based on this estimate, the gap theorem bounds the error in our X2C-CASCI ground-state energy for HBrTe to within 10 microhartree, which is well below any chemically meaningful threshold.

There is a potential workload imbalance associated with Algorithm 1: because the collection of categorical expansion vectors is distributed among computing nodes, the contraction workload of a given node is proportional to the number of nonzero categorical expansion coefficients it has. We alleviated some of the resulting computational delay by implementing passive one-sided MPI communication of categorical expansion vectors using remote memory access (RMA). This allows idle nodes to contract more categorical expansion vectors with their excitation lists without waiting for the corresponding busy nodes to broadcast them. Ultimately, overcoming this load-balancing issue requires dynamically redistributing categorical expansion vectors according to their sparsity, which changes during the iterative diagonalization.

As illustrated in Algorithm 1, the computational cost, both in memory and runtime, of the compression-compatible STP-DAS \sigma-build procedure increases with the density of the expansion vectors. Therefore, it is essential to maintain maximal compression in these vectors. To achieve this, we replace the traditional Davidson preconditioner with a compression-compatible alternative for generating new trial expansion vectors. This modification alters only the subspace expansion strategy in the Davidson algorithm, while preserving exact treatment of the full determinantal space. As a result, the computed matrix-vector product HC and the corresponding residual norm \|r\| = \|HC - \lambda C\| remain exact, unlike in selected CI and other truncated approaches, where both the Hamiltonian and CI vectors are explicitly approximated.

The compression-compatible STP-DAS \sigma-build algorithm significantly reduces the overall workload associated with the \sigma-build. The reduction in workload can be nonuniform: the contraction workload associated with a determinant \langle JC \mu \nu \kappa \lambda \rangle is proportional to the number of local addresses containing both nonzero categorical expansion vectors elements and Hamiltonian matrix elements. To improve the load-balance, we implemented dynamic SMP thread-level parallelism in the outermost loop (line 1 in Algorithm 1) instead of in the loop over determinants \langle JC \mu \nu \kappa \lambda \rangle (line 4 in Algorithm 1). Under dynamic parallelism, some threads execute many light contractions, while others execute fewer heavy contractions, resulting in a more uniform distribution of contraction workload. Such dynamic parallelism is ineffective in the loop over determinants \langle JC \mu \nu \kappa \lambda \rangle due to the small tensor product nature of the STP-DAS framework.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript reports a lareg configuration interaction calculation, claiming to include over quadrillion determinants. A true FCI calculation of this size is unprecedented. This is surely an impressive calculation, but the choice of system inflates the performance of this algorithm compared to a system of low symmetry and strong correlation. Also, the compression scheme that enables the speed up is not clearly described, therefore some claims are not supported. The current manuscript does not back up the claim of over a quadrillion determinants.\n1. Atoms have high symmetry, which causes many matrix and CI elements that are exactly zero. It appears that symmetry is not used, which stretches how impressive the calculations sound by introducing many determinants of different irreps not coupled to the ground state (many zeros that are easily compressed). There are not near 1 quadrillion coupled determinants. How would this algorithm perform in a system with no symmetry, where all quadrillion determinants may couple?\n2. Key to this work is lossless categorical compression of excitation lists and CI vectors, however, the compression algorithm is not given detailed description. A detailed description is needed to show it is truly lossless, help the reader understand how such efficient compression is possible, and for reproducibility.\n3. Related, Mg has only two valence electrons and is not strongly correlated. CISD is probably fine for this system, thus compression is easier. This is npt the usual case for use of FCI. Application to a strongly correlated system is needed to demonstrate the power of the algorithm for a true FCI application.", "reviewer_2": "The authors have reported the largest numerically exact FCI calculation till date involving one quadrillion (10^15) determinants enabled by an approach called lossless categorical compression within the small tensor product distributed active space (STP-DAS). They have performed the calculations to obtain relativistic FCI ground state of a magnesium atom involving over 10^15 complex valued 2-spinor determinants. The benchmarks shows the numerical exactness of the method with reduced resource utilization. They have eliminated the memory bottlenecks for storing excitation lists which was present in the original STP-DAS work (Ref. 30), while in this work they focused upon reducing the computing resources to store CI vectors by several orders of magnitude.\nThis manuscript shows extremely promising and exciting results, pushing the ultimate limits for state-of-the-art CI calculations. However, despite showing promising results, I would rate it more of a derivative work of the original STP-DAS framework, the only difference being the categorical compression approximation in the construction of sigma-build vectors and during Davidson’s pre-conditioner calculations. I think this is more suited for a more technically focused journal of the relevant field given the original idea got published previously and hence I do not recommend its publication in NatComm. However, I must clearly mention that I do not have any doubts regarding the scientific rigor of this work.\nAdditionally, I have a few suggestions that could enhance the quality of the manuscript:\n1. The manuscript contains relatively elaborate discussions on selected configuration interaction (SCI) methods but not on the main idea/ motive behind STP-DAS framework. The manuscript would greatly benefit from an discussion on STP-DAS.\n2. While STP-DAS has been developed by some of the authors in an earlier work (Ref. 30), this work is solely based on the “categorical compression” approximation of the same. The manuscript does not contain even a single line explaining what the authors actually mean by the term “categorical compression”. Although it is mentioned once in the supplementary file, but I would suggest the main idea should be contained in the main text, and not in the supporting document.\n3. The paper is written in an extremely technical language and lacks a clear explanation behind the approximation presented. It should be written in a more coherent manner for better readability. The authors should put equal stress upon the motive behind the categorical compression approximation and why it works along with the already written sections on ‘how it works’.\nOverall, this is a technically sound and robust work, but given the original idea got published previously, I would recommend that this article is not suited for Nat. Comm., but certainly is publishable in a more technical journal.", "reviewer_3": "The paper describes two very important method improvements for full configuration interaction calculations, “STP-DAS” and “categorically compressed CI vectors”, as well as good load balance needed for good HPC performance.\nHowever, their report of breaking the quadrillion barrier is only formally correct because they do not use symmetry. Their system is an atom, Magnesium, which has a lot of symmetry and therefore a major fraction of the quadrillion+ coefficients will be zero because of symmetry (J, MJ and parity for the relativistic calculation). Furthermore, spin-orbit coupling is not that strong in Mg, leading to many more coefficients being negligible because LS-coupling is almost ok, i.e. they will be excluded by the preconditioner in eqs. 1-2.\nIn order for me to accept a claim of breaking the quadrillion barrier, the CI must have quadrillion coefficients which we do not know is zero by symmetry. (If it had been a non-linear molecule I would have said OK, but for me the Mg calculation is blowing up the number of coefficients artificially to be able to say the quadrillion barrier is broken.)\nIt is excellent work and the paper should definitely be published somewhere, but not with this title. Or a real quadrillion calculation should be included in the paper.\nHans Jørgen Aa. Jensen", "reviewer_4": "The manuscript presents a significant extension of the previously proposed small-tensor-product distributed active space (STP-DAS) framework through the implementation of lossless categorical compression. This technical innovation allows configuration interaction (CI) calculations involving up to a quadrillion determinants — reportedly the largest number of determinants ever used in such a calculation. This represents a notable computational advancement and shows the authors' ability to push the boundaries of tensor factorization and distributed-memory computation in quantum chemistry. The methodology is sound, the technical details are well documented, and the manuscript is clearly written and should be published. However, before acceptance, I recommend that the following points be addressed.\nThe authors provide access to chronus github, but not explicitly to the code used in this work, which limits reproducibility and the ability of the community to build upon these results. At the very least, the authors should include a clear statement in the manuscript on how this code can be used and some examples of input. If the code cannot be made public, justification should be provided.\nWhile the abstract and main text mention benchmark calculations on N₂ and a carbon nanotube model system, there are no clear results shown for these systems either in the main text or the supplementary material. The authors should clarify exactly what was computed for these systems, and explicitly include representative data — such as energies or basis set convergence — in the paper or supporting information.\nThe calculation on the Mg atom (or Mg2+??? Not clear) appears to serve only as a demonstration that the method can handle a quadrillion determinants. However, the Mg atom is a closed-shell system with modest electron correlation, and it is unclear what is gained chemically or physically from such an extensive calculation on this target. The authors should explain why this example was chosen and discuss future applications that genuinely require this level of computational complexity. For example, do the authors envision studying strongly correlated bond dissociations, excited states with high multi-configurational character, or transition metal complexes? Without such a forward-looking discussion, the practical relevance of the current study is diminished.\nWhile the paper introduces technical improvements, the distinction between this work and the prior STP-DAS study by the same authors is not sufficiently emphasized. The authors should clearly articulate how these advancements go beyond their earlier framework. More importantly, they should elaborate on the broader implications for quantum chemistry. How does this work alter the practical frontier of CI-based electronic structure theory? Can the compression techniques introduced here accelerate real chemical discovery?\nAs it stands, the manuscript reads as a detailed and well-executed technical report. However, the novelty and potential impact of the work should be more clearly argued in order to appeal to a broader readership and to justify publication in a high-profile journal.\nMinor Comments:\n1. When describing the GASSCF method, the following reference should be included to provide appropriate historical context: Ma, D.; Li Manni, G.; Gagliardi, L. “The Generalized Active Space Concept in Multiconfigurational Self-Consistent Field Methods.” J. Chem. Phys. 2011, 135, 044128. https://doi.org/10.1063/1.3611401\n2. In the discussion of post-SCF methods, MC-PDFT should be cited alongside CASPT2 and NEVPT2, as it represents an increasingly adopted multireference method with favorable scaling and performance. A suitable reference is: Li Manni, G.; Carlson, R. K.; Luo, S.; Ma, D.; Olsen, J.; Truhlar, D. G.; Gagliardi, L. “Multiconfiguration Pair-Density Functional Theory.” J. Chem. Theory Comput. 2014, 10, 3669–3680. https://doi.org/10.1021/ct500483t\nIn summary, this is an impressive technical achievement in large-scale CI computations. However, the authors should better justify the scientific relevance of their results and clarify how this work significantly advances beyond previous efforts. With these revisions, the paper will be much stronger and suitable for publication."}, "author_response": "Response to Reviewer 1:\nReviewer 1: This manuscript reports a large configuration interaction calculation, claiming\nto include over quadrillion determinants. A true FCI calculation of this size is unprecedented.\nThis is surely an impressive calculation, but the choice of system inflates the performance of this\nalgorithm compared to a system of low symmetry and strong correlation. Also, the compression\nscheme that enables the speed up is not clearly described, therefore some claims are not supported.\nThe current manuscript does not back up the claim of over a quadrillion determinants.\nAuthor Response: We are grateful to the Reviewer for recognizing this as an impressive cal-\nculation of unprecedented scale. In response to the comments, we have extensively revised the\nmanuscript to include details of the algorithmic innovations and three additional benchmark sys-\ntems, each highlighting distinct challenges: low symmetry with strong spin–orbit coupling, pro-\nnounced static correlation, and significant dynamic correlation.\nWe also wish to emphasize that, while most of the community may not have access to the hundreds\nof compute nodes required for such extreme-scale runs, this work demonstrates the practical\nfeasibility of trillion-determinant calculations on just a few nodes and even on a laptop.\nWe hope that the revised discussion and added benchmarks satisfactorily address the Reviewer’s\nconcerns.\nReviewer 1: 1. Atoms have high symmetry, which causes many matrix and CI elements that are\nexactly zero. It appears that symmetry is not used, which stretches how impressive the calculations\nsound by introducing many determinants of different irreps not coupled to the ground state (many\nzeros that are easily compressed). There are not near 1 quadrillion coupled determinants. How\nwould this algorithm perform in a system with no symmetry, where all quadrillion determinants\nmay couple?\nAuthor Response: The Reviewer is correct that the magnesium atom is a highly symmetric\nsystem, which intrinsically reduces the dimensionality of the problem. The original choice of Mg\nwas not due to any limitation of the algorithm, but rather to our restricted access to supercom-\nputing resources at the time. After receiving your comments, we submitted a new proposal for\nadditional resources and were fortunate to be awarded extra computing time.\nWith this allocation, we have now computed the ground state of HBrTe (100 two-spinor orbitals,\n88 electrons, 1.05 × 1015 two-spinor determinants), a nonlinear molecule with strong spin-orbit\ncoupling. The calculation required only 34.5 hours over 9 iterations to fully converge the CI\nground state. The changes are highlighted below and included in Table 1 on Page 5. All memory\nanalyses in the manuscript have been updated to incorporate this new test case.\nPage 5: Before presenting methodological details and performance benchmarks, we highlight the\nlargest CASCI calculation to date for the ground state of HBrTe, enabled by the compression-\ncompatible STP-DAS framework to be introduced herein. HBrTe is a substituted form of a\nhydrogen chalcogenide where one of the hydrogens was substituted with a bromine atom to de-\ncrease the symmetry to the molecule. We performed a relativistic exact-two-component CASCI\n(X2C-CASCI) calculation (100 2-spinor orbitals, 88 electrons, 1.05×1015 2-spinor determinants)\nof the ground state of HBrTe using the compression-compatible STP-DAS framework. We also\nperformed an analogous calculation (with the same number of determinants) for the ground state\nof a magnesium atom (see Section S5 of the Supporting Information).\nverged in 9 iterations to microhartree precision (< 10−6 au) with a total runtime of 34.5 hours. In\neach iteration, an additional CI expansion vector was introduced to accelerate convergence, while\nthe compression algorithm dynamically adapted to the expanding vector space, leading to a grad-\nual increase in computational cost. A total of 9 expansion vectors were involved in the σ-build.\nDespite the enormous configuration space of over one quadrillion (1.05 × 1015) complex-valued\n2-spinor determinants, the average σ-build time per vector remained just 3.8 hours.\nReviewer 1: 2. Key to this work is lossless categorical compression of excitation lists and CI\nvectors, however, the compression algorithm is not given detailed description. A detailed descrip-\ntion is needed to show it is truly lossless, help the reader understand how such efficient compression\nis possible, and for reproducibility.\nAuthor Response: We appreciate the Reviewer for raising this important point. In response,\nwe have added the following statements to further illustrate the concept of ‘categorical compres-\nsion.’ We believe this clarification highlights how the compression scheme works. The numerical\nexactness of the algorithm is ultimately evident in the results of the CI calculations (the exact\nerror bounds on the HBrTe, Rb , and Xe calculations, as well as the calculations in Section S3\n4 2\nof the Supporting Information).\nPage 7ThecentralconceptoftheSTP-DASframeworkisthesystematicpartitioningofthefullCI\norbitalspaceintoacollectionofdistributedactivespaces. Withineachactivespace, configurations\nare further classified into categorical subspaces, rigorously defined by distinct electron occupation\npatterns, as illustrated in Figure 2(A).\nPage 8 With the STP-DAS framework, the memory bottleneck associated with storing CI expan-\nsion vectors is effectively eliminated by applying numerically exact categorical compression. The\nSTP-DAS CI expansion vectors take the form\n(cid:77)\nC = CB, (1)\nB\nwhere B is a category, defined by a unique electron occupation pattern within the distributed\nactive spaces.\nThe compression scheme, categorical compression (see Figure 2(D)), stores the CB’s as a com-\npressed sparse column (CSC) vector. In this format, each categorical expansion vector CB is\nrepresented by two equally sized arrays: one, VCB ≡ {C : C ̸= 0}, stores the values of the\nCSC KB KB\nnonzero coefficients, while the other, ACB ≡ {KB : C ̸= 0}, stores their corresponding local\nCSC KB\naddresses. The tensor-loop structure in the STP-DAS σ-build algorithm is reformulated in terms\nof categorically compressed local addresses together with their corresponding global phase factors\n(see Figure 2(E)).\nIn contrast to element-wise compression, categorical compression can eliminate all configurations\nwithin a category, i.e., skip an entire category at once. This approach is better able to preserve the\nvectorized structure compared to element-wise compression. The major difficulty of any sparse\nmatrix-vector product algorithm is the lack of a priori knowledge of the location of the nonzero\nelements. This difficulty leads to a major bottleneck rooted in the nonuniform accesses to memory.\nThe categorical compression localizes nonuniform memory accesses within a category, which is\nalways orders-of-magnitudes smaller than the size of the full CI space. This localized memory\naccess pattern is the central strength of the categorical compression scheme.\nσ-build algorithm to bypass entire categories of determinants while preserving both the vectorized\nstructure and the local addressing scheme of STP-DAS. Because the compression is fully lossless,\nthe omitted determinants have no impact on the resulting Ritz eigenvalue–eigenvector pair.\nReviewer 1: 3. Related, Mg has only two valence electrons and is not strongly correlated.\nCISD is probably fine for this system, thus compression is easier. This is not the usual case for\nuse of FCI. Application to a strongly correlated system is needed to demonstrate the power of the\nalgorithm for a true FCI application.\nAuthor Response: To demonstrate the performance of our algorithm on highly correlated sys-\ntems, we ran an X2C-CASCI calculation on Rb , which is a relativistic counterpart of the H\n4 4\nsystem. This calculation, which contains nearly 89 trillion complex-valued 2-spinor determinants,\nshowcases the ability of categorical compression within the STP-DAS framework to tackle highly\nstatically correlated systems.\nIn addition, we performed an X2C-CASCI calculation on a Xenon dimer, which is known to have\nsignificant dynamic correlation. This calculation contains 1.4 trillion determinants, and\ndemonstrates how strongly correlated systems can be treated within our framework with a high\ndegree of precision. Note that because the ground state of this system couples many determinants,\nthe diagonalization takes a long time. Therefore, because we are restricted to 48 hour-long runs\non NERSC’s Perlmutter supercomputer, we chose a smaller active space for this system. Without\nthis limitation, we would be able to enlarge this active space arbitrarily, as showcased by the\nHBrTe system. The changes from these two systems are highlighted below and included in Tables\n2 and 3 on Page 7.\nPage 6: Additionally, we performed X2C-CASCI calculations for the ground states of two highly\ncorrelated systems, showcasing the applicability of the compression-compatible STP-DAS frame-\nworktohighlycorrelatedsystems. Thesesystemswerechosenfromoppositeendsofthecorrelation\nspectrum from strongly statically correlated to strongly dynamically correlated. Square Rb ,\n4\nrelativistic analogue of H , displays strong static correlation and Xe is a dynamically correlated\n4 2\nnoble gas dimer.\nTables 2 and 3 summarize the results of the Rb and Xe calculations. The Rb calculation\n4 2 4\n(50 2-spinor orbitals, 28 electrons, 8.9×1013 2-spinor determinants) ran on the National Energy\nResearch Scientific Computing Center’s Perlmutter high-performance supercomputer with a total\nof 100 nodes (AMD EPYC 7763 Milan, 12,800 compute cores, 512 GB of RAM per node, 200\nGB/s NIC, 2 MPI ranks per node, and 64 SMP threads per MPI process) and took 6 iterations\nand 11.8 hours to converge. The Xe calculation (60 2-spinor orbitals, 12 electrons, 1.4 × 1012\n2\n2-spinor determinants) ran on the same platform with 256 nodes and took 7 iterations and 36.1\nhours to converge.\nTheground-stateenergiesoftheRb andXe moleculesfromthesecalculationsare−11916.1527249\n4 2\nau and −14889.6466956 au, accordingly. The gap theorem guarantees that these X2C-CASCI re-\nsults lie within 0.53×10−6 au and 17.56×10−6 au of the true ground state energies within the\ncorresponding active spaces and basis sets (see Section S1 of the Supporting Information).\nReviewer 2: The authors have reported the largest numerically exact FCI calculation till date\ninvolving one quadrillion (1015) determinants enabled by an approach called lossless categorical\ncompression within the small tensor product distributed active space (STP-DAS). They have per-\nformed the calculations to obtain relativistic FCI ground state of a magnesium atom involving\nover 1015 complex valued 2-spinor determinants. The benchmarks shows the numerical exactness\nof the method with reduced resource utilization. They have eliminated the memory bottlenecks for\nstoring excitation lists which was present in the original STP-DAS work (Ref. 30), while in this\nwork they focused upon reducing the computing resources to store CI vectors by several orders of\nmagnitude. This manuscript shows extremely promising and exciting results, pushing\nthe ultimate limits for state-of-the-art CI calculations.\nAuthor Response: We are delighted that the Reviewer recognized this work as a major break-\nthrough in electronic structure methodology. In response to the comments, we have extensively\nrevisedthemanuscripttoincludedetailsofthealgorithmicinnovationsandthreeadditionalbench-\nmarksystems,eachhighlightingdistinctchallenges: lowsymmetrywithstrongspin–orbitcoupling,\npronounced static correlation, and significant dynamic correlation.\nWe also wish to emphasize that, while most of the community may not have access to the hundreds\nof compute nodes required for such extreme-scale runs, this work demonstrates the practical\nfeasibility of trillion-determinant calculations on just a few nodes and even on a laptop.\nWe hope that the expanded discussion and more challenging examples will demonstrate to the\nReviewer that this work is well suited for Nature Communications themed collection on\n“Methodological Developments in Electronic Structure Theory”.\nReviewer 2: However, despite showing promising results, I would rate it more of a derivative\nwork of the original STP-DAS framework, the only difference being the categorical compression\napproximation in the construction of sigma-build vectors and during Davidson’s pre-conditioner\ncalculations. I think this is more suited for a more technically focused journal of the relevant field\ngiven the original idea got published previously and hence I do not recommend its publication in\nNatComm. However, I must clearly mention that I do not have any doubts regarding the scientific\nrigor of this work.\nAuthor Response: The STP-DAS algorithm was previously published, demonstrating that a\nbillion determinants can be computed on a laptop and a trillion on a supercomputer. In this\nwork, we extend both regimes by three orders of magnitude, achieving a trillion determinants on\na few nodes and a quadrillion on a supercomputer, enabled by the innovation of numerically exact\ncategorical compression.\nThe present work is not merely a derivative of our previous publication. The compression scheme\nintroduced here is innovative, as it is formulated in categorical space rather than relying solely on\nelement-wise compression. In practice, this allows the algorithm to skip entire categories as well as\nelementswithinthem,yieldingamorevectorizedstructurethanelement-wisecompressionthatcan\nfully exploit the STP-DAS framework. Implementing categorical compression, however, required a\ncomplete restructuring of the fundamental STP algorithm to accommodate this advance. Indeed,\nit represents a substantive methodological step forward, achieved through the contributions of a\nnew generation of researchers.\nWe hope that with this extensive revision, the Reviewer will recognize that this work is highly\n4\nappropriate for the Nature Communications themed collection on ‘Methodological Developments\nin Electronic Structure Theory.’\nReviewer 2: 1. The manuscript contains relatively elaborate discussions on selected configu-\nration interaction (SCI) methods but not on the main idea/ motive behind STP-DAS framework.\nThe manuscript would greatly benefit from an discussion on STP-DAS.\nAuthor Response: We thank the Reviewer for the suggestion. In this revision, we have revised\nthe following paragraph in the introduction:\nPage 4 A recently introduced CI matrix–vector product algorithm leverages the exact factoriza-\ntionoftheactivespaceintosmalltensorproductsofdistributedactivespaces—anapproachknown\nas the small-tensor product distributed active space (STP-DAS) framework, illustrated in Figure\n2(A). The STP-DAS algorithm reformulates the large CI matrix–vector product into a sequence of\nsmall tensor products, each embedded within a distributed active space, computed on-the-fly us-\ning string-based methods. This advance exploits the mathematical condition governing the phase\nrelationship between the global address and the local DAS address of any CI matrix element,\nenabling the use of only small local determinant address strings in the CI matrix–vector build and\novercoming the memory bottleneck associated with storing the full excitation list. This formu-\nlation enables extensive reuse of Hamiltonian excitation lists, leading to a dramatic reduction in\nmemory demands. For a CI problem involving one quadrillion (1015) determinants, the STP-DAS\nframework reduces the excitation list memory requirement from 12 exabytes (EB) to 25 gigabytes\n(GB), an 8-orders-of-magnitude reduction! In addition, by evenly distributing the computation\nof small tensor products, the STP-DAS algorithm achieves excellent load balance with minimal\nnode-to-node communication overhead, ensuring strong scalability across both single-node and\nlarge-scale parallel architectures.\nReviewer 2: 2. While STP-DAS has been developed by some of the authors in an earlier work\n(Ref. 30), this work is solely based on the “categorical compression” approximation of the same.\nThe manuscript does not contain even a single line explaining what the authors actually mean by\nthe term ”categorical compression”. Although it is mentioned once in the supplementary file, but\nI would suggest the main idea should be contained in the main text, and not in the supporting\ndocument.\nAuthor Response: We appreciate the Reviewer for raising this important question. In response,\nwe have added the following statements to further illustrate the concept of ‘categorical compres-\nsion.’ We believe this clarification also highlights that the present work is not merely a simple\nderivative of our previous study.\nPage 7ThecentralconceptoftheSTP-DASframeworkisthesystematicpartitioningofthefullCI\norbitalspaceintoacollectionofdistributedactivespaces. Withineachactivespace, configurations\nare further classified into categorical subspaces, rigorously defined by distinct electron occupation\npatterns, as illustrated in Figure 2(A).\nPage 8 With the STP-DAS framework, the memory bottleneck associated with storing CI expan-\nsion vectors is effectively eliminated by applying numerically exact categorical compression. The\nSTP-DAS CI expansion vectors take the form\n(cid:77)\nC = CB, (2)\nB\n5\nwhere B is a category, defined by a unique electron occupation pattern within the distributed\nactive spaces.\nThe compression scheme, categorical compression (see Figure 2(D)), stores the CB’s as a com-\npressed sparse column (CSC) vector. In this format, each categorical expansion vector CB is\nrepresented by two equally sized arrays: one, VCB ≡ {C : C ̸= 0}, stores the values of the\nCSC KB KB\nnonzero coefficients, while the other, ACB ≡ {KB : C ̸= 0}, stores their corresponding local\nCSC KB\naddresses. The tensor-loop structure in the STP-DAS σ-build algorithm is reformulated in terms\nof categorically compressed local addresses together with their corresponding global phase factors\n(see Figure 2(E)).\nIn contrast to element-wise compression, categorical compression can eliminate all configurations\nwithin a category, i.e., skip an entire category at once. This approach is better able to preserve the\nvectorized structure compared to element-wise compression. The major difficulty of any sparse\nmatrix-vector product algorithm is the lack of a priori knowledge of the location of the nonzero\nelements. This difficulty leads to a major bottleneck rooted in the nonuniform accesses to memory.\nThe categorical compression localizes nonuniform memory accesses within a category, which is\nalways orders-of-magnitudes smaller than the size of the full CI space. This localized memory\naccess pattern is the central strength of the categorical compression scheme.\nIn summary, the numerically exact categorical compression introduced here allows the STP-DAS\nσ-build algorithm to bypass entire categories of determinants while preserving both the vectorized\nstructure and the local addressing scheme of STP-DAS. Because the compression is fully lossless,\nthe omitted determinants have no impact on the resulting Ritz eigenvalue–eigenvector pair.\nReviewer 2: 3. The paper is written in an extremely technical language and lacks a clear\nexplanation behind the approximation presented. It should be written in a more coherent manner\nfor better readability. The authors should put equal stress upon the motive behind the categorical\ncompression approximation and why it works along with the already written sections on ‘how it\nworks’.\nAuthor Response: We greatly appreciate this suggestion to make the work more accessible to a\nbroader audience. In response, we have extensively revised the manuscript to emphasize not only\nhow it works but also why it works. We hope this revision meets the Reviewer’s high standards\nof scientific clarity and accessibility.\nWe wish to emphasize that this paper focuses on methodological development, in alignment with\nthe Nature Communications themed collection on ‘Methodological Developments in Electronic\nStructure Theory.’ While technical in nature, the manuscript uses benchmark examples to moti-\nvate the scientific community toward developing a new generation of CI codes—an endeavor that\nis by no means a task for non-experts.\nReviewer 3: The paper describes two very important method improvements for full configuration\ninteraction calculations, “STP-DAS” and “categorically compressed CI vectors”, as well as good\nload balance needed for good HPC performance.\nHowever, their report of breaking the quadrillion barrier is only formally correct because they do not\nuse symmetry. Their system is an atom, Magnesium, which has a lot of symmetry and therefore\na major fraction of the quadrillion+ coefficients will be zero because of symmetry (J, MJ and\nparity for the relativistic calculation). Furthermore, spin-orbit coupling is not that strong in Mg,\nleading to many more coefficients being negligible because LS-coupling is almost ok, i.e. they will\nbe excluded by the preconditioner in eqs. 1-2.\nIn order for me to accept a claim of breaking the quadrillion barrier, the CI must have quadrillion\ncoefficients which we do not know is zero by symmetry. (If it had been a non-linear molecule\nI would have said OK, but for me the Mg calculation is blowing up the number of coefficients\nartificially to be able to say the quadrillion barrier is broken.)\nIt is excellent work and the paper should definitely be published somewhere, but not with this title.\nOr a real quadrillion calculation should be included in the paper.\nAuthor Response: The original choice of Mg was not due to any limitation of the algorithm,\nbut rather to our restricted access to supercomputing resources at the time. After receiving your\ncomments, we submitted a new proposal for additional resources and were fortunate to be awarded\nextra computing time.\nIn this revision, we included an X2C-CASCI calculation for the ground state of HBrTe containing\n1.05×1015 2-spinor determinants in addition to the original Mg atom calculation. This system is\na nonlinear molecule and has minimal symmetry and strong spin-orbit coupling. The successful\nconvergence of this calculation showcases the performance of our algorithm on such systems.\nTo demonstrate the performance of our algorithm on highly correlated systems, we carried out\nX2C-CASCI calculations on Rb (89 trillion complex-valued 2-spinor determinants), the rela-\n4\ntivistic analogue of the H system, and on a xenon dimer (1.4 trillion complex-valued 2-spinor\n4\ndeterminants), a system well known for its significant dynamic correlation.\nThe changes are highlighted below and included in Table 1 on Page 5 for HBrTe. All memory\nanalyses in the manuscript have been updated to incorporate this new test case.\nPage 5: Before presenting methodological details and performance benchmarks, we highlight the\nlargest CASCI calculation to date for the ground state of HBrTe, enabled by the compression-\ncompatible STP-DAS framework to be introduced herein. HBrTe is a substituted form of a\nhydrogen chalcogenide where one of the hydrogens was substituted with a Bromine to decrease\nthe symmetry to the molecule. We performed a relativistic exact-two-component CASCI (X2C-\nCASCI) calculation (100 2-spinor orbitals, 88 electrons, 1.05 × 1015 2-spinor determinants) of\nthe ground state of HBrTe using the compression-compatible STP-DAS framework. We also\nperformed an analogous calculation (with the same number of determinants) for the ground state\nof a magnesium atom (see Section S5 of the Supporting Information).\nPage 5: Table 1 summarizes the results of the HBrTe calculation. The ground-state energy\nconverged in 9 iterations to microhartree precision (< 10−6 au) with a total runtime of 34.5 hours.\nIn each iteration, an additional CI expansion vector was introduced to accelerate convergence,\nwhile the compression algorithm dynamically adapted to the expanding vector space, leading to\n7\na gradual increase in computational cost. A total of 9 expansion vectors were involved in the\nσ-build. Despite the enormous configuration space of over one quadrillion (1.05×1015) complex-\nvalued 2-spinor determinants, the average σ-build time per vector remained just 3.8 hours.\nReviewer 4: The manuscript presents a significant extension of the previously proposed small-\ntensor-product distributed active space (STP-DAS) framework through the implementation of loss-\nless categorical compression. This technical innovation allows configuration interaction (CI) calcu-\nlations involving up to a quadrillion determinants — reportedly the largest number of determinants\never used in such a calculation. This represents a notable computational advancement and shows the\nauthors’ ability to push the boundaries of tensor factorization and distributed-memory compu-\ntation in quantum chemistry. The methodology is sound, the technical details are well documented,\nand the manuscript is clearly written and should be published. However, before acceptance, I rec-\nommend that the following points be addressed.\nAs it stands, the manuscript reads as a detailed and well-executed technical report. However, the\nnovelty and potential impact of the work should be more clearly argued in order to appeal to a\nbroader readership and to justify publication in a high-profile journal.\nIn summary, this is an impressive technical achievement in large-scale CI computations. However,\nthe authors should better justify the scientific relevance of their results and clarify how this work\nsignificantly advances beyond previous efforts. With these revisions, the paper will be much stronger\nand suitable for publication.\nAuthor Response: We are grateful for the Reviewer’s positive comments. In this revision,\nwe have extensively revised the manuscript and included additional examples: a quadrillion-\ndeterminant low-symmetry strong spin–orbit case and two highly correlated systems. In addi-\ntion, we have included the authors’ perspectives on how this work may influence the frontier of\nelectronic structure calculations and contribute to the development of quantum computing. We\nhope these additions further highlight the robustness, versatility, and broad applicability of the\nSTP-DAS framework.\nReviewer 4: The authors provide access to chronus github, but not explicitly to the code used in\nthis work, which limits reproducibility and the ability of the community to build upon these results.\nAt the very least, the authors should include a clear statement in the manuscript on how this code\ncan be used and some examples of input. If the code cannot be made public, justification should be\nprovided.\nAuthor Response: The code is currently accessible within the developers’ community. In con-\ntrast to other open-source codes where frequent releases are managed by the community, we only\nmake annual releases of Chronus. This choice reflects the need for extensive unit testing and, for\nthe STP-DAS module in particular, rigorous stress testing on multiple supercomputing facilities.\nThat said, we plan to release the STP-DAS module in December 2025.\nWith that said, we created a containerized version of the implementation we used to produce the\ndata in this manuscript, which is freely available at Docker Hub in the following organization\nand repository: uwligroup/chronusq:preview. We have added the following statement to the\nmanuscript:\nPage 10: The STP-DAS implementation currently accessible within the developers’ commu-\nnity. It will be included in the next official Chronus public release, following extensive testing\nacross multiple platforms to ensure robustness, scalability, and reproducibility. In the meantime,\nwe created an MPI-compatible, portable, and containerized version of Chronus Quantum con-\ntaining the compression-compatible STP-DAS implementation. This version has minimal over-\n9\nhead and is freely available at Docker Hub in the following organization and repository: uwli-\ngroup/chronusq:preview.\nReviewer 4: While the abstract and main text mention benchmark calculations on N and a\n2\ncarbon nanotube model system, there are no clear results shown for these systems either in the main\ntext or the supplementary material. The authors should clarify exactly what was computed for these\nsystems, and explicitly include representative data — such as energies or basis set convergence —\nin the paper or supporting information.\nAuthor Response: The results for N and a carbon nanotube model are provided in Table S2\n2\nof the Supporting Information to illustrate the STP-DAS convergence profiles. We have added\nadditional clarification statements in the Supporting Information.\nSupporting Information: The convergence of the change in energy (|δE|) and the norm of the\nresidual(∥r∥) with varying truncation level (ε) is explored.\nIt is evident in Table S2, where all three systems, regardless of the level of correlation, attain\neigenvalues within 10−7 E of their exact X2C-CASCI counterparts when the residual norms ∥r∥\nh\nare small. For loose thresholds, new trial expansion vectors quickly become linearly dependent\non the existing subspace expansion vectors, indicating that the span of the modified expansion\nvectors has grown to its limit before convergence. The resulting energies at these thresholds differ\nsignificantly from their exact counterparts, as expected. Furthermore, the residual norm remains\nfairly high, indicating that the Ritz vector is still far from the true eigenvector.\nReviewer 4: The calculation on the Mg atom (or Mg2+??? Not clear) appears to serve only\nas a demonstration that the method can handle a quadrillion determinants. However, the Mg\natom is a closed-shell system with modest electron correlation, and it is unclear what is gained\nchemically or physically from such an extensive calculation on this target. The authors should\nexplain why this example was chosen and discuss future applications that genuinely require this level\nof computational complexity. For example, do the authors envision studying strongly correlated\nbond dissociations, excited states with high multi-configurational character, or transition metal\ncomplexes? Without such a forward-looking discussion, the practical relevance of the current study\nis diminished.\nAuthor Response: In light of the Reviewer’s remarks, we have extended the benchmark calcu-\nlations to include more challenging systems, such as low-symmetry cases with strong spin–orbit\ncoupling, strongly correlated systems, and systems with pronounced dynamic correlation.\n• We have converged the ground state of a nonlinear HBrTe system containing 1.05×1015 two-\nspinor determinants, thereby demonstrating another quadrillion-determinant calculation in\na low-symmetry system with strong spin-orbit coupling.\n• To demonstrate the performance of the STP-DAS algorithm on highly correlated systems,\nwe ran an X2C-CASCI calculation on Rb (89 trillion 2-spinor determinants), which is a\n4\nrelativistic counterpart of the H system.\n4\n• We performed an X2C-CASCI calculation on a Xenon dimer, which is known to have signif-\nicant dynamic correlation. This calculation contains 1.4 trillion determinants, and demon-\nstrates how strongly correlated systems can be treated within our framework with a high\ndegree of precision.\nhave added the following statements in the Conclusion section.\nPage 10: The integration of categorical compression with STP-DAS marks a paradigm shift in\ntackling large-scale CI problems. As quantum chemistry continues to push the limits of system\ncomplexity, the ability to carry out quadrillion-determinant calculations within tractable resource\nboundsestablishesapowerfulfoundationforstudyinghighlycorrelated, multireference, relativistic\nsystems. While access to hundreds of compute nodes for quadrillion-determinant calculations may\nremain out of reach for most of the community, this work demonstrates the practical feasibility of\ntrillion-determinant calculations on a small cluster.\nFor transition-metal, rare-earth, and heavy-element complexes, such large-scale CI calculations\nenable predictive simulations of electronic structure properties (bond order, covalency, polariza-\ntion, etc.), spectroscopic observables (UV/Vis, X-ray, etc.), and reaction pathways, with the full\norbital space consisting of both metal and ligand orbitals, treated on an equal footing.\nThe ability to simulate a full CI space of 100 orbitals on a classical computer not only challenges\ncurrent notions of quantum supremacy, but also establishes a robust platform for developing and\nbenchmarking new quantum algorithms aimed at achieving chemical accuracy.\nReviewer 4: While the paper introduces technical improvements, the distinction between this\nwork and the prior STP-DAS study by the same authors is not sufficiently emphasized. The\nauthors should clearly articulate how these advancements go beyond their earlier framework. More\nimportantly, they should elaborate on the broader implications for quantum chemistry. How does\nthis work alter the practical frontier of CI-based electronic structure theory? Can the compression\ntechniques introduced here accelerate real chemical discovery?\nAuthor Response: The main innovation of this work is the categorical compression technique,\nwhich advances the STP-DAS framework to the quadrillion-determinant regime on supercomput-\ners and the trillion-determinant regime on just a few nodes. To emphasize this point, we have\nadded content addressing not only how categorical compression works, but also why it is criti-\ncal. Furthermore, we have expanded the Conclusion to include the authors’ perspective on the\npotential impact of this work on advancing the frontier of CI-based electronic structure theory.\nPage 8 With the STP-DAS framework, the memory bottleneck associated with storing CI expan-\nsion vectors is effectively eliminated by applying numerically exact categorical compression. The\nSTP-DAS CI expansion vectors take the form\n(cid:77)\nC = CB, (3)\nB\nwhere B is a category, defined by a unique electron occupation pattern within the distributed\nactive spaces.\nThe compression scheme, categorical compression (see Figure 2(D)), stores the CB’s as a com-\npressed sparse column (CSC) vector. In this format, each categorical expansion vector CB is\nrepresented by two equally sized arrays: one, VCB ≡ {C : C ̸= 0}, stores the values of the\nCSC KB KB\nnonzero coefficients, while the other, ACB ≡ {KB : C ̸= 0}, stores their corresponding local\nCSC KB\naddresses. The tensor-loop structure in the STP-DAS σ-build algorithm is reformulated in terms\nof categorically compressed local addresses together with their corresponding global phase factors\n(see Figure 2(E)).\nwithin a category, i.e., skip an entire category at once. This approach is better able to preserve the\nvectorized structure compared to element-wise compression. The major difficulty of any sparse\nmatrix-vector product algorithm is the lack of a priori knowledge of the location of the nonzero\nelements. This difficulty leads to a major bottleneck rooted in the nonuniform accesses to memory.\nThe categorical compression localizes nonuniform memory accesses within a category, which is\nalways orders-of-magnitudes smaller than the size of the full CI space. This localized memory\naccess pattern is the central strength of the categorical compression scheme.\nIn summary, the numerically exact categorical compression introduced here allows the STP-DAS\nσ-build algorithm to bypass entire categories of determinants while preserving both the vectorized\nstructure and the local addressing scheme of STP-DAS. Because the compression is fully lossless,\nthe omitted determinants have no impact on the resulting Ritz eigenvalue–eigenvector pair.\nReviewer 4: Minor Comments:\n1. When describing the GASSCF method, the following reference should be included to provide\nappropriate historical context: Ma, D.; Li Manni, G.; Gagliardi, L. “The Generalized Active\nSpace Concept in Multiconfigurational Self-Consistent Field Methods.” J. Chem. Phys. 2011,\n135, 044128. https://doi.org/10.1063/1.3611401\n2. In the discussion of post-SCF methods, MC-PDFT should be cited alongside CASPT2 and\nNEVPT2, as it represents an increasingly adopted multireference method with favorable scaling\nand performance. A suitable reference is: Li Manni, G.; Carlson, R. K.; Luo, S.; Ma, D.; Olsen,\nJ.; Truhlar, D. G.; Gagliardi, L. “Multiconfiguration Pair-Density Functional Theory.” J. Chem.\nTheory Comput. 2014, 10, 3669–3680. https://doi.org/10.1021/ct500483t\nAuthor Response: We would like to thank the Reviewer for providing us with the appropriate\nreferences for the GASSCF and the MC-PDFT methods. In response, we cited these references\nat the appropriate places.\nPage 4: Applying more computationally demanding methods such as multiconfigurational self-\nconsistentfield(MCSCF)6−8,12,13,16,41−50,multireferenceconfigurationinteraction(MRCI)6,19,42−48,\nandmany-bodyperturbationtheory(MRPT2, CASPT2, NEVPT2, MC-PDFT)13,14,22,49−55 istyp-\nically required to achieve qualitative and quantitative agreement with experiment.\n12"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Most of my comments are solved, and this is very nice work. However, even for HBrTe (Cs symmetry), there are less than one quadrillion determinants. The claim of 1 quadrillion determinants should be removed.", "reviewer_2": "I am more or less satisfied with the responses from the authors. They have addressed all my comments and included the necessary sections in the revised manuscript. Particularly what is impressive this time is that they have performed calculations on a rather non-trivial system HBrTe in this revised manuscript which shows the true power of the categorical-compression architecture in STP-DAS framework. I recommend the publication of the work in Nature Communications themed collection on \"Methodological Developments in Electronic Structure Theory\".", "reviewer_3": "In the revised manuscript and supplementary material all my comments have been addressed very satisfactorily, now a “real quadrillion” calculation is included. The paper is also well-written, I only found one small error: in the “Author Contribution Statements” is written: “A. Shayit and A. Shayit acquired computing resources.”. I recommend that the manuscript is accepted for publication now.\nHans Jørgen Aa. Jensen", "reviewer_4": "The authors adequately addressed the comments of all reviewers and I am in favor of publication."}, "author_response": "Reviewer 1: Most of my comments are solved, and this is very nice work. However, even for HBrTe (Cs symmetry), there are less than one quadrillion determinants. The claim of 1 quadrillion determinants should be removed.\nAuthorResponse: We note that while the HBrTe calculation reported here is formally conducted\nin a quadrillion-dimensional Hilbert space. The excitation list constructed without any knowledge\nof the symmetry of the desired state. In this sense, all determinants are inherently included in our\ncalculation. This is consistent with the counting reported in previous work by other groups.\nMoreover, the determinant space of the computed state is equivalent to an non-relativistic quadrillion-\ndimensional Hilbert space. Because we conducted a relativistic two-component calculation, the\ndouble that of an analogous non-relativistic calculation. Additionally, the complex nature of the\nof a non-relativistic calculation consisting of over twice as many determinants.\nLastly, we note that we selected HBrTe because of the concrete suggestion made by reviewer\n3 during the initial review (NCOMMS-25-40041) of our original submission (If it had been a\nnon-linear molecule I would have said OK). We believe that for the reason mentioned above,\nreviewer 3 now acknowledge our calculation is a “real quadrillion calculation”.\nPage 5: The excitation list is generated without assuming any symmetry of the target state.\nConsequently, the calculation is formally performed in a quadrillion-determinant space, with all\nis twice that of an analogous non-relativistic calculation. Moreover, the complex arithmetic make\nthe computational cost (in FLOP count) equivalent to that of a non-relativistic calculation with\nmore than twice as many determinants."}]

Decision letter