SqDRIFT algorithm for ground state estimation

Usage estimate: 180 seconds on a Heron r3 processor (NOTE: This is an estimate only. Your runtime might vary.)

Learning outcomes

• Learn how to create smaller depth circuits as compared to Trotterization
• Walk through an end-to-end workflow for ground state estimation using qDRIFT and SQD
• Learn how to use qiskit-fermions in tandem with other Qiskit addons to implement such a workflow

Prerequisites

• Read the Sample-based quantum diagonalization (SQD) overview
• Read the Sample-based Krylov Quantum Diagonalization (SKQD) lesson

Background

SqDRIFT is a variant of SKQD that replaces the need to choose an ansatz from which to sample bitstrings with an ensemble of time-evolution circuits constructed directly from the target Hamiltonian. This is achieved by subsampling smaller time-evolution operators from the Hamiltonian based on its coefficients, which is known as the qDRIFT Trotterization method.

Both implementations use Qiskit Fermions to represent the molecular Hamiltonian and construct circuits for the qDRIFT algorithm. The Python implementation uses fermionic transpiler passes, while the C++ implementation explicitly samples groups and maps them to qubit operators.

Let the Hamiltonian be of the form:

where, without loss of generality, we require $c_i > 0$ and that the largest eigenvalue of $h_i$ be equal, in absolute value, to $1$. Any signed or complex prefactor is absorbed into $h_i$, so the coefficients $c_i$ are strictly positive weights while the $h_i$ carry the direction of each term. Here $N$ is the number of terms (or, after grouping, the number of groups) in the Hamiltonian; it is a property of the Hamiltonian and is distinct from the number of operators sampled into a single circuit, written $n$ below.

The qDRIFT algorithm then realizes, for the target time $t$, some operator $V_k$, where $k$ goes from $1 \cdots K$ and signifies the $k_{th}$ SqDRIFT circuit, defined as:

Here $n$ is the number of sampled operators per circuit and $K$ is the number of circuits in the ensemble. The product runs over the $n$ draws, not over all $N$ Hamiltonian terms, and because the terms are drawn with replacement, the same $h_i$ can appear more than once in a single $V_k$.

The quantity:

is the $L_1$ norm of the coefficients, so each of the $n$ steps evolves for the same duration $\lambda t / n$ regardless of which term was drawn. The uniformity of the step angle is the characteristic feature of qDRIFT: a coefficient influences the result through how often its term is drawn, not through how far that term is rotated. The indices are sampled from the distribution:

so the series $(k_1, \ldots, k_n)$ is a random sequence of term indices drawn from this distribution. Since the $c_i$ are positive and sum to $\lambda$, this is a normalized probability distribution, and the expectation of the resulting channel over the random draws approximates evolution under $H$, with an error that decreases as $n$ grows. Note that the approximation error depends on $\lambda$ rather than on the number of terms $N$.

(The SqDRIFT paper writes the number of terms as $\mathcal{N}$ and the sequence length as $N$; we use $N$ and $n$ here to keep the two clearly distinct.)

This tutorial shows how to generate an ensemble of such randomized circuits. After we have created these circuits, similar to how we create a Krylov subspace for different operators, we sample bitstrings from multiple such operators with different time parameters. This ensures a higher overlap between the ground state vectors and sampled bitstrings.

Requirements

Before starting this tutorial, make sure you have installed

• A Python (>=3.10) virtual environment
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (Note that the name is plural)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

You can install all required packages with:

pip install "qiskit~=2.5" "qiskit-fermions==0.1.0" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy

Setup

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

# Qiskit core
from qiskit import transpile
from qiskit.primitives import BitArray

# Qiskit Aer
from qiskit_aer import AerSimulator

# IBM Quantum Compute Service
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler

# Qiskit Fermions
from qiskit_fermions.operators.library import FCIDump
from qiskit_fermions.operators import FermionOperator
from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms
from qiskit_fermions.operators.terms.grouping import (
    group_terms_by_electronic_structure,
)
from qiskit_fermions.operators.terms.ordering import canonical_order
from qiskit_fermions.circuit import FermionicCircuit
from qiskit_fermions.circuit.library import Evolution
from qiskit_fermions.transpiler import FermionicPassManager
from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
from qiskit_fermions.transpiler.passes import QDriftTrotterization
from qiskit_fermions.circuit.library import InitializeModes

# Qiskit addon SQD
from qiskit_addon_sqd.fermion import (
    diagonalize_fermionic_hamiltonian,
    SCIResult,
)

Implementation

Step 1: Map classical inputs to a quantum problem

The Python steps below use the Aer simulator. The hardware example follows afterward.

Reading and Preparing the FCIDump

For this tutorial, we will load up the electronic structure Hamiltonian for nitrogen (N2). There are other ways to create fermionic operators as well. Refer to the documentation at qiskit_fermions.operators.library.

About this FCIDump. The file N2_sto_3g describes a nitrogen molecule ($N_2$) in the minimal STO-3G basis at an interatomic separation of 1.09 $\AA$, the experimental equilibrium bond length. Its header declares NORB=10, NELEC=14, and MS2=0: 10 spatial orbitals (hence 20 spin orbitals, and 20 qubits under Jordan-Wigner), 14 electrons in a spin singlet, so seven $\alpha$ and seven $\beta$ electrons. All orbitals are given symmetry label 1, that is, no point-group symmetry is exploited. Being a full-space STO-3G dump, no orbitals are frozen and the correlation space is small enough that an exact FCI reference energy can be computed classically for comparison, as shown in the next code block.

An equivalent file can be regenerated with PySCF:

from pyscf import gto, scf, tools

mol = gto.M(atom="N 0 0 0; N 0 0 1.09", basis="sto-3g", symmetry=False)
mf = scf.RHF(mol).run()
tools.fcidump.from_scf(mf, "N2_sto_3g")

Because the integrals depend on the converged SCF orbitals, a regenerated file might differ from the shipped one in orbital phase or ordering; the total energies are unaffected.

Obtaining the file. Find the FCIDump in this GitHub repository. You can run the code below to fetch it into the location the rest of the tutorial expects.

First we use the cisolver provided by pyscf to get the reference energy. This is the true ground state energy of the molecule we are working with. For this we will first declare norb and nelec, which are the number of orbitals and the number of electrons, respectively. Then we declare h1e and h2e, which are the one- and two-electron integrals respectively. All of these will later be used for SQD as well.

import os
from urllib.request import urlopen

# The FCIDump is stored with this tutorial in the Qiskit documentation repository.
FCIDUMP_URL = "https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g"
FCIDUMP_PATH = "fcidump_files/N2_sto_3g"

if not os.path.exists(FCIDUMP_PATH):
    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)
    with urlopen(FCIDUMP_URL) as response:
        contents = response.read()
    with open(FCIDUMP_PATH, "wb") as f:
        f.write(contents)
    print(f"Downloaded FCIDump to {FCIDUMP_PATH}")
else:
    print(f"Using existing FCIDump at {FCIDUMP_PATH}")

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

Loading the Hamiltonian

With the necessary data ready, we read the Hamiltonian from the FCI file in a format that is compatible with qiskit-fermions

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

Fermionic workflows with qiskit-fermions

We will first map the Hamiltonian into a fermionic circuit model using qiskit-fermions, which provides transpiler passes and gates specific to fermionic circuits. These will later be used before Qiskit's traditional transpiler passes for this workflow.

Term grouping

To ensure the reproducibility of results, we first use canonical_order to sort the terms based only on their structure. The order of the operators in the canon list is therefore fixed. This ensures reproducibility of created operators because the QDriftTrotterization pass that we will use in the future samples random indices to create the qDRIFT operators.

In this step, we exploit the many symmetries that are present in the electronic structure Hamiltonian by grouping related terms with identical coefficients. While doing so changes the operator coefficient distribution which the qDRIFT protocol samples from, this does not affect its convergence guarantees. Crucially, grouping terms related by symmetry results in a favorable cancellation of Pauli terms and in an overall shorter circuit depth when time-evolving a state under their action.

qiskit-fermions provides the group_terms_by_electronic_structure function that does this grouping for us.

Note that the group_terms_by_electronic_structure assumes normal ordering terms.

Filtering diagonal terms

We remove the diagonal terms from the Hamiltonian used to generate the circuits, so that the $n$ qDRIFT sampling slots are spent on terms that move population between configurations. Such terms are best filtered out of the Hamiltonian at this point, before the Evolution gate is constructed in the next step.

The terms in question are the ones that are diagonal in the occupation-number basis, that is, the products of number operators $a^\dagger_i a_i$. Three kinds of term fall under this description:

• the constant energy offset, a product of zero number operators, whose time evolution contributes only a global phase;
• the individual number operators $n_i$, whose time evolution reduces to single-qubit $Z$ rotations;
• the higher-order products such as $n_i n_j$.

On their own, none of these move population between occupation-number configurations; they act only on the phases of the configurations already present. They are not inert, however: those relative phases feed into the interference generated by the excitation terms later in the circuit, so filtering them changes the evolution that is actually generated and can change the sampling distribution. This is a deliberate approximation in the circuit-generation step, made to focus sampling on excitation terms, rather than a step that leaves the sampled distribution untouched. Unlike the symmetry grouping above, which leaves the qDRIFT convergence guarantees intact, this filter changes the operator being evolved. The circuits therefore no longer approximate evolution under the full Hamiltonian, and the qDRIFT error bounds apply to the filtered operator rather than the original one. This is acceptable here because the circuits are only a sampling heuristic used to propose configurations: no term is lost from the energy estimate itself, since the filter applies only to the Hamiltonian used to build the circuits, while the classical diagonalization later uses the full Hamiltonian, diagonal terms included. SQD's accuracy depends on that classical step, which remains variational in the sampled subspace regardless of how the configurations were proposed.

The filter_diagonal_terms() function removes such terms from an operator in place. It identifies them from their normal-ordered structure — the multiset of creation modes matching the multiset of annihilation modes — so it is only valid on an operator that is already normal-ordered. This assumption is not checked at runtime.

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

5060

Now that we have grouped the terms in the Hamiltonian, we will decide on the following parameters to generate the ensemble of circuits:

• The number of circuits to generate: num_circuits
• The length of each circuit in terms of excitation groups: num_exc
• The factor for the different evolution times: times

Creating fermionic circuits

We will now create fermionic circuits for each of the time-steps. Each circuit will consist of a single evolution gate, with the evolution time we declared earlier. The evolution operator is the Hamiltonian. Later we run transpiler passes on these circuits to create qDRIFT circuits.

Ansatz preparation

We prepare the Hartree-Fock state using the InitializeModes class. For nitrogen, the process is simply applying X gates to the first num_elec_a qubits and then to the num_elec_b qubits, both of which equal to seven for nitrogen. This state represents the seven $\alpha$ and seven $\beta$ electrons of nitrogen.

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate


init_circuits = []

hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

Step 2: Optimize problem for quantum hardware execution

Now that we have our circuits, we will first use the passes available in qiskit-fermions to perform fermionic level optimizations, followed by transpiling our circuit for the backend of choice. Since this is a simulator experiment, we will first do this for the AerSimulator.

Weight calculation for each group

In this step, we perform the qDRIFT sampling of terms stochastically with probabilities proportional to their coefficients in the Hamiltonian. The qDRIFT transpiler pass does this for us. We can now create shallower circuits that can be executed on the hardware more efficiently despite limited qubit connectivity, even when the Hamiltonian contains long-range couplings and higher-than-quadratic terms. After term grouping, it samples the operators based on their weights. For each operator $h_i$, the weight $W_{h_i}$ is defined as follows:

Because the terms were grouped in Step 1, each $h_i$ here is a whole group: $c_i$ is the mean absolute coefficient of the terms in group $i$, and each term in the group is evolved with its coefficient reduced to its sign.

Fermionic and hardware-native optimizations

The function generate_preset_jw_pass_manager() returns a MultiStagePassManager that takes a FermionicCircuit and produces an optimized final circuit that we can transpile to run on our hardware. We replace its default optimization stage with a FermionicPassManager containing our QDriftTrotterization pass:

• The QDriftTrotterization pass uses the weight-calculation and sampling internally to generate the circuits that we will use for sampling
• The RelabelModes pass is another optimization pass that can be used to permute the fermionic modes to optimize connectivity across qubits and reduce gate depth; read more in the API reference

The remaining stages of the MultiStagePassManager run automatically and handle the full fermion-to-qubit mapping:

• F2QLayout: The preset pass manager applies the TrivialF2QLayout pass, which trivially maps $n$ fermionic bits to $n$ qubits.
• F2QSynth: A transpilation pass to map fermion-based circuit instructions to qubit-based ones.

qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

400

Now that we are done with the fermionic-level optimizations, we can transpile the circuits for execution on the simulator.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Step 3: Execute using Qiskit primitives

Now that we have our circuits, we can run them using Qiskit primitives on the AerSimulator. We will combine all the counts from different circuits. We convert them to boolean vectors before finally post-processing with SQD.

print(
    f"Executing {len(transpiled_circuits)} circuits with {shots} shots each..."
)

job = simulator.run(transpiled_circuits, shots=shots)
result = job.result()

all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]

print(len(all_counts), "length before post processing")

Executing 400 circuits with 100 shots each...
400 length before post processing

Step 4: Post-process and return result in desired classical format

Using bitstrings for SQD

We can now run the diagonalization scheme on the selected bitstrings to find the lowest eigenvalue that will correspond to the ground state energy of the molecule. We create a callback function, declare initial occupancies, and set the parameters before finally running the diagonalization scheme. The callback function is used to print the current iteration and the current eigenvalue estimate at each iteration.

Finally, to get the ground state estimate, we add the nuclear_repulsion_energy to the resultant energy.

Note: The subspace dimension is not fixed across iterations, even on the noiseless simulator — each subsample draws a different set of configurations, and the recovery step reshapes the pool between iterations, so the reported dimension varies from one subsample to the next. Noiseless sampling does not by itself pin the selected-subspace dimension. The hardware run, however, tends to give systematically larger subspaces, because noisy shots break particle-number symmetry and configuration recovery turns them into additional basis vectors. Because of that, we will also introduce another step for pruning bitstrings in the hardware section.

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

result_history = []


def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )


# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")


energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

40000
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64767025226178
		Subspace dimension: 5538
	Subsample 1
		Energy: -107.64772799119115
		Subspace dimension: 5670
	Subsample 2
		Energy: -107.64765512281548
		Subspace dimension: 5767
Iteration 2
	Subsample 0
		Energy: -107.64795948524682
		Subspace dimension: 6080
	Subsample 1
		Energy: -107.64806617355072
		Subspace dimension: 6300
	Subsample 2
		Energy: -107.64802260640258
		Subspace dimension: 6308
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712
 0.99301927 0.0373282  0.0373266  0.00944508]
Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737
 0.99302145 0.03733536 0.03733399 0.0094585 ]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6480661736 Ha
Error:            1.1811817564e-04 Ha

Hardware example

This example uses 20 qubits (10 spatial orbitals). That choice is a convenience for a tutorial that should run quickly, not a hard ceiling on the method.

The cost of the classical step is not set by the qubit count directly. SQD diagonalizes the Hamiltonian projected onto the subspace spanned by the sampled configurations, so what drives the classical cost is the dimension of that selected subspace — governed here by samples_per_batch, num_batches, and how many distinct configurations the circuits actually produce — together with the sparse linear algebra needed to apply the projected Hamiltonian. The full CI space grows combinatorially with orbitals and electrons, but the selected subspace is a small, tunable slice of it, and we control its size directly. Consequently, the number of qubits and the classical difficulty can be varied somewhat independently: a wider orbital space sampled into a modest subspace can be cheaper than a smaller system diagonalized over a very large one.

In practice, then, the feasible system size depends on the subspace dimension you need for the accuracy you want and on the memory and cores available to the eigensolver. Larger orbital spaces typically do require a larger subspace to reach chemical accuracy, and that is what eventually motivates distributed resources — see qiskit-addon-sqd-hpc for scaling this step out. Rather than assuming a fixed cutoff, the practical approach is to watch the reported subspace dimension and the energy convergence across iterations and increase the subspace size until the energy stops improving or you exhaust available memory.

Note: Due to sampling error from the noise in the hardware, the subspace created for diagonalization in the hardware run will be larger than what we get when using the simulator. While it increases the dimension of the subspace we want to diagonalize, the workflow still gives us an accurate answer due to the robustness of SQD towards noise.

Pruning of spurious strings

Here we can choose to perform an additional step. When we have all the bitstrings from the circuit executions, we can either filter out the invalid bitstrings before running SQD, or move forward without pruning. Skipping the pruning is generally preferable for hardware runs, because it leaves the symmetry-broken shots available to configuration recovery, which can repair them into valid configurations and thereby widen the subspace instead of discarding those shots outright.

Since nitrogen can only have seven $\alpha$ and seven $\beta$ electrons, any bitstrings that have more or fewer than seven 1s in the first and the second half of the output can be discarded. We define a function that checks if the bitstrings are valid, and if not, discards them. Once we filter out the spurious bitstrings, the rest are sent into the diagonalization scheme. Use the PRUNE flag below to switch between the two behaviors.

Keep in mind that pruning is only one of several choices that shape the final subspace, alongside the number of circuits, the set of evolution times, and diagonal-term filtering. Comparing a pruned run against an unpruned one is only informative if everything else is held fixed; the C++ tab discusses this in more detail, since it postselects rather than recovers and also differs in those other parameters.

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []
hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

# Calculate weights for sampling (one per group)
qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService(channel="ibm_quantum_platform")

# Select backend (choose based on qubit requirements)
backend = service.least_busy(
    operational=True,
    simulator=False,
    min_num_qubits=2 * norb,
)

print(f"Selected backend: {backend.name} ({backend.num_qubits} qubits)")

# Transpile for hardware
transpiled_circuits = transpile(
    sqdrift_circuits,
    backend=backend,
    optimization_level=3,
    seed_transpiler=42,
)

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-SqDRIFT"]

job = sampler.run(transpiled_circuits, shots=shots)
result = job.result()

# Extract counts from SamplerV2 results
all_counts = [pub_result.data.meas.get_counts() for pub_result in result]

# Set to True to filter out bitstrings that violate electron-number conservation
PRUNE = False


def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )


if PRUNE:
    all_counts_filtered = []
    for counts in all_counts:
        filtered_count = {}
        for key in counts:
            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):
                continue
            elif key not in filtered_count.keys():
                filtered_count[key] = counts[key]
            else:
                filtered_count[key] += counts[key]
        all_counts_filtered.append(filtered_count)
    all_counts = all_counts_filtered

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print("Electron configuration:")
print(f"  Total electrons: {nelec}")
print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

result_history = []


def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )


# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")


energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha
5060
400

Selected backend: ibm_aachen (156 qubits)
40000
Electron configuration:
  Total electrons: 14
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64593072647523
		Subspace dimension: 7221
	Subsample 1
		Energy: -107.6458270048177
		Subspace dimension: 7209
	Subsample 2
		Energy: -107.64007673117075
		Subspace dimension: 7138
Iteration 2
	Subsample 0
		Energy: -107.64757372124944
		Subspace dimension: 9009
	Subsample 1
		Energy: -107.64674060104392
		Subspace dimension: 8245
	Subsample 2
		Energy: -107.64731360491942
		Subspace dimension: 8178
Iteration 3
	Subsample 0
		Energy: -107.64765518770588
		Subspace dimension: 8835
	Subsample 1
		Energy: -107.64767975712016
		Subspace dimension: 8649
	Subsample 2
		Energy: -107.64761634415606
		Subspace dimension: 8648
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295
 0.99298797 0.03728154 0.03728186 0.00938359]
Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174
 0.99298424 0.03728121 0.03728169 0.00939159]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6476797571 Ha
Error:            5.0453460619e-04 Ha

Usage estimate: 180 seconds on a Heron r3 processor (NOTE: This is an estimate only. Your runtime might vary.)

Learning outcomes (C++)

• Learn how to create smaller depth circuits as compared to Trotterization
• Walk through an end-to-end workflow for ground state estimation using qDRIFT and SQD
• Learn how to use qiskit-fermions in tandem with other Qiskit addons to implement such a workflow

Prerequisites (C++)

• Read the Sample-based quantum diagonalization (SQD) overview
• Read the Sample-based Krylov Quantum Diagonalization (SKQD) lesson

Background (C++)

SqDRIFT is a variant of SKQD that replaces the need to choose an ansatz from which to sample bitstrings with an ensemble of time-evolution circuits constructed directly from the target Hamiltonian. This is achieved by subsampling smaller time-evolution operators from the Hamiltonian based on its coefficients, which is known as the qDRIFT Trotterization method.

Both implementations use Qiskit Fermions to represent the molecular Hamiltonian and construct circuits for the qDRIFT algorithm. The Python implementation uses fermionic transpiler passes, while the C++ implementation explicitly samples groups and maps them to qubit operators.

Let the Hamiltonian be of the form:

where, without loss of generality, we require $c_i > 0$ and that the largest eigenvalue of $h_i$ be equal, in absolute value, to $1$. Any signed or complex prefactor is absorbed into $h_i$, so the coefficients $c_i$ are strictly positive weights while the $h_i$ carry the direction of each term. Here $N$ is the number of terms (or, after grouping, the number of groups) in the Hamiltonian; it is a property of the Hamiltonian and is distinct from the number of operators sampled into a single circuit, written $n$ below.

The qDRIFT algorithm then realizes, for the target time $t$, some operator $V_k$, where $k$ goes from $1 \cdots K$ and signifies the $k_{th}$ SqDRIFT circuit, defined as:

Here $n$ is the number of sampled operators per circuit and $K$ is the number of circuits in the ensemble. The product runs over the $n$ draws, not over all $N$ Hamiltonian terms, and because the terms are drawn with replacement, the same $h_i$ can appear more than once in a single $V_k$.

The quantity:

is the $L_1$ norm of the coefficients, so each of the $n$ steps evolves for the same duration $\lambda t / n$ regardless of which term was drawn. The uniformity of the step angle is the characteristic feature of qDRIFT: a coefficient influences the result through how often its term is drawn, not through how far that term is rotated. The indices are sampled from the distribution:

so the series $(k_1, \ldots, k_n)$ is a random sequence of term indices drawn from this distribution. Since the $c_i$ are positive and sum to $\lambda$, this is a normalized probability distribution, and the expectation of the resulting channel over the random draws approximates evolution under $H$, with an error that decreases as $n$ grows. Note that the approximation error depends on $\lambda$ rather than on the number of terms $N$.

(The SqDRIFT paper writes the number of terms as $\mathcal{N}$ and the sequence length as $N$; we use $N$ and $n$ here to keep the two clearly distinct.)

This tutorial shows how to generate an ensemble of such randomized circuits. After we have created these circuits, similar to how we create a Krylov subspace for different operators, we sample bitstrings from multiple such operators with different time parameters. This ensures a higher overlap between the ground state vectors and sampled bitstrings.

Requirements (C++)

This implementation runs on IBM Quantum® hardware and uses an external SBD executable for classical diagonalization. It does not include the Python simulator example.

Follow the C++ build instructions to install the compiler, CMake, Rust, Boost, and the Qiskit dependencies. Install Rust with rustup: the Qiskit dependencies pin their own Rust version, which rustup downloads automatically, but a Rust installed from Homebrew or a Linux package manager might be too old. The final SBD step additionally requires MPI, OpenMP, and BLAS/LAPACK. Windows is not currently supported by this implementation.

Use the complete C++ source and CMake project from the repository. Keep the adjacent fcidump_files directory when downloading the example. From a checkout of the documentation repository, build with the following commands. The first build downloads and compiles all dependencies and takes several minutes.

cd docs/tutorials/assets/sqdrift/cpp
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --config Release

The build writes the SqDRIFT executable to docs/tutorials/assets/sqdrift/cpp (the project root), not to build/. Run ./SqDRIFT from the project root so that it can find ../fcidump_files/N2_sto_3g. All other commands in this tab also run from the project root unless they say otherwise. If ./SqDRIFT fails with a missing shared library error, follow the platform-specific library-path instructions in the build guide; on macOS, the executable usually runs without them.

Setup (C++)

The snippets below are excerpts from the complete SqDRIFT.cpp program. Build and run that file as a single program.

Before running SqDRIFT, ensure your IBM Quantum API token is configured for IBM Cloud Platform.

If you have Python with qiskit-ibm-runtime installed, the easiest way to create a correctly formatted config file is:

from qiskit_ibm_runtime import QiskitRuntimeService
QiskitRuntimeService.save_account(channel="ibm_quantum_platform", token="YOUR_IBM_QUANTUM_TOKEN_HERE")

This writes ~/.qiskit/qiskit-ibm.json (Linux/macOS) or %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) with the right shape. If you prefer to create the file manually:

{
  "default-ibm-quantum-platform": {
    "channel": "ibm_quantum_platform",
    "token": "YOUR_IBM_QUANTUM_TOKEN_HERE"
  }
}

Get your API key from IBM Quantum Platform.

Note: A missing or malformed config file causes an unhelpful Rust unwrap() panic at startup rather than a clean error. Double-check the key name (default-ibm-quantum-platform) and channel value (ibm_quantum_platform) if you encounter a problem.

Implementation (C++)

Step 1: Map classical inputs to a quantum problem (C++)

Load the molecular Hamiltonian

char filename[] = "../fcidump_files/N2_sto_3g";
QfFCIDump* fcidump = qf_fcidump_from_file(filename);
if (fcidump == nullptr) {
    std::cerr << "Failed to load FCIDump!" << std::endl;
    return 1;
}

uint32_t norb = qf_fcidump_norb(fcidump);
uint32_t nelec = qf_fcidump_nelec(fcidump);
uint32_t n_alpha = nelec / 2;
uint32_t n_beta  = nelec - n_alpha;
uint32_t num_modes = 2 * norb;
uint32_t num_qubits = num_modes;
std::cout << "✓ Loaded N2 molecule (" << norb << " orbitals, "
          << num_qubits << " qubits, " << nelec << " electrons)" << std::endl;

Output:

✓ Loaded N2 molecule (10 orbitals, 20 qubits, 14 electrons)

Loads the $N_2$ molecule Hamiltonian from an FCIDump file. The number of qubits equals twice the number of orbitals (spin-up and spin-down).

Create and normal-order the Hamiltonian

// Create and normal order Hamiltonian
QfFermionOperator* hamiltonian = qf_ferm_op_from_fcidump(fcidump);
QfFermionOperator* normal_ordered = qf_ferm_op_normal_ordered(hamiltonian, nullptr);

Converts the FCIDump data into a fermionic operator and applies normal ordering to simplify the operator structure.

Group terms by electronic structure

// Group terms by electronic structure
QfExitCode exit_code = qf_ferm_op_group_terms_by_electronic_structure(
    normal_ordered, num_modes, false
);
if (exit_code != QfExitCode_Success) {
    std::cerr << "Failed to group terms!" << std::endl;
    return 1;
}

uint32_t num_groups = qf_ferm_op_num_groups(normal_ordered);
std::cout << "Grouped into " << num_groups << " groups" << std::endl;

// Split into group operators
QfFermionOperator** group_ops = new QfFermionOperator*[num_groups];
qf_ferm_op_split_out_groups(normal_ordered, nullptr, 0, group_ops);

Output:

Grouped into 1902 groups

Groups related Hamiltonian terms by electronic structure before constructing their evolution circuits.

Calculate sampling weights, then normalize and map to qubit operators

Order matters here. Normalization rescales every coefficient to unit magnitude in place, so the sampling weights must be computed first. If normalization runs first, every $|c_j|$ reads back as exactly 1.0 and each group's weight collapses into a bare count of its terms — the sampling distribution would then be driven by group size rather than by the physical coefficients.

Each group's weight is the mean absolute coefficient of its terms. The Python tab's QDriftTrotterization pass uses the same weight for a grouped Hamiltonian. Using the sum instead would make $\lambda$ scale with group size, which in turn skews the evolution time $t_k = \lambda\tau/n$.

// Calculate sampling weights BEFORE normalization.
std::vector<double> weights(num_groups);
double total_weight = 0.0;

for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    double abs_sum = 0.0;
    for (uint64_t j = 0; j < num_terms; j++) {
        abs_sum += std::sqrt(coeffs[j].re * coeffs[j].re +
                             coeffs[j].im * coeffs[j].im);
    }
    weights[i] = (num_terms > 0) ? (abs_sum / static_cast<double>(num_terms)) : 0.0;
    total_weight += weights[i];
}

std::cout << "Total weight (λ): " << total_weight << std::endl;

// Normalize the evolution operators term-by-term, now that the sampling
// weights above have already captured the original magnitudes.
for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    for (uint64_t j = 0; j < num_terms; j++) {
        const double magnitude = std::sqrt(coeffs[j].re * coeffs[j].re +
                                           coeffs[j].im * coeffs[j].im);
        if (magnitude > 0.0) {
            coeffs[j].re /= magnitude;
            coeffs[j].im /= magnitude;
        }
    }
}

QkObs** qubit_ops = new QkObs*[num_groups];
for (uint32_t i = 0; i < num_groups; i++) {
    QfExitCode jw_exit = qf_ferm_op_jordan_wigner(group_ops[i], num_qubits, &qubit_ops[i]);
    if (jw_exit != QfExitCode_Success) {
        std::cerr << "Failed to map group " << i << std::endl;
        return 1;
    }
}
std::cout << "Mapped all " << num_groups << " normalized groups to qubit operators" << std::endl;

Output:

Total weight (λ): 337.261
Mapped all 1902 normalized groups to qubit operators

The weight $\lambda_i = \frac{1}{|G_i|}\sum_{j \in G_i} |c_{ij}|$ is the mean absolute coefficient over the terms of group $G_i$, read from the original coefficients. Normalization then rescales each fermionic term to unit magnitude before the Jordan-Wigner transformation, preserving only its phase/sign in the evolved operator. This matches qDRIFT's normalized-term evolution: coefficient magnitudes determine sampling probabilities via $\lambda$, while the circuit evolution uses the normalized term so large coefficients are not counted twice.

Sample operator sets

// SqDRIFT Sampling: Create operator sets
    const int num_circuits = 100;    // Number of circuits to create
    const int ops_per_circuit = 10;  // Operators per circuit
    const double time_step = 1;    // Time step for evolution
    
    std::cout << "\n Generating " << num_circuits << " circuits with " 
              << ops_per_circuit << " operators each..." << std::endl;
    
    std::mt19937 gen(42);
    std::discrete_distribution<> dist(weights.begin(), weights.end());
    
    // Store sampled operator indices
    std::vector<std::vector<int>> operator_sets(num_circuits);
    for (int i = 0; i < num_circuits; i++) {
        operator_sets[i].resize(ops_per_circuit);
        for (int j = 0; j < ops_per_circuit; j++) {
            operator_sets[i][j] = dist(gen);
        }
    }

Output:

 Generating 100 circuits with 10 operators each...

Draws the operator sets for 100 random circuits, each containing 10 groups sampled from the weighted distribution. The next step turns these sets into the ensemble of qDRIFT circuits.

Step 2: Build the qDRIFT circuits (C++)

In the C++ version, this step builds the circuits from the sampled operator sets. Transpilation for the target backend happens in Step 3, immediately before each circuit is submitted.

    // Create Suzuki-Trotter circuits by composing all operators
    std::vector<Qiskit::circuit::QuantumCircuit> circuits;
    circuits.reserve(num_circuits);

    const uint32_t trotter_order = 1;
    const uint32_t trotter_reps = 1;
    const bool preserve_order = false;
    const bool insert_barriers = false;

    // Gate name mapping is static — compute once before all loops
    auto name_map = Qiskit::circuit::get_standard_gate_name_mapping();

    for (int circuit_idx = 0; circuit_idx < num_circuits; circuit_idx++) {
        // Create initial circuit with state preparation using C++ API
        Qiskit::circuit::QuantumCircuit qc(num_qubits, num_qubits);

        // State preparation: apply X gates to spin-up qubits 0..n_alpha-1
        // and spin-down qubits norb..norb+n_beta-1 (derived from FCIDump NELEC)
        for (uint32_t q = 0; q < n_alpha && q < num_qubits; q++) {
            qc.x(q);
        }
        for (uint32_t q = norb; q < norb + n_beta && q < num_qubits; q++) {
            qc.x(q);
        }

        // For each operator in this circuit's set
        for (int op_idx = 0; op_idx < ops_per_circuit; op_idx++) {
            int group_idx = operator_sets[circuit_idx][op_idx];
            double evolution_time = (total_weight * time_step) / ops_per_circuit;
            
            // Create Suzuki-Trotter evolution circuit
            QkCircuit* evolution_circuit = qk_circuit_library_suzuki_trotter(
                qubit_ops[group_idx], trotter_order, trotter_reps,
                evolution_time, preserve_order, insert_barriers
            );
            
            if (evolution_circuit != nullptr) {
                // Get number of instructions
                uint32_t num_ops = qk_circuit_num_instructions(evolution_circuit);
                
                // Manually append each instruction
                for (uint32_t i = 0; i < num_ops; i++) {
                    QkCircuitInstruction inst;
                    qk_circuit_get_instruction(evolution_circuit, i, &inst);

                    // Prepare qubit vector
                    std::vector<uint32_t> qubits(inst.num_qubits);
                    for (uint32_t j = 0; j < inst.num_qubits; j++) {
                        qubits[j] = inst.qubits[j];
                    }

                    // Get operation kind
                    QkOperationKind kind = qk_circuit_instruction_kind(evolution_circuit, i);

                    // Get mutable pointer to circuit
                    QkCircuit* mutable_circuit = qc.get_rust_circuit().get();

                    // Append based on operation type
                    if (kind == QkOperationKind_Gate) {
                        // Convert gate name string to QkGate enum
                        std::string gate_name(inst.name);
                        auto gate_it = name_map.find(gate_name);
                        if (gate_it != name_map.end()) {
                            QkGate gate_enum = gate_it->second.gate_map();
                            qk_circuit_parameterized_gate(
                                mutable_circuit,
                                gate_enum,           // Use QkGate enum, not string
                                qubits.data(),
                                inst.params
                            );
                        } else {
                            std::cerr << "Unknown gate in evolution circuit: " << gate_name << std::endl;
                            qk_circuit_instruction_clear(&inst);
                            qk_circuit_free(evolution_circuit);
                            return 1;
                        }
                    } else if (kind == QkOperationKind_Barrier) {
                        qk_circuit_barrier(
                            mutable_circuit,
                            qubits.data(), 
                            inst.num_qubits
                        );
                    } else if (kind == QkOperationKind_Reset) {
                        qk_circuit_reset(
                            mutable_circuit,
                            qubits[0]
                        );
                    }

                    qk_circuit_instruction_clear(&inst);
                }

                qk_circuit_free(evolution_circuit);
            }
        }

        circuits.push_back(qc);

        // Display progress
        if ((circuit_idx + 1) % 10 == 0) {
            std::cout << "  Created " << (circuit_idx + 1) << "/" 
                    << num_circuits << " circuits" << std::endl;
        }
    }

    std::cout << "✓ Created all " << num_circuits << " Suzuki-Trotter circuits" << std::endl;

Output:

  Created 10/100 circuits
  Created 20/100 circuits
  Created 30/100 circuits
  Created 40/100 circuits
  Created 50/100 circuits
  Created 60/100 circuits
  Created 70/100 circuits
  Created 80/100 circuits
  Created 90/100 circuits
  Created 100/100 circuits
✓ Created all 100 Suzuki-Trotter circuits

The complete program then prints circuit statistics and the operator indices of circuit 0, which are omitted here.

For each sampled operator, creates a Suzuki-Trotter evolution circuit with time $t_k = \frac{\lambda \tau}{n}$, where $n$ is the number of operators sampled per circuit and $\tau$ is the time step. Because each sampled group was normalized term-by-term before mapping, this evolution applies only the phase/sign of each term during the rotation; the magnitudes contribute through the sampling distribution only. The Hartree-Fock initial state is prepared by applying X gates to qubits 0-6 and 10-16 (seven electrons in each spin sector). Each evolution circuit's instructions are manually appended to the main circuit using the C++ API.

Step 3: Execute using Qiskit primitives (C++)

std::cout << "\n Connecting to IBM Quantum Cloud..." << std::endl;

// Initialize IBM Quantum Compute Service using qiskit-cpp
Qiskit::service::QiskitRuntimeService service;
std::cout << "✓ Connected to IBM Quantum" << std::endl;

// Get backend
const std::string backend_name = "ibm_fez";
auto backend = service.backend(backend_name);

if (backend.name().empty()) {
    std::cerr << "Backend " << backend_name
              << " is not available to this account." << std::endl;
    return 1;
}

std::cout << "✓ Selected backend: " << backend.name() << std::endl;

// Create sampler primitive
const int32_t shots = 100;  // Number of shots per circuit
Qiskit::primitives::BackendSamplerV2 sampler(backend, shots);

std::cout << "\n Submitting " << num_circuits << " circuits ("
            << shots << " shots each)..." << std::endl;

// Submit one job per circuit.
std::vector<std::shared_ptr<Qiskit::primitives::BasePrimitiveJob>> jobs;

for (int i = 0; i < num_circuits; i++) {
    for (uint32_t q = 0; q < num_qubits; q++) {
        circuits[i].measure(q, q);
    }
    Qiskit::circuit::QuantumCircuit transpiled_qc =
        Qiskit::compiler::transpile(circuits[i], backend, 2, 1.0, 42);
    std::vector<Qiskit::primitives::SamplerPub> pubs;
    pubs.push_back(Qiskit::primitives::SamplerPub(transpiled_qc, shots));
    jobs.push_back(sampler.run(pubs));
}

std::cout << "\n⏳ Waiting for " << num_circuits << " jobs to complete..." << std::endl;

// Collect results — poll each job and extract bitstrings
std::vector<boost::dynamic_bitset<>> all_bitstrings;

for (int i = 0; i < num_circuits; i++) {
    while (!jobs[i]->in_final_state()) {
        sleep(5);
    }
    if (jobs[i]->done()) {
        auto result = jobs[i]->result();
        if (result.size() > 0) {
            auto& pub_result = result[0];
            auto& bit_array  = pub_result.data();
            for (size_t s = 0; s < bit_array.num_shots(); s++) {
                auto sample = bit_array[s].to_string();
                boost::dynamic_bitset<> bs(num_qubits);
                for (uint32_t bit = 0; bit < num_qubits; bit++) {
                    bs[bit] = (sample[num_qubits - 1 - bit] == '1');
                }
                all_bitstrings.push_back(bs);
            }
        }
    } else {
        std::cout << "  Circuit " << (i+1) << " did not complete successfully" << std::endl;
    }
}

std::cout << "\n✓ Collected " << all_bitstrings.size()
            << " total bitstrings from all circuits" << std::endl;

Connects to IBM Quantum, selects the backend by name, transpiles each circuit for that backend, and submits one job per circuit: 100 separate Sampler jobs of 100 shots each on a 20-qubit register. It then polls each job until it finishes and collects the bitstrings from result[0]. Every job uses QPU time and waits in the queue separately, so expect this step to take much longer than a single batched job.

The hardware results quoted on this page come from a run on ibm_pittsburgh (100 circuits × 100 shots), with backend_name changed from the default ibm_fez. The saved samples were postprocessed with the corrected determinant bit ordering shown in Step 4. Hardware noise varies between runs, so your counts and final energy will differ.

Output:

 Connecting to IBM Quantum Cloud...
✓ Connected to IBM Quantum
✓ Selected backend: ibm_pittsburgh

 Submitting 100 circuits (100 shots each)...

⏳ Waiting for 100 jobs to complete...

✓ Collected 10000 total bitstrings from all circuits

Changing the backend: Change backend_name in SqDRIFT.cpp to target a different device; service.backends() lists what your account can reach.

Step 4: Post-process and return result in desired classical format (C++)

Postselect bitstrings and write CI strings

// Postselect bitstrings by Hamming weight (derived from FCIDump NELEC)
std::vector<double> bitstring_weights(all_bitstrings.size(), 1.0);
auto [filtered_bitstrings, filtered_weights] = Qiskit::addon::sqd::postselect_bitstrings(
    all_bitstrings,
    bitstring_weights,
    Qiskit::addon::sqd::MatchesRightLeftHamming<uint32_t>(n_alpha, n_beta)
);
std::ignore = filtered_weights;

std::cout << "✓ Postselected " << filtered_bitstrings.size()
          << " bitstrings with Hamming weight ("
          << n_alpha << "," << n_beta << ")" << std::endl;

// Convert bitstrings to CI strings using SQD addon
auto ci_strings = Qiskit::addon::sqd::bitstrings_to_ci_strings_symmetrize_spin(
    filtered_bitstrings,
    std::nullopt  // No dimension limit
);

std::cout << "✓ Generated " << ci_strings.size() << " CI strings" << std::endl;

// Write CI strings to file for SBD
std::ofstream alpha_file("alphadets_from_sqd.txt");
for (const auto& ci_string : ci_strings) {
    std::string bitstr;
    // SBD expects the highest orbital index first (bit zero on the right).
    boost::to_string(ci_string, bitstr);
    alpha_file << bitstr << "\n";
}
alpha_file.close();

std::cout << "✓ Wrote " << ci_strings.size()
          << " CI strings to alphadets_from_sqd.txt" << std::endl;
std::cout << " Ready for SBD diagonalization!" << std::endl;

Output:

✓ Postselected 1137 bitstrings with Hamming weight (7,7)
✓ Generated 49 CI strings
✓ Wrote 49 CI strings to alphadets_from_sqd.txt
 Ready for SBD diagonalization!

Postselects bitstrings with the correct Hamming weight (n_alpha spin-up, n_beta spin-down electrons, derived from the FCIDump NELEC field), converts them to Configuration Interaction (CI) strings using spin symmetrization, and writes them to alphadets_from_sqd.txt for subsequent Selected Basis Diagonalization.

Clean up resources

// Cleanup
// circuits vector will be automatically cleaned up (RAII)
// Free qubit operators
for (uint32_t i = 0; i < num_groups; i++) {
    qk_obs_free(qubit_ops[i]);
}
delete[] qubit_ops;
// Free fermionic operators
for (uint32_t i = 0; i < num_groups; i++) {
    qf_ferm_op_free(group_ops[i]);
}
delete[] group_ops;
qf_ferm_op_free(normal_ordered);
qf_ferm_op_free(hamiltonian);
qf_fcidump_free(fcidump);
return 0;

Properly frees all allocated memory and resources. The C++ std::vector<QuantumCircuit> is automatically cleaned up via RAII. Manually frees qubit operators, fermionic operators, and the FCIDump data.

Now that we have created the basis for projecting our Hamiltonian over, we can proceed with the diagonalization process to obtain the ground state estimate. Run the following commands from the project root, docs/tutorials/assets/sqdrift/cpp.

Build the SBD diag binary

CMake fetches SBD as source-only into deps/sbd/ (the main build intentionally skips compiling it). The SBD app uses a plain Makefile, and SBD itself requires MPI, OpenMP, and BLAS/LAPACK — so plain make is not sufficient until those are installed and the app's Configuration file matches your toolchain.

The Makefile reads its compiler and link flags from the Configuration file in the same directory. The version shipped upstream is tuned for one specific macOS/Homebrew layout and will not build as-is on a stock machine — most notably -I/opt/homebrew/opt/llvm/include does not contain omp.h (it lives in a version-specific clang resource directory), and Apple's clang++ rejects -fopenmp outright:

clang++: error: unsupported option '-fopenmp'

The following recipe is verified working on macOS (Apple Silicon, Darwin 25.6, Open MPI 5.x + Homebrew LLVM).

macOS (Apple Silicon) — tested

# 1. Install MPI and an OpenMP-capable compiler.
#    Apple clang cannot do OpenMP, so Homebrew's LLVM provides both
#    clang++ and libomp. BLAS/LAPACK come from Apple's Accelerate
#    framework, which is already part of macOS — no openblas needed.
brew install open-mpi llvm

Then replace deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization/Configuration with:

# Path to the SBD library
SBD_PATH=../..

# MPI C++ compiler. OMPI_CXX points Open MPI's wrapper at Homebrew's clang++,
# which (unlike Apple clang) supports -fopenmp.
CCCOM=OMPI_CXX=/opt/homebrew/opt/llvm/bin/clang++ mpicxx

# Build flags. omp.h lives in a version-specific clang resource dir, so let
# Homebrew's clang++ find it itself rather than hardcoding an -I path.
CCFLAGS= -std=c++17 -fopenmp -O3

# Link flags: OpenMP runtime from Homebrew LLVM, BLAS/LAPACK from the
# Accelerate framework that ships with macOS (no openblas install needed).
SYSLIB= -L/opt/homebrew/opt/llvm/lib -lomp -framework Accelerate

# 2. Build
cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization
make
cd ../../../..

A successful build prints the two compile/link lines and leaves a diag binary next to the Makefile. One ld: warning: ignoring duplicate libraries: '-lomp' is harmless.

Intel macOS: Homebrew's prefix is /usr/local rather than /opt/homebrew, so substitute it in both preceding paths.

Linux (Ubuntu/Debian) — untested

The following prerequisites cover MPI, OpenMP, and BLAS/LAPACK. GCC supports -fopenmp natively, so no compiler override is needed. This has not been verified on Linux:

sudo apt install -y libopenmpi-dev libomp-dev libblas-dev liblapack-dev

SBD_PATH=../..
CCCOM=mpicxx
CCFLAGS= -std=c++17 -fopenmp -O3
SYSLIB= -llapack -lblas

Run SqDRIFT to generate CI strings

If you already ran ./SqDRIFT through Steps 1–4, alphadets_from_sqd.txt already exists in the project root and you can skip this step. Otherwise, run it from the project root:

./SqDRIFT
# Writes alphadets_from_sqd.txt in the project root

Run the diagonalization

cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization

# The FCIDump lives in the shared sqdrift/fcidump_files/ directory, one level
# above this tutorial's project root, hence five levels up from here.
ln -sf ../../../../../fcidump_files/N2_sto_3g fcidump.txt

# alphadets_from_sqd.txt is written by ./SqDRIFT into the project root itself,
# which is four levels up.
ln -sf ../../../../alphadets_from_sqd.txt alphadets.txt

./diag \
--fcidump fcidump.txt \
--adetfile alphadets.txt \
--method 0 \
--iteration 100 \
--block 10 \
--tolerance 1e-8 \
| tee ../../../../sbd_output.txt

cd ../../../..

Parameters:

• --fcidump: Path to the FCIDump file containing the Hamiltonian
• --adetfile: Path to the alpha determinants file (CI strings from SqDRIFT)
• --method 0: Davidson diagonalization method
• --iteration 100: Maximum number of Davidson iterations
• --block 10: Block size for Davidson algorithm
• --tolerance 1e-8: Convergence tolerance for energy

 Elapsed time for helper construction 0.001512 (sec)
 Elapsed time for init 1e-06 (sec)
 Davidson iteration 0.0 (tol=0.1857493937033264): -107.493531425221
 Davidson iteration 0.1 (tol=0.0455699976495391): -107.5098150493379 -105.93368397386
 Davidson iteration 0.2 (tol=0.009291923190816917): -107.5108356028703 -106.1245977311231 -105.6823366883346
 Davidson iteration 0.3 (tol=0.001488674426419356): -107.5108711084406 -106.3219246016818 -105.7376912083761 -105.4392974334491
 Davidson iteration 0.4 (tol=0.0003188697654737439): -107.5108722237947 -106.5935152393201 -106.1108017566682 -105.4932621297472
 Davidson iteration 0.5 (tol=3.080273116687244e-05): -107.5108722587463 -106.7245930484005 -106.2313070085243 -105.6423237790271
 Davidson iteration 0.6 (tol=2.981303155550973e-06): -107.5108722590151 -106.8008116251279 -106.330411722324 -105.69489630447
 Davidson iteration 0.7 (tol=4.565075388900609e-07): -107.5108722590186 -106.8090569201597 -106.3345525819528 -106.0578325555555
 Davidson iteration 0.8 (tol=4.870689101320012e-08): -107.5108722590186 -106.8098939189322 -106.3973756614146 -106.1617813840531
 Davidson iteration 0.9 (tol=6.224330513865847e-09): -107.5108722590187 -106.8159869231416 -106.4315891597535 -106.2222021964013
 Elapsed time for davidson 0.031276 (sec)
 Elapsed time for diagonalization 0.031286 (sec)
 Elapsed time for mult 0.002511 (sec)
 Energy = -107.5108722590187
 Elapsed time for measurement 4.6e-05 (sec)
 Sample-based diagonalization: Energy = -107.5108722590187
 Sample-based diagonalization: density = [1.999989340109303,1.99999662423008,1.999199121034164,1.981735389885314,1.999825046204946,1.999847567448985,1.995347419300126,0.01136512984757364,0.009844780833945199,0.002849581105572152
 Sample-based diagonalization: carryover bitstrings = [], size = 0

The Energy = -107.5108722590187 line is the ground-state estimate over the sampled subspace.

How this differs from the Python example

The two examples differ in several respects at once, not only in how they treat bad-symmetry bitstrings, so they are not expected to produce identical energies:

Both versions use the same Hamiltonian, group its terms by electronic structure, weight each group by the mean absolute coefficient of its terms, start from the Hartree-Fock state, and sample 10 groups and 100 shots per circuit.

Each difference affects the sampled subspace, and therefore the energy, independently. The larger ensemble and the second, longer evolution time both broaden the configurations the Python example explores; sampling at more than one time gives SqDRIFT its Krylov-like coverage of the low-energy space. Filtering diagonal terms changes which operators the sampler can draw and their weights, so the two runs sample different distributions.

Postselection is a strict filter, and a shot with the wrong particle number is dropped outright. In the run above that discarded most of the data, 1137 of 10000 shots survived the (7,7) Hamming-weight check (about 11%), collapsing into just 49 distinct $\alpha$-determinants. Configuration recovery instead uses average orbital occupancies from a prior diagonalization to flip bits and repair such shots, recycling what postselection would discard. That feedback makes it inherently iterative: recover_configurations in qiskit-addon-sqd-hpc requires an avg_occupancies argument that exists only after a diagonalization has run.

This C++ run yields a subspace of 49 $\alpha$-determinants and an energy of -107.5109 Ha. The full $\alpha$-determinant space for seven electrons in 10 orbitals is only $\binom{10}{7} = 120$ strings, so feeding all 120 to the same diag binary gives the exact FCI energy -107.6482 Ha. The run therefore sits about 0.1373 Ha above exact while spanning 49 of the 120 $\alpha$-determinants. Since SQD is variational in the selected subspace, a smaller subspace can only raise the energy, so the incomplete subspace is the most direct reading of the gap.

This exact-FCI check is only possible because the example is small; it is not part of the SqDRIFT workflow, and for system sizes where SQD is actually needed, no such reference exists. It serves here only to quantify how much of the space this run recovered.

A comparison against the Python example is not controlled, so a difference should not be attributed to the absence of configuration recovery alone; note that circuit count, the single evolution time, diagonal-term filtering, and hardware noise all differ at once and all shape the subspace.

To close the gap, here are some possible next steps to try after completing this tutorial: raise shots, num_circuits, or ops_per_circuit; sample at more than one time_step; and add a recovery loop on top of Qiskit::addon::sqd::recover_configurations.