Hamiltonian simulation-based quantum-selected configuration interaction for large-scale electronic structure calculations with a quantum computer

Abstract

Quantum-selected configuration interaction (QSCI) is an approach for quantum chemical calculations using current quantum computers. In conventional QSCI, Slater determinants used for the wave function expansion are sampled by iteratively performing approximate wave function preparation and subsequent measurement in the computational basis, and then the subspace Hamiltonian matrix is diagonalized on a classical computer. In this approach, preparation of a high-quality approximate wave function is necessary to accurately compute total energies. Here we propose a Hamiltonian simulation-based QSCI (HSB-QSCI) to avoid this difficulty, by sampling the Slater determinants from quantum states generated by the real-time evolution of approximate wave functions. We provide numerical simulations for the lowest spin-singlet and triplet states of oligoacenes (benzene, naphthalene, and anthracene), phenylene-1,4-dinitrene, and hexa-1,2,3,4,5-pentaene. We found that the HSB-QSCI is applicable not only to molecules where the Hartree–Fock provides a good approximation of the ground state, but also to strongly correlated systems where preparing a high-quality approximate wave function is hard. Hardware demonstrations of the HSB-QSCI are also reported for carbyne molecules expressed by up to 36 qubits, using an IBM quantum processor. The HSB-QSCI captures more than 99.18% of the correlation energies in the active space by considering about 1% of all the Slater determinants in 36 qubit systems, illustrating the ability of the proposed method to efficiently consider important electronic configurations.

---

1 Introduction

Among the diverse topics in the field of quantum computing, quantum chemical calculations of atoms and molecules have attracted attention as a promising application of quantum computers. A method to perform the full-configuration interaction (full-CI) calculation using the quantum phase estimation (QPE) algorithm was proposed in 2005,¹ and proof-of-principle experiments for the full-CI/STO-3G of the H₂ molecule were reported in 2010.^2,3 The quantum circuit for QPE-based full-CI is usually very deep, so highly accurate QPE demonstrations were limited to small models with a few qubits.^4–8 Recently a QPE-type algorithm was demonstrated in models with up to 33 qubits,⁹ however, it is still challenging to compute total energies of larger molecules with chemical precision on quantum computers available today, where chemical precision is defined as an error from the full-CI energy [or complete active space configuration interaction (CAS-CI) energy when the active space approximation is employed] less than 1.0 kcal mol⁻¹.

To reduce the computational load of quantum computers, several quantum–classical hybrid algorithms have been proposed. In 2014, the variational quantum eigensolver (VQE) was proposed for quantum chemical calculations using noisy intermediate-scale quantum (NISQ) devices.¹⁰ In VQE, the quantum state corresponding to an approximate wave function is generated using a parameterized quantum circuit, and the energy expectation value is computed by repeatedly performing state preparation and measurements. The parameters are then optimized to minimize the energy expectation value on a classical computer. These steps are iterated until convergence is reached. VQE has been extensively studied both theoretically and experimentally,¹¹ but it often suffers from challenges related to the high sampling cost required to achieve chemical precision¹² and issues with variational optimization in the presence of barren plateaus.¹³

As an alternative approach to quantum chemical calculations using current quantum devices, the quantum-selected configuration interaction (QSCI) method has been proposed.¹⁴ In the QSCI, quantum computers are used to sample Slater determinants that are important contributors to the ground state wave function. This can be done by running a quantum circuit to prepare an approximate wave function and measuring the quantum state in the computational basis. The subspace Hamiltonian matrix is then constructed from the Slater determinants selected via quantum computation and diagonalized on a classical computer. Note that the concept of the selected CI has been well investigated in classical computation, and various selected CI methods have been developed, such as the CI using a perturbative selection made iteratively (CIPSI),¹⁵ heat-bath CI (HCI),¹⁶ and adaptive sampling CI (ASCI).¹⁷ These methods are often combined with multi-reference perturbation theory.^18,19 In the quantum domain, QSCI deals with the variational optimization of the wave function in the selected subspace, and the effects of the remaining dynamical correlations can be taken into account by combining it with, for example, auxiliary-field quantum Monte Carlo (AFQMC)^20,21 or general multi-configurational quasi-degenerate perturbation theory (GMC-QDPT),²² although the perturbative energy correction from the QSCI is outside the scope of this study.

In the QSCI approach, the quantum computer is anticipated to efficiently sample a polynomial number of important Slater determinants from an exponentially large Hilbert space, and the nature of the quantum states used for sampling determines the accuracy of the calculation. At the stage of preparing an approximate wave function in QSCI, we may employ VQE or the adiabatic state preparation method.¹⁴ To advance the former technique, combining an adaptive strategy for ansatz construction (ADAPT²³) to enhance the efficiency of VQE,²⁴ and using a simpler cost function for VQE parameter optimization²⁵ have also been proposed. However, these approaches may still face challenges in variational optimizations when implemented on quantum hardware. It should also be noted that in the latter approach,²⁵ the VQE wave function with optimized parameters does not necessarily correspond to the minimum with respect to the given Hamiltonian. Another strategy to prepare the quantum state for QSCI is to perform the coupled cluster singles and doubles (CCSD) calculation on a classical computer and embed the CCSD wave function using the local unitary cluster Jastrow (LUCJ) ansatz.²⁶ By using the LUCJ ansatz and introducing an error mitigation technique called self-consistent configuration recovery (SCCR), which is inspired by the structure of chemistry problems, Robledo-Moreno and coworkers reported the QSCI calculations on the IBM superconducting quantum processor and supercomputer “Fugaku” for the nitrogen molecule (N₂) and the iron–sulfur clusters [2Fe–2S] and [4Fe–4S] with up to 77 qubits.²⁷ This approach has also been applied to interaction energy²⁸ and excited state calculations.^29,30 The combination of the QSCI with density matrix embedding theory has also been reported.³¹ While this approach appears to be successful, it is not clear whether it is superior to the method of sampling important Slater determinants from the CCSD wave function on a classical computer.

In this work, we report a Hamiltonian simulation-based quantum state preparation for QSCI (HSB-QSCI), in which important Slater determinants are sampled from the quantum states after the real-time evolutions of initial wave functions.³² Our proposed approach has several advantages over existing methods: (1) a simple approximate wave function, such as the HF wave function, can be used as the initial wave function, and no variational optimization in VQE is required. (2) Different samples can be obtained by changing the duration of the evolution time, and higher-order excitations can be considered even for short-time evolution. (3) The accuracy of the Hamiltonian simulation does not need to be very high. As we will show later, the first-order Trotter decomposition with the time length of a single Trotter step Δt = 1 in atomic unit is sufficient to collect important Slater determinants for the organic molecules under study. In addition, Hamiltonian term truncation based on operator locality in the localized molecular orbital basis can help to reduce the computational cost. (4) Feedback of the QSCI results into quantum computations is possible; that is, Slater determinants found to be important in the QSCI wave function can be added to the initial wave functions. (5) Although short-time evolution can be accurately simulated on a classical computer, long-time evolution poses significant challenges for classical simulation.

Note that methods of sampling basis states for Hamiltonian matrix diagonalization based on Hamiltonian simulation have already been investigated in the framework of quantum Krylov subspace (QKS) algorithms.^33,34 In the QKS method, the Hamiltonian matrix H_kl = 〈Φ_k|H|Φ_l〉 and the overlap matrix S_kl = 〈Φ_k|Φ_l〉 are computed using a quantum computer via a Hadamard test, and the generalized eigenvalue problem Hc = ScE is solved on a classical computer. Here, |Φ_k〉 = U^k|Φ₀〉, where U = e^−iHΔt and |Φ₀〉 is the initial wave function. Since the Taylor expansion of the time evolution operator yields polynomials of H as follows, time evolution is expected to help finding electron configurations that are important for the wave function expansion.

Unlike HSB-QSCI, QKS is applicable to systems in which the correlated wave function is supported by an exponentially large number of Slater determinants. However, the overlap matrix S becomes ill-conditioned when the subspace basis is nearly linearly dependent,^34,35 and the sampling cost in QKS becomes large in order to avoid numerical instability. In contrast, HSB-QSCI requires only measurements in the computational (Pauli-Z) basis and therefore the computational cost on a quantum computer is significantly reduced. HSB-QSCI is also expected to be numerically more stable, because Slater determinants are orthonormal, while subspace Hamiltonian diagonalization on a classical computer can be a bottleneck. In the context of the Monte Carlo simulation, a quantum dynamics Hamiltonian Monte Carlo (qdHMC) method is also proposed,³⁶ in which the proposal step is performed on a quantum computer by sampling after time evolution.

As a proof-of-concept demonstration of the HSB-QSCI method, we first carried out numerical simulations for the spin-singlet ground state of the H₂O molecule, and compared the results obtained using different sampling strategies (randomly sampling Slater determinants preserving spatial symmetry, M_S, and n_elec, where M_S and n_elec represent the spin magnetic quantum number and the number of electrons, respectively, VQE-based QSCI,¹⁴ and HCI). Then, we performed numerical simulations for the spin-singlet ground (S₀) state and the first excited spin-triplet (T₁) state of oligoacenes [benzene (1; 12 qubits), naphthalene (2; 20 qubits), and anthracene (3; 28 qubits)], phenylene-1,4-dinitrene (4; 20 qubits), and hexa-1,2,3,4,5-pentaene (5; 20 qubits) in both planar and twisted geometries. Our numerical simulations showed that the HSB-QSCI can calculate the total energy with chemical precision by sampling the Slater determinants with about 10–20 steps of Hamiltonian simulations. We also report hardware demonstrations of the HSB-QSCI for the S₀ and T₁ states of three carbyne (one-dimensional sp-hybridized carbon atom chain) molecules including 5, octa-1,2,3,4,5,6,7-heptaene (6; 28 qubits), and deca-1,2,3,4,5,6,7,8,9-nonaene (7; 36 qubits) using an IBM quantum processor ibm_kawasaki, with the aid of matrix product operator (MPO)-based classical compression of the quantum circuit for time evolution. The HSB-QSCI energies calculated on ibm_kawasaki agree with the CAS-CI values within 0.6 kcal mol⁻¹ of error in 5. In 6, chemical precision for the total energy was not achieved with the real quantum device, but the HSB-QSCI predicted the singlet–triplet energy gap and the energy difference between planar and twisted geometries in chemical precision. The energy difference between planar and twisted geometries in the S₀ state of 7 was also calculated with an error of 0.29 kcal mol⁻¹. In 7, the HSB-QSCI captures more than 99.18% of correlation energies by considering only about 1% of all the Slater determinants.

2 Theory

In the ab initio molecular orbital theory, the wave function taking into account the electron correlation is expressed by a linear combination of the Slater determinants as follows:

Here, |ψ_HF〉 is the HF determinant and |ψ_n〉 are electronically excited determinants from the HF, and c_n are the corresponding expansion coefficients. In the CI method, the number of Slater determinants included in the wave function expansion increases rapidly with the number of molecular orbitals and electrons, and the excitation order. To calculate the total energy with high accuracy and low computational cost, methods to perform the CI expansion with selected Slater determinants have been investigated.^{15–17,37,38} The QSCI method uses a quantum computer to select important Slater determinants through the approximate wave function preparation and the subsequent computational-basis measurement. In the HSB-QSCI method, important Slater determinants are sampled from the measurement of the quantum state obtained from the time evolution of an approximate wave function.

A schematic view of the HSB-QSCI method is shown in Fig. 1. In the HSB-QSCI, the electronic configurations (Slater determinants) used as the basis for the Hamiltonian matrix are sampled from the time-evolved wave functions {|Φ_k〉}:

|Φ_k〉 = e^−iHkΔt|Φ₀〉. (3)

Fig. 1 Schematic view of the HSB-QSCI algorithm.

To simulate the real time evolution on a quantum computer, the second quantized Hamiltonian given as

is transformed to a qubit Hamiltonianusing a fermion-qubit mapping method. Here, h_pq and g_pqrs in eqn (4) are one- and two-electron integrals, respectively, and a^†_p and a_p are the creation and annihilation operators, respectively, acting on the p-th spin orbital. P_j is a tensor product of Pauli operators referred to as a Pauli string, and w_j is its corresponding coefficient computed from h_pq and g_pqrs. In this work, we used the Jordan–Wigner transformation (JWT)³⁹ as the fermion-qubit mapping, where each qubit stores the occupation number of the corresponding spin orbital: 1 if the spin orbital is occupied, and 0 otherwise. In the JWT, the number of qubits required for wave function mapping is equal to the number of spin orbitals in the active space. The quantum circuit for the time evolution operator U = e^−iHΔt can be constructed using conventional approximate approaches such as the Trotter–Suzuki decomposition^40,41 and qDRIFT,⁴² or prepared through classical optimization⁹ or via variational quantum algorithms.⁴³

The initial wave function |Φ₀〉 should have an overlap with the target electronic state. The HF wave function |ψ_HF〉 is a reasonable choice for |Φ₀〉 in the electronic ground state calculations of typical closed-shell singlet molecules in their equilibrium geometries, but |Φ₀〉 need not be a single Slater determinant. By defining the initial wave function as

where |Ψ_l〉 is the l-th eigenfunction of the Hamiltonian, the quantum state after the time evolution is written as follows:

It is clear that the total evolution time length kΔt controls the magnitude of the interferences between the eigenstates, and the measurements of the quantum states in the computational basis with different k can yield different Slater determinants. The information of the sampled Slater determinants is transferred to the classical computer, and the Hamiltonian matrix with the selected configurations is constructed and then diagonalized to obtain the QSCI wave functions and energies.

The effectiveness of sampling in the HSB-QSCI is mainly controlled by the evolution time length Δt, the number of time steps k, and the initial wave function |Φ₀〉. Note that from eqn (7), we expect that kΔt must be set longer for systems with smaller energy gaps with the excited states. In the HSB-QSCI, we can reconstruct |Φ₀〉 using the information of the QSCI wave function by constructing a multiconfigurational wave function consisting of several Slater determinants that have large contributions in the QSCI wave function in the previous step. Performing Hamiltonian simulations with different initial wave functions and merging the measurement results to perform the QSCI is another option. The availability of such feedback from the subspace Hamiltonian diagonalization part on a classical computer to the state preparation part on a quantum computer is one of the important features of the HSB-QSCI method.

Note that the wave functions are simultaneous eigenfunctions of the Hamiltonian and electron spin S² operator in the non-relativistic regime, but Slater determinants are not always eigenfunctions of the S² operator (spin eigenfunctions). In fact, open-shell Slater determinants with spin-β unpaired electron(s) are not spin eigenfunctions. Here we have assumed that n_α ≥ n_β, where n_α and n_β are the numbers of spin-α and spin-β electrons, respectively. In order to make the QSCI wave function being spin eigenfunctions, we introduced a symmetry completion step before constructing the subspace Hamiltonian matrix. In the symmetry completion step, open-shell Slater determinants that are missing to span the spin eigenfunction are added. For example, in the spin-singlet state calculations, if the bit string corresponding to the Slater determinant |2ααββ0〉 is sampled (2, α, β, and 0 mean that the corresponding molecular orbital is doubly occupied, singly occupied by a spin-α electron, singly occupied by a spin-β electron, and unoccupied, respectively), then we add |2αβαβ0〉, |2αββα0〉, |2ββαα0〉, |2βαβα0〉, and |2βααβ0〉, if they are not sampled. Similarly for the spin-triplet state, if |2αααβ0〉 is measured, then we add |2ααβα0〉, |2αβαα0〉, and |2βααα0〉.

In a case when the Hamiltonian simulations are carried out on a current noisy quantum hardware, the measurements of the quantum state give bit strings that correspond to the Slater determinants with different numbers of electrons. The SCCR proposed recently²⁷ can be adopted to recover the Slater determinants with a correct number of electrons. In the SCCR, if the bit string obtained from the measurement has the wrong number of electrons, bits are flipped to recover the desired number of electrons, and the probability of bit flipping is determined from the occupation number of the molecular orbital calculated from the QSCI wave function of the previous step. This process is iterated until convergence or the iteration step reaches the predefined maximum steps. The SCCR is performed before symmetry completion.

3 Computational conditions

To demonstrate the HSB-QSCI, we performed numerical simulations on a classical computer. The target systems are the S₀ state of the H₂O molecule, the S₀ and T₁ states of oligoacenes (1, 2, and 3), 4 as a representative system of the open-shell singlet ground state, and 5, as shown in Fig. 2. We also performed hardware demonstrations using an IBM Quantum processor for carbyne molecules 5, 6, and 7. For the H₂O molecule, we used the same geometry as in the original QSCI paper.¹⁴ Geometry optimizations of oligoacenes and carbynes in planar geometry were performed at the B3LYP^44,45/6-31G*^46,47 level of theory. In all DFT calculations, we used default settings in Gaussian 16.⁴⁸ The SCF convergence thresholds were set to 10⁻⁸ and 10⁻⁶ for the RMS and maximum density changes, respectively, and 10⁻⁶ Hartree for the energy changes. For geometry optimizations, the convergence criteria were 0.000450 and 0.000300 Hartree Bohr⁻¹ for the maximum and RMS forces, and 0.001800 and 0.001200 Bohr or Radian for the maximum and RMS displacements, respectively. We used an UltraFine grid with a pruned (99 590) grid, having 99 radial shells and 590 angular points per shell. The twisted geometry of carbynes is generated by rotating one of the CH₂ moieties 90 degrees from the equilibrium geometry. The geometry optimization of 4 was carried out at the CASSCF(10e,10o)/6-31G* level. Here, (Ne, Mo) represents the active space consisting of N electrons and M molecular orbitals. For the CASSCF geometry optimization, we used a gradient convergence threshold of 10⁻⁴ Hartree Bohr⁻¹, which is the default setting in the GAMESS-US⁴⁹ software. Convergence is achieved when the largest gradient component is below the specified threshold and the root mean square gradient is less than 1/3 of the threshold. Cartesian coordinates of all molecules are given in the SI.

Fig. 2 Target molecules being studied.

Since the computational cost of numerical simulations of quantum circuits scales exponentially with the number of qubits, we adopted the active space approximation. For the H₂O molecule, we used the (6e,5o) active space, consistent with the QSCI study by Kanno and coworkers.¹⁴ We used the (6e,6o), (10e,10o), and (14e,14o) active spaces for 1, 2, and 3, respectively, which consist of the valence π and π* orbitals. The one- and two-electron integrals of H₂O and oligoacenes are calculated at the RHF/STO-3G⁵⁰ level. The active space of 4 contains the valence π and π* orbitals and the in-plane 2p orbital of the nitrogen atoms, and the size is (10e,10o). The CASSCF(10e,10o)/6-31G* molecular orbitals are used to construct the second quantized Hamiltonian of 4. In the case of carbyne molecules 5, 6, and 7 with (10e,10o), (14e,14o), and (18e,18o) active spaces, respectively, the RHF/6-31G* calculation is performed and then the occupied molecular orbitals are localized using the Pipek–Mezey method.⁵¹ For the virtual orbitals, we first formed the singular value decomposition (SVD) quasi-atomic external orbitals using an SVD with respect to the accurate atomic minimal basis functions, and then formed the ordered external orbitals using exchange integrals (using a keyword EXTLOC = ATMNOS in GAMESS-US software⁴⁹). Localizing the occupied and unoccupied orbitals separately is essential to maintain the invariance of the HF wave function. The active orbitals of all molecules are provided in Fig. S1–S10 in the SI. One- and two-electron atomic orbital integrals are obtained from GAMESS-US, and the AO–MO integral transformations were performed using our own Python3 program. The one- and two-electron MO integrals in the FCIDUMP format are provided in the SI. The qubit Hamiltonian is then generated using the OpenFermion library.⁵² All the DFT calculations were done with the Gaussian 16 package,⁴⁸ and the RHF and CASSCF calculations were carried out with the GAMESS-US software.⁴⁹

As we mentioned in the previous section, the time evolution operation in the HSB-QSCI does not have to be exact, and we can introduce various approximations. In this work, we investigated the effect of Hamiltonian truncation on the sampled Slater determinants in the Hamiltonian simulation and on the QSCI energies. We investigated Hamiltonian truncation based on the locality of the qubit Hamiltonian terms. In this approach, we first generate the localized molecular orbitals (LMOs), and then rearrange the LMOs according to their relative spatial distances. During the qubit Hamiltonian construction using the reordered LMOs, we calculate the locality of the Pauli string (number of Pauli-X, Y, and Z operators in the Pauli string). If the locality of the Pauli string exceeds the threshold, then we exclude it from the qubit Hamiltonian. We investigated this method in 5, keeping in mind to perform the Hamiltonian simulations on a superconducting quantum device with MPO-based classical compression of the time evolution quantum circuit. The procedure of this step is illustrated in Fig. 3.

Fig. 3 Procedure of orbital transformations for the Hamiltonian term truncations based on the locality of the Pauli strings, in the case of 5 in its planar geometry. Red arrows indicate the electron occupancies in the RHF wave function.

To sample the Slater determinants based on the real-time evolution of the initial wave function, we set the evolution time for U to be Δt = 1 in atomic unit and k = 1, 2, …, 10, unless otherwise noted. To construct the quantum circuit for U, we adopt the first-order Trotter decomposition with a single Trotter slice. Magnitude ordering⁵³ is used for the Trotterized term ordering. The quantum circuit simulations were done using the qsim library,⁵⁴ which allows us to use GPGPU. The subspace Hamiltonian matrix construction and diagonalization was done using the PyCI library,⁵⁵ and the SCCR for hardware demonstrations was carried out using the qiskit-addon-sqd library.⁵⁶ It should be noted that the subspace Hamiltonian diagonalization can also be done using qiskit-addon-sqd, but it uses a selected CI kernel_fixed_space subroutine in PySCF⁵⁷ where the spin α and β occupancies are stored separately, and the Slater determinants constructed from all possible combinations of spin α and β strings are used as the basis for the wave function expansion. In this implementation, Slater determinants not sampled by the quantum computer can be included in the subspace Hamiltonian diagonalization, resulting in a quadratic increase in the dimension of the Hamiltonian matrix. This can be seen as a kind of subspace enlargement at the expense of compactness of the QSCI wave function. We also performed a subspace Hamiltonian diagonalization using PySCF, and the results are given in the SI, Section III. Quantum circuit simulations and the subspace Hamiltonian diagonalization in the QSCI part were performed on the NVIDIA DGX H100 system.

4 Results and discussion

Sections 4.1–4.4 are devoted to showing numerical simulations for the S₀ state of H₂O, S₀ and T₁ states of oligoacenes, phenylene-1,4-dinitrene, and hexa-1,2,3,4,5-pentaene, respectively. Section 4.5 gives a hardware demonstration for the S₀ and T₁ states of carbyne molecules.

4.1 The S₀ state of the H₂O molecule

In this subsection, we discuss the performance of the sampling strategies for Slater determinants in the QSCI, using the H₂O molecule. We compare the HSB-QSCI, QSCI with randomly sampled Slater determinants that preserve spatial symmetry, M_S, and n_elec, VQE-based QSCI,¹⁴ and HCI. HSB-QSCI calculations were performed with k = 10, and the number of shots per time step was set to 1 × 10⁴. For the random sampling-based QSCI, we always included the HF configuration in the subspace Hamiltonian, and excited determinants belonging to the same irreducible representation of the S₀ state and having the correct M_S and n_elec. The number of excited determinants satisfying these conditions is 99. We performed five independent simulations. The HCI is performed using the PyCI package.⁵⁵ The results are summarized in Fig. 4. In HSB-QSCI with k = 1, the energy error is already less than 6 × 10⁻⁴ Hartree (0.36 kcal mol⁻¹). Chemical precision is achieved with only one time evolution step, considering only 20 Slater determinants. This is in contrast to the VQE-based QSCI,¹⁴ which requires about 50 VQE iterations and 32 Slater determinants to achieve an energy error below 1 × 10⁻³ Hartree. Note that they also reported that 16 Slater determinants are enough to achieve chemical precision, after VQE convergence. It is worth noting that the probability of sampling the HF configuration in the first time step of HSB-QSCI is approximately 0.942. In contrast, in the perfect state preparation limit of VQE-based sampling, this probability is calculated to be 0.978. This indicates that HSB-QSCI has a higher chance of sampling excited determinants than the VQE-based approach. Of course, the quantum circuit is deeper in HSB-QSCI than in the VQE-based one. However, improvements in the overall number of required shots and the compactness of the QSCI wave function are significant. HCI yields slightly lower energy than HSB-QSCI with 20 Slater determinants, ΔE_HCI–CAS-CI = 5 × 10⁻⁵ Hartree = 0.032 kcal mol⁻¹. However, as the number of Slater determinants increases, the difference in energy errors between HSB-QSCI and HCI diminishes, and the two methods give identical energies when 28 Slater determinants are used. It is important to note that HSB-QSCI yields energies comparable to HCI, although starting from an approximate wave function in HSB-QSCI may lead to sampling Slater determinants that are more relevant to excited states than to the ground state. The comparison with the random sampling-based implementation is even more illuminating. In this case, the energy error remains large even when 80% of the Slater determinants are included in the subspace Hamiltonian diagonalization.

Fig. 4 HSB-QSCI, random sampling-based QSCI, and HCI results of the H₂O molecule.

4.2 The S₀ and T₁ states of oligoacenes

There is no doubt that aromatic rings are one of the most important molecular skeletons in chemistry, and the study of the electronic structures of oligoacenes is very important. It is known that the zigzag edge of graphene fragments exhibits strong open-shell characters,⁵⁸ and the contribution of the HF electronic configuration to the full-CI wave function decreases for larger oligoacenes. Electronic excited states of oligoacenes are also of interest because of their potential for various applications such as singlet fission photovoltaics.⁵⁹

The results of the HSB-QSCI simulations of oligoacenes are shown in Fig. 5. The RHF and ROHF-like single configurational wave functions are used as the starting wave functions for the Hamiltonian simulations of the S₀ and T₁ states, respectively, and the number of shots for each time step was set as 1 × 10⁵. The number of Slater determinants in the CAS-CI wave function is 104 (S₀) and 61 (T₁) for 1, 15 912 (S₀) and 11 076 (T₁) for 2, and 2 945 056 (S₀) and 2 255 121 (T₁) for 3, in the D_2h point group. Our numerical simulations revealed that the Slater determinants that are important to describe the target electronic states are efficiently sampled from the Hamiltonian simulations. In 1, the system size is very small, and all possible Slater determinants are sampled up to the fifth steps. In 2, chemical precision (ΔE < 1.0 kcal mol⁻¹) was achieved with three steps of the Hamiltonian simulations. In 3, we found that 10 steps of time evolution with 1 × 10⁵ shots are not enough to achieve chemical precision. We expect that the number of important Slater determinants increases with [scr O, script letter O](M⁴), which is the same as the scaling of Hamiltonian terms, when the target state wave function is well approximated by the single determinant as in the ground state of oligoacenes and the energy gap between the ground and the first excited states is unchanged. Considering the fact that the energy gap is smaller in 3 than in 2, we need more than 3.81 times more shots than in the calculation of 2 to obtain the HSB-QSCI energy of 3 with the same accuracy. We have also studied the dependence of the HSB-QSCI energy on the number of shots and the number of time steps in 3. By increasing the number of shots from 1 × 10⁵ and simulating more time steps, the number of Slater determinants sampled from the Hamiltonian simulation increases, and the HSB-QSCI energy decreases systematically (see Fig. S14 in the SI). Although we need more steps and shots to achieve chemical precision for larger systems, it should be noted that up to 7-electron excited determinants were sampled in a Hamiltonian simulation with kΔt = 1 in 3. It is known that higher-order excitations are important to compute the energy in a quantitative manner, especially in strongly correlated systems.¹⁹ The ability of Hamiltonian simulations to capture higher-order excitations is remarkable and is an important feature for applications to large and strongly correlated systems. These results illustrate the ability of the Hamiltonian simulation to sample a polynomial number of important Slater determinants from exponentially large Hilbert space. Note that a similar trend has also been reported by other groups.^60,61 When the kernel_fixed_space subroutine in PySCF is used for subspace Hamiltonian diagonalization, chemical precision was achieved with only three steps of time evolutions with 1 × 10⁵ shots in 3. However, the HSB-QSCI wave function considers about 57% and 46% of the Slater determinants, hence the HSB-QSCI wave function is less compact (see Fig. S11 in the SI).

Fig. 5 HSB-QSCI results of oligoacenes. The number of Slater determinants included in the Hamiltonian diagonalization is given in (a), (b) and (c) for 1, 2, and 3, respectively. Red and blue indicate the S₀ and T₁ states, respectively. Dotted lines represent the number of Slater determinants in the CAS-CI wave function. The difference of the HSB-QSCI energy from the CAS-CI values in units of kcal mol⁻¹ are given in (d)–(f) on the logarithmic scale.

4.3 The S₀ and T₁ states of phenylene-1,4-dinitrene

4 has two major resonance structures, the diradical with a quinonoid skeleton as shown in Fig. 2, and the dinitrene with an aromatic ring. Its electronic ground state is an open-shell spin-singlet state with two unpaired electrons in in-plane 2p orbitals of the nitrogen atoms. The thermally excited spin-triplet state of 4 has been observed by ESR spectroscopy,^62–64 and the singlet–triplet energy gap is determined to be −0.82 kcal mol⁻¹ from the Curie analysis of the ESR signal.⁶⁵ The experimentally determined zero-field splitting parameter in the excited spin-triplet state is |D| = 0.171 cm⁻¹,⁶² which is about 10 times larger than the value estimated from a point-dipole approximation. This deviation was explained by the lack of electron correlation effect in the point-dipole approximation,⁶⁶ insisting that sophisticated consideration of electronic correlation is essential to describe its electronic structure.

In the CASSCF optimized orbital basis, the CAS-CI wave function of the S₀ state is mainly described by the HF configuration (|ψ_HF〉) and the HOMO–LUMO two-electron excited configuration from the HF configuration (|ψ_2e〉) (see Table S7 in the SI). In the HSB-QSCI simulations of the S₀ state, we performed Hamiltonian simulations with three different initial wave functions, |ψ_HF〉, |ψ_2e〉, and the CAS-CI(2e,2o) wave function |Ψ_CAS22〉 = 0.7138|ψ_HF〉 − 0.7003|ψ_2e〉, to investigate the dependence of the initial wave function on the energies and the convergence behavior of the HSB-QSCI method. The QSCI simulations were performed under four different conditions: (1) use the |ψ_HF〉 results with 1 × 10⁵ shots for each time step, (2) use the |ψ_2e〉 with 1 × 10⁵ shots, (3) use the samples from |Ψ_CAS22〉 with 1 × 10⁵ shots, and (4) merge the samples from |ψ_HF〉 and |ψ_2e〉 with 5 × 10⁴ shots for each step. For the T₁ state, we used the ROHF-like single determinant as the initial wave function for the Hamiltonian simulation. The HSB-QSCI results are summarized in Fig. 6. Our numerical simulations suggest that almost the same number of Slater determinants are sampled when |ψ_HF〉 and |ψ_2e〉 are used. Interestingly, the number of Slater determinants sampled from the Hamiltonian simulation is smaller when |Ψ_CAS22〉 is used as the initial wave function. The accuracy of the HSB-QSCI energy with |Ψ_CAS22〉 is about an order of magnitude worse than the HF-reference one, but still chemical precision was achieved in 2 steps of time evolution. Thus, using the CAS-CI wave function with a smaller active space can help generate a more compact HSB-QSCI wave function for strongly correlated systems. However, preparing the CAS-CI reference state itself may become a bottleneck when the number of quasi-degenerate orbitals is large. One possible solution to this challenge is to adopt a technique for constructing a multiconfigurational wave function without performing post-HF calculations.⁶⁷ The singlet–triplet energy gap is calculated to be −1.11 kcal mol⁻¹, in good agreement with the CAS-CI value (−1.06 kcal mol⁻¹) and the experimental one.⁶⁵ Using PySCF for the subspace Hamiltonian diagonalization slightly increases the number of Slater determinants and yields lower HSB-QSCI energies. The trends of the initial wave function dependence are similar between PyCI- and PySCF-based implementations (see Fig. S12 in the SI).

Fig. 6 HSB-QSCI results of 4. (a) The number of Slater determinants sampled from Hamiltonian simulations. The red and blue dotted lines indicate the number of Slater determinants in the CAS-CI wave function. (b) The difference of the HSB-QSCI energy from the CAS-CI values in units of kcal mol⁻¹, on the logarithmic scale.

4.4 The S₀ and T₁ states of hexa-1,2,3,4,5-pentaene

Our last example of numerical simulations is the S₀ and T₁ states of 5. From the DFT calculations, the carbyne with =CH₂ termination is predicted to have a helical π conjugation in the non-planar geometries.^68,69 In the S₀ state, the planar structure is the energy minimum, but in the T₁ state, the planar geometry is a saddle point and the twisted geometry becomes stable. As a result, the singlet–triplet energy gap strongly depends on the dihedral angle between two =CH₂ terminations. In this study, we calculated the S₀ and T₁ states of the geometries with the dihedral angles of 0° and 90°.

Because carbyne is a one-dimensional molecule, using LMOs and adopting Hamiltonian truncation by maximum locality may be a good option to reduce the computational cost of Hamiltonian simulation. In this work, we investigated Hamiltonian term truncation by operator locality of the Pauli string in the qubit Hamiltonian. By reordering the LMOs by relative distances before constructing a qubit Hamiltonian, operator locality-based truncation is roughly equivalent to Hamiltonian term truncation based on spatial distances. To assess the effect of Hamiltonian truncation with the locality, we first calculated the number of Hamiltonian terms in the second quantized Hamiltonian and the fidelity of the ground state wave function, |〈Ψ₀(truncated)|Ψ₀(full)〉|², where |Ψ₀(truncated)〉 and |Ψ₀(full)〉 are the ground state wave functions with the truncated and untruncated Hamiltonians, respectively. The results are summarized in Tables 1 and 2 for planar and twisted geometries, respectively.

Table 1 Number of terms in the second quantized Hamiltonian and fidelity of the electronic ground state wave function of 5 in planar geometry

Maximum locality m	#Terms	Fidelity
20	20 072	1.0
18	19 676	0.999816
16	18 596	0.998764
14	17 592	0.998571
12	16 248	0.996394
10	12 708	0.992327
8	7836	0.946569
6	4840	0.941506
4	3048	0.912734

---

Table 2 Number of terms in the second quantized Hamiltonian and fidelity of the electronic ground state wave function of 5 in twisted geometry

Maximum locality m	#Terms	Fidelity
20	40 200	1.0
18	39 772	0.999576
16	38 244	0.997395
14	35 232	0.996860
12	30 640	0.995021
10	24 660	0.987597
8	17 772	0.898829
6	10 744	0.891525
4	4632	0.843284

---

Note that the number of Hamiltonian terms is larger in the twisted geometry than in the planar one, which is due to the difference in the spatial distribution of π orbitals. In planar geometry, the one-electron MO integrals between the nearest neighbor π orbitals (e.g., π₁ and π₂ in Fig. 3) are zero, while in twisted geometry, they are non-zero due to helical π conjugations. As can be seen from Tables 1 and 2, the number of Hamiltonian terms is reduced by about 40% when considering up to 10-local terms, with the fidelity of the ground state wave function being about 0.99. The rapid decrease of the fidelity for the maximum locality smaller than 10 can be explained by the fact that the π_x–π_x and π_y–π_y interactions (here we assumed that the one-dimensional carbon atom chain is parallel to the z-axis) with the next nearest neighbor C=C bond (for example, π₁–π₃ in Fig. 3) described as 10-local operators in the qubit Hamiltonian.

The results of the HSB-QSCI simulations are shown in Fig. 7. The number of shots for each time step is 1 × 10⁵. In the T₁ state calculations, we set the starting wave function to carry two unpaired electrons in the central π bonds (π₃ and in Fig. 3). It is clear that the number of Slater determinants sampled from the Hamiltonian simulations and the convergence behavior of the HSB-QSCI energies are almost the same for m ≥ 10, where m is the maximum locality of the qubit Hamiltonian terms. These results are consistent with the trend of the S₀ state fidelity value. The singlet–triplet energy gap ΔE_S–T = E_S − E_T calculated using the real-time evolution with the untruncated Hamiltonian (m = 20) is −53.96 and −20.33 kcal mol⁻¹ for planar and twisted geometries, respectively, and the difference from the CAS-CI singlet–triplet energy gap is less than 0.1 kcal mol⁻¹ (ΔE_S–T(CAS-CI) = −53.93 and −20.26 kcal mol⁻¹ for planar and twisted geometry, respectively). The number of sampled Slater determinants is larger for twisted geometry than for planar geometry, reflecting the spatial distribution of the π orbitals, as discussed above. When the maximum locality is set to be m = 6 or 8, the error in the HSB-QSCI energy becomes larger. However, our numerical simulations revealed that chemical precision can be achieved even with the maximum locality m = 6, by increasing the number of shots to more than 3 × 10⁶ and setting the number of maximum time steps to k = 20 (see Fig. S15 in the SI). In this case, we observed an abrupt improvement in energy convergence. Since we did not observe such behavior when the untruncated Hamiltonian was used for the time evolution, we suspect that Hamiltonian truncation is responsible for this. When a Hamiltonian truncated to maximum locality is used for time evolution, the probability of measuring excited determinants with long-range excitations becomes smaller, because long-range interaction terms are discarded by the truncation. Consequently, the use of the truncated Hamiltonian in the time evolution operator causes an imbalance in the description of short-range and long-range interactions, which is likely responsible for the observed behavior.

Fig. 7 HSB-QSCI results of 5 in its planar and twisted geometries with different maximum locality values m. The number of Slater determinants included in the Hamiltonian diagonalization is given in (a)–(d). The red and blue dotted lines represent the number of Slater determinants in the CAS-CI wave function of the S₀ and T₁ states, respectively. The difference of the HSB-QSCI energy from the CAS-CI values in units of kcal mol⁻¹ are given in (e)–(h) on the logarithmic scale.

Note that when using the PySCF kernel_fixed_space subroutine for subspace Hamiltonian diagonalization, chemical precision was achieved with up to five steps of time evolution with 10⁵ shots (see Fig. S13 in the SI), but this is a consequence of the quadratic increase of the Slater determinants in PySCF.

It should be emphasized that the ROHF-like single configurational wave function in the LMO basis used in the T₁ state calculation has a rather small overlap with the CAS-CI wave function, because the LMOs were generated according to the occupation numbers of the RHF wave function for the S₀ state. The squared overlap values |〈Φ₀|Ψ_CAS-CI〉|² in the T₁ state are 0.1691 and 0.1240 for planar and twisted geometries, respectively (see Tables S10 and S12 in the SI for the CAS-CI wave function). The fact that the HSB-QSCI provides accurate energy with such a small overlap is promising because the overlap between the HF and the full-CI wave functions decreases with system size, and the preparation of a sophisticated approximate wave function becomes challenging for larger molecules.

4.5 Hardware demonstration of the HSB-QSCI calculations of the S₀ and T₁ states of carbyne molecules

Since the Hamiltonian term truncation based on the locality of the qubit Hamiltonian terms works excellently in 5, we performed the HSB-QSCI calculations of 5, 6, and 7 on an IBM Quantum Eagle processor, namely ibm_kawasaki. To reduce the quantum circuit depth of the Hamiltonian simulation, we performed the MPO-based classical compression of the quantum circuit. This approach was demonstrated in our previous study on the quantum phase difference estimation-based energy gap calculations.⁹ The details of the MPO-based quantum circuit optimization are given in ref. 9. In this study, we used a Hamiltonian truncated by the maximum locality (m = 10), and the 10-layer brick wall type quantum circuit is generated for the time evolution operator of Δt = 1 in atomic unit. The number of sweeps on the optimizations was set to 10 000. In constructing the MPOs, the second-order Trotter decomposition with a single Trotter slice was employed, and the singular value cutoff for the SVD was set to 10⁻⁴. The order of terms in the Trotter decomposition is set as follows. First, we generate the truncated second quantized Hamiltonian by checking the one- and two-electron integrals in the lexicographic order: q after p in h_pq, and p, r, q, s, in that order in g_pqrs. Then the truncated second-quantized Hamiltonian is transformed into a qubit Hamiltonian while retaining the order of the terms, and the Trotterized time evolution operator is constructed subsequently.

We first checked the Hamming weight distribution of the bit strings obtained in the measurements and compared it to the uniform distribution. The Hamming weight distributions at kΔt = 10 are plotted in Fig. 8. We confirmed that the Hamming weight distributions of the bit strings obtained from the Hamiltonian simulation with kΔt = 10 are significantly different from the uniform distribution. However, due to hardware noise, only 11–23% of the measured bit strings have correct Hamming weights (10, 14, and 18 for 5, 6, and 7, respectively).

Fig. 8 Hamming weight distribution of the bit strings obtained from measurements after Hamiltonian simulation of carbyne molecules with kΔt = 10, using ibm_kawasaki. (a) 5, (b) 6, and (c) 7.

The HSB-QSCI results of 5 are summarized in Fig. 9. Here, the Hamiltonian simulations were carried out with k = 10, and the number of shots was set from 1 × 10⁵ to 1 × 10⁶. The depth of the quantum circuit after transpilation with optimization level = 3 in qiskit is 21, regardless of the evolution time steps. The SCCR²⁷ was adopted to recover configurations with a correct number of electrons, before constructing the subspace Hamiltonian matrix. In the first five to eight steps of time evolution, the number of Slater determinants sampled increases rapidly, but thereafter the number of additionally sampled Slater determinants becomes gradual. When the HSB-QSCI is performed up to ten steps of time evolution and 1 × 10⁵ shots for each time step, the energy error is about 2 kcal mol⁻¹. Increasing the number of shots further improves the energy, and the error in the total energy with 1 × 10⁶ shots were 0.37, 0.56, 0.33, and 0.43 kcal mol⁻¹ for planar (S₀), planar (T₁), twisted (S₀), and twisted (T₁), respectively, achieving chemical precision.

Fig. 9 HSB-QSCI results of 5 obtained with a Hamiltonian truncated by the maximum locality m = 10 and ibm_kawasaki, with SCCR. The number of Slater determinants included in the Hamiltonian diagonalization are given in (a)–(d), and the difference of the HSB-QSCI energy from the CAS-CI values in units of kcal mol⁻¹ are given in (e)–(h) on the logarithmic scale.

The dimensions of the Hamiltonian matrix and the total energies of the HSB-QSCI and the CAS-CI methods are summarized in Table 3. We used k = 10 and performed 3 × 10⁶ shots for each time step in the HSB-QSCI for 6 and 7. In 6, the HSB-QSCI considered about 50–64% of the entire Slater determinant, and it calculated the total energy with an error of about 2.0–2.9 kcal mol⁻¹. Although chemical precision for the total energy was not achieved, the HSB-QSCI provides the singlet–triplet energy gap with an error of 0.93 and 0.12 kcal mol⁻¹ for planar and twisted geometries, respectively, of 6. The energy difference between the planar and twisted geometries is also calculated with an error of about 0.72 and 0.09 kcal mol⁻¹ for the S₀ and T₁ states, respectively. These results exemplify the ability of the HSB-QSCI to accurately calculate energy differences between two electronic states and two geometries. We expect that the total energy can be further improved by increasing the number of shots in the Hamiltonian simulation. In 7, the error of the HSB-QSCI energy is about 14 kcal mol⁻¹ for the S₀ state, and it increases to 23.36 and 25.79 kcal mol⁻¹ for planar and twisted geometries, respectively, of the T₁ state. Again, chemical precision for the total energy was not achieved in 7. In addition, the HSB-QSCI overestimated the singlet–triplet energy gap, about 9–12 kcal mol⁻¹. As we pointed out above, the CAS-CI wave function of the T₁ state in the localized orbital basis shows an inherent multiconfigurational character, and the single Slater determinant used as the initial wave function of the time evolution has small overlap with the CAS-CI wave function. As a result, it is expected that the number of time steps or the number of shots will have to be increased to account for more Slater determinants in order to accurately describe the electronic structure of the T₁ state. Nevertheless, it is worth noting that the HSB-QSCI is able to compute the energy difference between the planar and twisted geometries of the S₀ state with chemical precision. Note that when PySCF is used for the subspace Hamiltonian diagonalization, the HSB-QSCI provided the total energy of the S₀ state with about 0.3 kcal mol⁻¹ of an error, by considering 46–48% of the Slater determinants, but it failed to compute the total energy of the T₁ state with chemical precision (see Table S13 in the SI for details).

Table 3 The HF energy, number of Slater determinants and total energies of the CAS-CI, and the number of selected Slater determinants and the error of the HSB-QSCI energy calculated using ibm_kawasaki

System	HF	CAS-CI		HSB-QSCI
	E/Hartree	#Dets	E/Hartree	%#Dets^a	ΔE/kcal mol^−1	%ECorr^b
a Percentage of the number of Slater determinants calculated as . b Percentage of the correlation energy in the active space considered in the HSB-QSCI, calculated as .
5 (planar, S0)	−226.4832222338	63 504	−229.5024030030	82.21	0.37	99.98
5 (planar, T1)	−226.4299826264	44 100	−229.4164659619	90.91	0.56	99.97
5 (twisted, S0)	−226.4033087159	63 504	−229.4452184415	85.96	0.33	99.98
5 (twisted, T1)	−226.4440971899	44 100	−229.4129399110	92.84	0.43	99.99
6 (planar, S0)	−301.1974823957	11 778 624	−305.2359857689	52.20	2.00	99.92
6 (planar, T1)	−301.1536301615	9 018 009	−305.1666727092	63.80	2.93	99.88
6 (twisted, S0)	−301.1365655879	11 778 624	−305.1941677868	49.57	2.71	99.89
6 (twisted, T1)	−301.1634490372	9 018 009	−305.1698864405	60.33	2.84	99.89
7 (planar, S0)	−375.9128437860	2 363 904 400	−380.9720622554	0.95	14.26	99.55
7 (planar, T1)	−375.8222817220	1 914 762 564	−380.9133659533	0.77	23.36	99.27
7 (twisted, S0)	−375.8635891161	2 363 904 400	−380.9394312061	1.09	13.97	99.56
7 (twisted, T1)	−375.8817559215	1 914 762 564	−380.9186151034	0.87	25.79	99.18

---

Table 3 also summarizes the percentages of the number of Slater determinants and the correlation energies covered by the HSB-QSCI method in 5–7. Here, we calculated the correlation energies as the energy difference between HF (RHF and ROHF for the S₀ and T₁ state, respectively, in the canonical orbital basis) and the HSB-QSCI energies. It is worth emphasizing that the HSB-QSCI covers more than 99.18% of the correlation energies in the active space in 7, by considering only ca. 1% of all the Slater determinants.

5 Conclusion

In this study, we proposed an HSB-QSCI, which uses a quantum computer to perform Hamiltonian simulation and subsequent measurement of the quantum state in the computational basis to sample the important Slater determinants used to span the subspace Hamiltonian to be diagonalized on a classical computer. Compared to the reported VQE-based state preparation approaches for the QSCI,^14,24,25 the HSB-QSCI is free from variational optimizations of the quantum states, and it works with the simple initial wave function like the HF. Since long-time evolution is in general difficult to simulate on a classical computer, we expect this approach to become powerful as the system size increases. The fact that higher-order excitations can be incorporated even in a short time evolution also makes HSB-QSCI suitable for simulation of large and strongly correlated systems. Proof-of-concept simulations were performed for the S₀ and T₁ states of organic molecules 1–5, and we achieved chemical precision in all molecules studied. The quantum hardware demonstrations of the HSB-QSCI for the S₀ and T₁ states of three carbyne molecules 5–7 were also reported using ibm_kawasaki with 20, 28, and 36 qubits, respectively, with the aid of MPO-based classical optimization of the quantum circuit for the time evolution operator. The HSB-QSCI achieved chemical precision for the total energy in 5, for the singlet–triplet energy difference in 6, and for the energy difference between planar and twisted geometries of the S₀ state in 7. In 7, by considering about 1% of the Slater determinants, the HSB-QSCI captured more than 99.18% of correlation energies in the active space. The largest system studied in this work is 7 with 36 qubits, which is the largest system capable of performing the CAS-CI calculation in a typical computational environment without a supercomputer.^57,70 We expect that the HSB-QSCI can handle the system larger than 36 qubits, by introducing batch-based implementations in the Hamiltonian matrix diagonalization part, and calculating the variance of the energy expectation value ΔH = 〈Ψ|H²|Ψ〉 − 〈Ψ|H|Ψ〉² as explored in ref. 27. Note that, as discussed by Reinholdt and coworkers in their recent paper,⁷¹ the QSCI wave function is in some cases less compact than that obtained by classical heuristics. From eqn (7), the probability that a particular Slater determinant is measured in the HSB-QSCI depends not only on the CI coefficient of the ground state, but also on the gap to the excited state and the CI coefficient of the excited state, and screening out unimportant Slater determinants from the subspace is very important in the application to larger systems. Applying the auxiliary-field quantum Monte Carlo (AFQMC) method with the QSCI wave function as the input can be another direction towards larger-scale quantum chemical calculations with a quantum computer, as suggested in recent studies.^20,21 Application of the HSB-QSCI to larger systems is currently in progress.

Author contributions

K. Sugisaki and N. Yamamoto planned and conducted the project. K. Sugisaki and S. Kanno developed the simulation programs and performed numerical calculations. K. Sugisaki carried out hardware runs. All the authors discussed the results and wrote the paper.

Conflicts of interest

The author has no conflicts to disclose.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The data supporting this article have been included as part of the SI. Supplementary information: CAS-CI active space, CAS-CI wave function, HSB-QSCI results calculated using kernel_fixed_space subroutine in PySCF, shot number dependence on the HSB-QSCI simulations of 3, shot number dependence on the HSB-QSCI simulations of 5 with the maximum locality m = 6, Cartesian coordinates in the xyz format, and one- and two-electron integral files in the FCIDUMP format. See DOI: https://doi.org/10.1039/d5cp02202a.

Acknowledgements

This work was supported by the Quantum-LEAP Flagship Program (Grant No. JPMXS0120319794 and JPMXS0118067285) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. K. S. acknowledges support from Center of Innovations for Sustainable Quantum AI (JPMJPF2221) from Japan Science and Technology Agency (JST), Japan and Grants-in-Aid for Scientific Research C (21K03407) and for Transformative Research Area B (23H03819) from Japan Society for the Promotion of Science (JSPS), Japan. A part of this work was performed for the Council for Science, Technology and Innovation (CSTI), Cross-ministerial Strategic Innovation Promotion Program (SIP), “Promoting the application of advanced quantum technology platforms to social issues” (Funding agency: QST). The authors thank Takashi Abe, Yoshiharu Mori, and Hajime Nakamura for useful discussions. Part of the calculations was performed on the Mitsubishi Chemical Corporation (MCC) high-performance computer (HPC) system “NAYUTA”, where “NAYUTA” is a nickname for MCC HPC and is not a product or service name of MCC. We acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors, and do not reflect the official policy or position of IBM or the IBM Quantum team.

Notes and references

1. A. Aspuru-Guzik, A. D. Dutoi, P. J. Love and M. Head-Gordon, Science, 2005, 309, 1704–1707.
2. B. P. Lanyon, J. D. Whitfield, G. G. Gillett, M. E. Goggin, M. P. Almeida, I. Kassal, J. D. Biamonte, M. Mohseni, B. J. Powell, M. Barbieri, A. Aspuru-Guzik and A. G. White, Nat. Chem., 2010, 2, 106–111.
3. J. Du, N. Xu, X. Peng, P. Wang, S. Wu and D. Lu, Phys. Rev. Lett., 2010, 104, 030502.
4. Y. Wang, F. Dolde, J. Biamonte, R. Babbush, V. Bergholm, S. Yang, I. Jakobi, P. Neumann, A. Aspuru-Guzik, J. D. Whitfield and J. Wrachtrup, ACS Nano, 2015, 9, 7769–7774.
5. P. J. J. O'Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik and J. M. Martinis, Phys. Rev. X, 2016, 6, 031007.
6. R. Santagati, J. Wang, A. A. Gentile, S. Paesani, N. Wiebe, J. R. McClean, S. Morley-Short, P. J. Shatbolot, D. Bonneau, J. W. Silverstone, D. P. Tew, X. Zhou, J. L. O'Brien and M. G. Thompson, Sci. Adv., 2018, 4, eaap9646.
7. N. S. Blunt, L. Caune, R. Izsak, E. T. Campbell and N. Holzmann, PRX Quantum, 2023, 4, 040341.
8. K. Yamamoto, S. Duffield, Y. Kikuchi and D. Muñoz Ramo, Phys. Rev. Res., 2024, 6, 013221.
9. S. Kanno, K. Sugisaki, H. Nakamura, H. Yamauchi, R. Sakuma, T. Kobayashi, Q. Gao and N. Yamamoto, Proc. Natl. Acad. Sci. U. S. A., 2025, 122, e2425026122.
10. A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik and J. L. O'Brien, Nat. Commun., 2014, 5, 4213.
11. J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth and J. Tennyson, Phys. Rep., 2022, 986, 1–128.
12. J. F. Gonthier, M. D. Radin, C. Buda, E. J. Doskocil, C. M. Abuan and J. Romero, Phys. Rev. Res., 2022, 4, 033154.
13. J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush and H. Neven, Nat. Commun., 2018, 9, 4812.
14. K. Kanno, M. Kohda, R. Imai, S. Koh, K. Mitarai, W. Mizukami and Y. O. Nakagawa, Quantum-selected confiugration interaction: classical diagonalization of Hamiltonians in subspaces selected by quantum computers, arXiv, 2023, preprint, arXiv:2302.11320, DOI: 10.48550/arXiv.2302.11320.
15. B. Huron, J. P. Malrieu and P. Rancurel, J. Chem. Phys., 1973, 58, 5745–5759.
16. A. A. Holmes, N. M. Tubman and C. J. Umrigar, J. Chem. Theory Comput., 2016, 12, 3674–3680.
17. N. M. Tubman, C. D. Freeman, D. S. Levine, M. Head-Gordon and K. B. Whaley, J. Chem. Theory Comput., 2020, 16, 2139–2159.
18. Y. Garniron, A. Scemama, E. Giner, M. Caffarel and P.-F. Loos, J. Chem. Phys., 2018, 149, 064103.
19. S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi and C. J. Umrigar, J. Chem. Theory Comput., 2017, 13, 1595–1604.
20. Y. Yoshida, L. Erhart, T. Murokoshi, R. Nakagawa, C. Mori, T. Miyanaga, T. Mori and W. Mizukami, Auxiliary-field quantum Monte Carlo method with quantum selected configuration interaction, arXiv, 2025, preprint, arXiv:2502.21081, DOI: 10.48550/arXiv.2502.21081.
21. D. Danilov, J. Robledo-Moreno, K. J. Sung, M. Motta and J. Shee, Enhancing the accuracy and efficiency of sample-based quantum diagonalization with phaseless auxiliary-field quantum Monte Carlo, arXiv, 2025, preprint, arXiv:2503.05967 DOI:10.48550/arXiv.2503.05967.
22. S. Shirai, S.-Y. Tseng, H. Iwakiri, T. Horiba, H. Hirai and S. Koh, ACS Omega, 2025, 10, 39736–39750.
23. H. R. Grimsley, S. E. Economou, E. Barnes and N. J. Mayhall, Nat. Commun., 2019, 10, 3007.
24. Y. O. Nakagawa, M. Kamoshita, W. Mizukami, S. Sudo and Y. Ohnishi, J. Chem. Theory Comput., 2024, 20, 10817–10825.
25. L. Nützel, A. Gresch, L. Hehn, L. Marti, R. Freund, A. Steiner, C. D. Marciniak, T. Eckstein, N. Stockinger, S. Wolf, T. Monz, M. Kühn and M. J. Hartmann, Quantum Sci. Technol., 2025, 10, 015066.
26. M. Motta, K. J. Sung, K. B. Whaley, M. Head-Gordon and J. Shee, Chem. Sci., 2023, 14, 11213–11227.
27. J. Robledo-Moreno, M. Motta, H. Haas, A. Javadi-Abhari, P. Jurcevic, W. Kirby, S. Martiel, K. Sharma, S. Sharma, T. Shirakawa, I. Sitdikov, R.-Y. Sun, K. J. Sung, M. Takita, M. C. Tran, S. Yunoki and A. Mezzacapo, Sci. Adv., 2025, 11, eadu9991.
28. D. Kaliakin, A. Shajan, J. Robledo Moreno, Z. Li, A. Mitra, M. Motta, C. Johnson, A. A. Saki, S. Das, I. Sitdikov, A. Mezzacapo and K. M. Merz Jr., Accurate quantum-centric simulations of supramolecular interactions, arXiv, 2024, preprint, arXiv:2410.09209, DOI: 10.48550/arXiv.2410.09209.
29. S. Barison, J. Robledo Moreno and M. Motta, Quantum Sci. Technol., 2025, 10, 025034.
30. I. Liepuoniute, K. D. Doney, J. Robledo Moreno, J. A. Job, W. S. Friend and G. O. Jones, J. Chem. Theory Comput., 2025, 21, 5062–5070.
31. A. Shajan, D. Kaliakin, A. Mitra, J. Robledo Moreno, Z. Li, M. Motta, C. J. A. A. Saki, S. Das, I. Sitdikov, A. Mezzacapo and K. M. Merz Jr., J. Chem. Theory Comput., 2025, 21, 6801–6810.
32. After posting this manuscript to arXiv, Mikkelsen and Nakagawa, and Yu and coworkers independently submitted their own manuscripts to arXiv, presenting the same idea explored in this study, M. Mikkelsen and Y. O. Nakagawa, arXiv, 2024, preprint, arXiv:2412.13839 DOI:10.48550/arXiv.2412.13839; J. Yu, et al., arXiv, 2025, preprint, arXiv:2501.09702 DOI:10.48550/arXiv.2501.09702.
33. C. L. Cortes and S. K. Gray, Phys. Rev. A, 2022, 105, 022417.
34. Z. Zhang, A. Wang, X. Xu and Y. Li, Quantum, 2023, 8, 1438.
35. G. Lee, S. Choi, J. Huh and A. F. Izmaylov, Digital Discovery, 2025, 4, 954–969.
36. O. Lockwood, P. Weiss, F. Aronshtein and G. Verdon, Phys. Rev. Res., 2024, 6, 033142.
37. J. B. Schriber and F. A. Evangelista, J. Chem. Phys., 2016, 144, 161106.
38. P. M. Zimmerman, J. Chem. Phys., 2017, 146, 104102.
39. P. Jordan and E. Wigner, Z. Phys., 1928, 47, 631–651.
40. H. F. Trotter, Proc. Am. Math. Soc., 1959, 10, 545–551.
41. M. Suzuki, Prog. Theor. Phys., 1976, 56, 1454–1469.
42. Y. Ouyang, D. R. White and E. T. Campbell, Quantum, 2020, 4, 235.
43. H. Kurogi, K. Endo, Y. Sato, M. Sugawara, K. Wada, K. Sugisaki, S. Kanno, H. C. Watanabe and H. Nakano, Optimizing a parameterized controlled gate with Free Quaternion Selection, arXiv, 2024, preprint, arXiv:2409.13547, DOI: 10.48550/arXiv.2409.13547.
44. A. D. Becke, J. Chem. Phys., 1993, 98, 1372–1377.
45. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
46. R. Ditchfield, W. J. Hehre and J. A. Pople, J. Chem. Phys., 1971, 54, 724–728.
47. W. J. Hehre, R. Ditchfield and J. A. Pople, J. Chem. Phys., 1972, 56, 2257–2261.
48. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Revision C.01, Gaussian Inc., Wallingford CT, 2016.
49. G. M. J. Barca, C. Bertoni, L. Carrington, D. Datta, N. De Silva, J. Emiliano Deustua, D. G. Fedorov, J. R. Gour, A. O. Gunina, E. Guidez, T. Harville, S. Irle, J. Ivanic, K. Kowalski, S. S. Leang, H. Li, W. Li, J. J. Lutz, I. Magoulas, J. Mato, V. Mironov, H. Nakata, B. Q. Pham, P. Piecuch, D. Poole, S. R. Pruitt, A. P. Rendell, L. B. Roskop, K. Ruedenberg, T. Sattasathuchana, M. W. Schmidt, J. Shen, L. Slipchenko, M. Sosonkina, V. Sundriyal, A. Tiwari, J. L. Galvez Vallejo, B. Westheimer, M. Wloch, P. Xu, F. Zahariev and M. S. Gordon, J. Chem. Phys., 2020, 152, 154102.
50. W. J. Hehre, R. F. Stewart and J. A. Pople, J. Chem. Phys., 1969, 51, 2657–2664.
51. J. Pipek and P. G. Mezey, J. Chem. Phys., 1989, 90, 4916–4926.
52. J. R. McClean, N. C. Rubin, K. J. Sung, I. D. Kivlichan, X. Bonet-Monroig, Y. Cao, C. Dai, E. Schuyler Fried, C. Gidney, B. Gimby, P. Gokhale, T. Häner, T. Hardikar, V. Havlìcek, O. Higgott, C. Huang, J. Izaac, Z. Jiang, X. Liu, S. McArdle, M. Neeley, T. O'Brien, B. O'Gorman, I. Ozfidan, M. D. Radin, J. Romero, N. P. D. Sawaya, B. Senjean, K. Setia, S. Sim, D. S. Steiger, M. Steudtner, Q. Sun, W. Sun, D. Wang, F. Zhang and R. Babbush, Quantum Sci. Technol., 2020, 5, 034014.
53. A. Tranter, P. J. Love, F. Mintert and P. V. Coveney, J. Chem. Theory Comput., 2018, 14, 5617–5630.
54. Quantum AI team and collaborators qsim, 2020 DOI:10.5281/zenodo.4023103.
55. M. Richer, G. Sánchez-Díaz, M. Martínez-González, V. Chuiko, T. D. Kim, A. Tehrani, S. Wang, P. B. Gaikwad, C. E. V. de Moura, C. Masschelein, R. A. Miranda-Quintana, A. Gerolin, F. Heider-Zadeh and P. W. Ayers, J. Chem. Phys., 2024, 161, 132502.
56. A. A. Saki, S. Barison, B. Fuller, J. R. Garrison, J. R. Glick, C. Johnson, A. Mezzacapo, J. Robledo-Moreno, M. Rossmannek, P. Schweigert, I. Sitdikov and K. J. Sung, Qiskit addon: sample-based quantum diagonalization, 2024, https://github.com/Qiskit/qiskit-addon-sqd.
57. Q. Sun, X. Zhan, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, Z. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov and G. K.-L. Chan, J. Chem. Phys., 2020, 153, 024109.
58. W. Zeng and J. Wu, Chem, 2021, 7, 358–386.
59. A. J. Baldacchino, M. I. Collins, M. P. Nielsen, T. W. Schmidt, D. R. McCamey and M. J. Y. Tayebjee, Chem. Phys. Rev., 2022, 3, 021304.
60. M. Mikkelsen and Y. O. Nakagawa, Quantum-selected configuration interaction with time-evolved state, arXiv, 2024, preprint, arXiv:2412.13839 DOI:10.48550/arXiv.2412.13839.
61. J. Yu, J. Robledo Moreno, J. Iosue, L. Bertels, D. Claudino, B. Fuller, P. Groszkowski, T. S. Humble, P. Jurcevic, W. Kirby, T. A. Maier, M. Motta, B. Pokharel, A. Seif, A. Shehata, K. J. Sung, M. C. Tran, V. Tripathi, A. Mezzacapo and K. Sharma, Quantum-centric algorithm for sample-based Krylov diagonalization, arXiv, preprint, arXiv:2501.09702, DOI: 10.48550/arXiv.2501.09702.
62. B. Singh and J. S. Brinen, J. Am. Chem. Soc., 1971, 93, 540–542.
63. A. Nicolaides, H. Tomioka and S. Murata, J. Am. Chem. Soc., 1998, 120, 11530–11531.
64. S. Nimura, O. Kikuchi, T. Ohana, A. Yabe and M. Kaise, Chem. Lett., 2016, 25, 125–126.
65. A. S. Ichimura, K. Sato, T. Kinoshita, T. Takui, K. Itoh and P. M. Lahti, Mol. Cryst. Liq. Cryst., 1995, 271/272, 279–288.
66. K. Sugisaki, K. Toyota, K. Sato, D. Shiomi, M. Kitagawa and T. Takui, EPR of Free Radicals in Solids I: Trends in Methods and Applications, Springer, Dordrecht, 2nd edn, 2013, pp. 363–392.
67. K. Sugisaki, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi and T. Takui, ACS Cent. Sci., 2019, 5, 167–175.
68. M. Liu, V. I. Artyukhov, H. Lee, F. Xu and B. I. Yakobson, ACS Nano, 2013, 7, 10075–10082.
69. M. Liu, V. I. Artyukhov, H. Lee, F. Xu and B. I. Yakobson, ACS Nano, 2017, 11, 5186.
70. H. Gao, S. Imamura, A. Kasagi and E. Yoshida, J. Chem. Theory Comput., 2024, 20, 1185–1192.
71. P. Reinholdt, K. M. Ziems, E. R. Kjellgren, S. Coriani, S. P. A. Sauer and J. Kongsted, J. Chem. Theory Comput., 2025, 21, 6811–6822.

---