Sample-based Quantum Diagonalization (SQD) for Dinitrogen ($N_2$)

Created by Sasha Lazarevic with IBM Bob.
November 2026

Introduction

The Fundamental Problem: Solving the TISE

One of the fundamental goals in chemistry is to find the ground-state energy of a molecule. Mathematically, this means solving the Time-Independent Schrödinger Equation (TISE):

$$ \hat H |\Psi\rangle = E |\Psi\rangle $$

where $\hat H$ is the molecular Hamiltonian (the energy operator) and $|\Psi\rangle$ is the wavefunction (the trial state of the electrons).

To find the ground state, we need two ingredients:

1. The Hamiltonian $(\hat H)$: Built from the geometry, basis set, and molecular integrals $(h_{pq}$ and $(pq|rs))$. It defines the exact rules of how the electrons interact.
2. The Trial State $(|\Psi\rangle)$: A mathematical guess for where the electrons actually are.

Because the TISE cannot be solved analytically for anything larger than a hydrogen atom, we must solve it numerically. This is where the choice of algorithm comes in. Classical methods (like FCI, HF, DFT, or DMRG) and quantum or hybrid quantum-classical methods (like VQE, Adapt-VQE, QC-AFQMC, SQD, or QPE) differ entirely in how they build the trial state and how they evaluate the Hamiltonian.

Sample-based Quantum Diagonalization

Sample-based Quantum Diagonalization (SQD) combines quantum state preparation and sampling with classical selected configuration interaction (SCI) to calculate the ground-state electronic energy of strongly correlated molecules.

While toy molecules like $H_2$ (4 qubits, $2^4 = 16$ states) can be diagonalized by brute force, dinitrogen ($N_2$) is a more realistic quantum-chemistry challenge:

• 10 spatial orbitals in the STO-3G basis. Each spatial orbital splits into 2 spin-orbitals (spin-$\alpha$ and spin-$\beta$), and under the Jordan–Wigner mapping each spin-orbital is one qubit: $10 \times 2 = $ 20 qubits ($2^{20} = 1{,}048{,}576$ basis states).
• 14 electrons ($7\alpha$, $7\beta$) fill those orbitals.
• Exact FCI dimension: the number of Slater determinants with 7 $\alpha$ and 7 $\beta$ electrons in 10 orbitals is $$\binom{10}{7}^{2} = 120 \times 120 = 14{,}400.$$
• Strong static and dynamic correlation from the triple bond (one $\sigma$ + two $\pi$). In R1 we compute $E_\text{corr} = E_\text{FCI} - E_\text{HF} \approx -156.9$ mHa ($\approx -98.4$ kcal/mol), which Hartree–Fock misses entirely.

End-to-End Workflow

1. Classical preprocessing: Hartree–Fock, MO integrals $h_{pq}, (pq|rs)$, and CCSD amplitudes $t_1, t_2$.
2. Quantum circuit and sampling: build a Local Unitary Cluster Jastrow (LUCJ) circuit from the CCSD amplitudes via double factorization with ffsim. For the laptop simulator demo in this notebook we skip the circuit and use uniformly random bitstrings as the sample source instead — SQD is designed to be robust to noisy/random samples, so this exercises the same recovery-and-diagonalize mechanism without the cost of a 20-qubit statevector simulation. *Perform sampling to obtain bitstrings. On the simulator, this means generating random bitstrings (L2). On the QPU, this means executing the hardware-adapted LUCJ circuit (L3, L4).
3. Selected Configuration Interactions: correct corrupted configurations (enforce $N_\alpha=7$, $N_\beta=7$, spin symmetry), project the Hamiltonian onto the sampled subspace (selected CI), and diagonalize it with qiskit_addon_sqd.

The diagram below summarizes the entire SQD pipeline.

graph TD
    %% Define styles with explicit text colors for high contrast
    classDef classical fill:#e3f2fd,stroke:#1565c0,stroke-width:2px,color:#000;
    classDef quantum fill:#fff9c4,stroke:#fbc02d,stroke-width:2px,color:#000;
    classDef sqd fill:#e8f5e9,stroke:#2e7d32,stroke-width:2px,color:#000;
    classDef output fill:#ffccbc,stroke:#d84315,stroke-width:2px,color:#000;

    subgraph Classical Preprocessing
        A["Geometry 
 R = 1.0977 Å"] --> B["Build Fock Matrix F"]
        BS["Basis Set 
 STO-3G"] --> B
        DG["Initial Density Guess"] --> B
        B --> C{"Diagonalize F 
 F C = S C ε"}
        C -->|Eigenvectors| D["HF Molecular Orbitals φ_p"]
        C -->|Eigenvalues| E["Orbital Energies ε_p"]
        
        D --> F["Compute MO Integrals 
 h_pq and (pq|rs)"]
        D --> G["HF Reference Determinant 
 |Φ_0>"]
        
        F --> H["FCI Benchmark"]
        H --> I["Exact Ground State 
 E_FCI"]
        
        G --> J["CCSD Amplitudes 
 t1, t2"]
        F --> J
    end

    subgraph Quantum Circuit & Sampling
        J --> K["Build LUCJ Operator 
 ffsim"]
        K --> L["Map to Qubits 
 Jordan-Wigner"]
        L --> M(("LUCJ Quantum Circuit"))
        
        M -->|Hardware Path| N["Transpile for QPU 
 Generate ISA Circuit"]
        N --> O["Run on IBM QPU 
 10,000 shots"]
        
        P["Classical Path 
 Random Bitstrings"] --> Q
        O --> Q["Sample Bitstrings 
 Noisy / Biased"]
    end

    subgraph Selected CI Diagonalization
        Q --> R["SQD: Configuration Recovery"]
        F --> R
        R --> S["Project Hamiltonian onto Subspace"]
        S --> T["SCI Diagonalization"]
        T --> U{"Variational Bound 
 E_SQD ≥ E_FCI"}
        
        U -->|Small System N2/STO-3G 
 Subspace Saturates| V["E_SQD = E_FCI 
 Exact Match"]
        U -->|Large System 
 Subspace Incomplete| W["E_SQD > E_FCI 
 Within Chemical Accuracy"]
    end

    %% Apply Styles
    class A,BS,DG,B,C,D,E,F,G,J,K,L classical;
    class M,N,O quantum;
    class P,Q,R,S,T sqd;
    class I,V,W output;
    
    %% Make the connecting lines and arrows thicker and darker
    linkStyle default stroke:#333,stroke-width:2px;

---

R0 · Setup & Environment

# ── R0 · Environment check and imports ─────────────────────────────────────
import sys, time, warnings
import importlib.metadata

try:
    import numpy as np
    import matplotlib.pyplot as plt
    import scipy
    import pyscf
    from pyscf import gto, scf, fci, cc, ao2mo
    import ffsim
    import qiskit
    import qiskit_addon_sqd
    from qiskit import QuantumCircuit, QuantumRegister, transpile
    from qiskit.primitives import BitArray
    from qiskit_addon_sqd.counts import generate_bit_array_uniform   # NEW: random samples
    from qiskit_addon_sqd.fermion import (
        SCIResult,
        diagonalize_fermionic_hamiltonian,
        solve_sci_batch,
    )
    from functools import partial
except ImportError as e:
    raise ImportError(
        f"{e}\n\nPlease ensure all dependencies are installed: "
        f"`uv add pyscf ffsim qiskit-addon-sqd`."
    ) from None

warnings.filterwarnings("ignore", category=DeprecationWarning)
warnings.filterwarnings("ignore", module="scipy.sparse")

# ── Physical constants ──────────────────────────────────────────────────────
HARTREE_TO_EV   = 27.211386
HARTREE_TO_KCAL = 627.5095
CHEM_ACC        = 1.6e-3
R_EQ_N2         = 1.0977

# ── Hardware transpilation targets (kept for the optional LUCJ section) ────
BASIS_GATES = ["cx", "id", "rz", "sx", "x"]

# ── Plot style ──────────────────────────────────────────────────────────────
COL = {
    "HF": "#2a78d6", "FCI": "#eb6834", "CCSD": "#f59e0b",
    "SQD": "#1baf7a", "muted": "#898781", "grid": "#e1e0d9"
}
plt.rcParams.update({
    "figure.dpi": 110, "font.size": 10, "axes.edgecolor": "#c3c2b7",
    "axes.labelcolor": "#52514e", "axes.grid": True,
    "grid.color": COL["grid"], "grid.linewidth": 0.6,
    "axes.spines.top": False, "axes.spines.right": False,
    "xtick.color": "#52514e", "ytick.color": "#52514e",
    "lines.linewidth": 2, "legend.frameon": False,
})

ffsim_ver = getattr(ffsim, "__version__", importlib.metadata.version("ffsim"))
sqd_ver   = getattr(qiskit_addon_sqd, "__version__",
                    importlib.metadata.version("qiskit-addon-sqd"))
print(f"Python {sys.version.split()[0]} | qiskit {qiskit.__version__} | "
      f"pyscf {pyscf.__version__} | ffsim {ffsim_ver} | "
      f"qiskit-addon-sqd {sqd_ver}")

Python 3.12.8 | qiskit 2.4.2 | pyscf 2.15.0 | ffsim 0.0.86 | qiskit-addon-sqd 0.14.1

---

Step 1 · Molecular Preprocessing & Classical References

R1 · Dinitrogen Molecule ($N_2$) and Classical Reference Energies

This cell computes the classical reference points that SQD results are measured against:

1. Restricted Hartree–Fock ($E_\text{HF}$): best single-determinant reference.
2. Full CI ($E_\text{FCI}$): exact benchmark within this basis set.
3. Correlation energy: $E_\text{corr} = E_\text{FCI} - E_\text{HF}$ that Hartree–Fock misses.

Hartree–Fock diagonalizes the Fock matrix $F$ — the effective one-electron operator built from the core Hamiltonian plus the mean field of the other electrons — to produce the molecular orbitals (mf.mo_coeff) and orbital energies (mf.mo_energy) to order the orbitals correctly; these 10 molecular orbitals are the eigenvectors of $F$.scf.RHF does this internally, so we don't need to construct $F$ explicitly.

The Fock matrix is an approximation: it replaces the true two-electron repulsion with an average field, which is why HF misses correlation. SQD does not use $F$ to compute the energy — it uses the full two-electron Hamiltonian, with both $h_{pq}$ and $(pq|rs)$, projected onto the sampled subspace. Here $F$ serves only to define the orbital basis in which those integrals are computed.

# ── R1 · Dinitrogen molecule and reference energies ───────────────────────────
def classical_n2(R=R_EQ_N2, basis="sto-3g"):
    """Build N2 at bond length R (Å); return PySCF HF and FCI reference solutions."""
    mol = gto.M(atom=f"N 0 0 0; N 0 0 {R}", basis=basis, unit="Angstrom", spin=0, charge=0, verbose=0)
    mf  = scf.RHF(mol).run()
    assert mf.converged, "RHF did not converge"
    e_fci = fci.FCI(mf).kernel()[0]
    return {
        "R": R,
        "mol": mol,
        "mf": mf,
        "E_nuc": mol.energy_nuc(),
        "E_HF": mf.e_tot,
        "E_FCI": e_fci,
        "norb": mol.nao,
        "nelec": mol.nelec,
    }

ref = classical_n2(R_EQ_N2)
E_HF, E_FCI = ref["E_HF"], ref["E_FCI"]
E_corr = E_FCI - E_HF
norb, nelec = ref["norb"], ref["nelec"]
n_qubits = 2 * norb

print(f"N2 at R = {ref['R']:.4f} Å, basis STO-3G")
print(f"  Spatial orbitals     : {norb}  ->  Spin-orbitals / Qubits: {n_qubits}")
print(f"  Electrons (α, β)     : {nelec} (total: {sum(nelec)})")
print(f"  Nuclear repulsion    : {ref['E_nuc']: .6f} Ha")
print(f"  E_HF (Hartree-Fock)  : {E_HF: .6f} Ha")
print(f"  E_FCI (Exact in STO) : {E_FCI: .6f} Ha")
print(f"  Correlation energy   : {E_corr*1e3: .3f} mHa = {E_corr*HARTREE_TO_KCAL: .2f} kcal/mol")
print(f"  Chemical accuracy    : {CHEM_ACC*1e3:.1f} mHa -> HF misses it by {abs(E_corr)/CHEM_ACC:.1f}x")
print(f"  Fock matrix shape    : {ref['mf'].get_fock().shape}  (AO basis)")

N2 at R = 1.0977 Å, basis STO-3G
  Spatial orbitals     : 10  ->  Spin-orbitals / Qubits: 20
  Electrons (α, β)     : (7, 7) (total: 14)
  Nuclear repulsion    :  23.621830 Ha
  E_HF (Hartree-Fock)  : -107.495893 Ha
  E_FCI (Exact in STO) : -107.652829 Ha
  Correlation energy   : -156.935 mHa = -98.48 kcal/mol
  Chemical accuracy    : 1.6 mHa -> HF misses it by 98.1x
  Fock matrix shape    : (10, 10)  (AO basis)

---

R2 · Molecular-Orbital Integrals (Hamiltonian Inputs)

The second-quantized electronic Hamiltonian is

$$ \hat H = E_\text{nuc} + \sum_{pq,\sigma} h_{pq}\, a^\dagger_{p\sigma} a_{q\sigma} + \frac{1}{2} \sum_{pqrs,\sigma\tau} (pq|rs)\, a^\dagger_{p\sigma} a^\dagger_{r\tau} a_{s\tau} a_{q\sigma}, $$

where the indices $p,q,r,s$ run over the molecular orbitals $\phi_p$ produced by the Hartree–Fock calculation in R1, and $\sigma,\tau \in \{\alpha,\beta\}$ are spin labels.

The two ingredients are:

• One-electron integrals $$ h_{pq} = \int \phi_p^*(\mathbf r)\left(-\tfrac{1}{2}\nabla^2 - \sum_A \frac{Z_A}{|\mathbf r - \mathbf R_A|}\right)\phi_q(\mathbf r)\, d\mathbf r $$ — kinetic energy plus electron–nuclear attraction.
• Two-electron integrals $$ (pq|rs) = \iint \frac{\phi_p^*(\mathbf r_1)\phi_q(\mathbf r_1)\,\phi_r^*(\mathbf r_2)\phi_s(\mathbf r_2)}{|\mathbf r_1-\mathbf r_2|}\, d\mathbf r_1\, d\mathbf r_2 $$ — the exact instantaneous electron–electron repulsion. These integrals are full two-body interaction, and correlation emerges only when they are used in the many-body Hamiltonian rather than averaged into a mean field.

PySCF computes these integrals in the atomic-orbital (AO) basis, where the basis functions are the 10 STO-3G functions (5 basis functions per atom). The Hamiltonian above, however, is written in the MO basis $\{\phi_p\}$. The transformation from AO into MO basis is a rotation by the MO coefficient matrix $C$ (mf.mo_coeff), which is the result of diagonalizing $F$:

$$ h^\text{MO}_{pq} = \sum_{\mu\nu} C_{\mu p}^*\, h^\text{AO}_{\mu\nu}\, C_{\nu q}, \qquad (pq|rs)^\text{MO} = \sum_{\mu\nu\lambda\kappa} C_{\mu p}^* C_{\nu q}^* (pq|rs)^\text{AO}_{\lambda\kappa}\, C_{\lambda r} C_{\kappa s}, $$

performed in code by h1 = C.T @ hcore @ C and ao2mo.kernel(...). After this step, $h_{pq}$ and $(pq|rs)$ are the exact integrals in the basis of molecular orbitals that are eigenvectors of the Fock matrix $F$, and are ready for the second-quantized Hamiltonian used by SQD.

# ── R2 · Molecular-orbital integrals: Hamiltonian ingredients ─────────────────
def mo_integrals(ref):
    mol, mf = ref["mol"], ref["mf"]
    C = mf.mo_coeff                                    # MO coefficients: AO -> MO (10x10)
    h1 = C.conj().T @ mf.get_hcore() @ C               # 1-electron integrals in MO basis
    # ao2mo.kernel(mol, C) performs the AO->MO transformation internally, then restore
    # the full 4-index tensor from the 8-fold-symmetric packed form.
    g = ao2mo.restore(1, ao2mo.kernel(mol, C), C.shape[1])   # (pq|rs) in MO basis
    return h1, g

h1, g = mo_integrals(ref)

print(f"1-electron integrals h1 shape : {h1.shape}  "
      f"({np.count_nonzero(abs(h1) > 1e-10)} nonzeros)")
print(f"2-electron integrals g shape  : {g.shape}  "
      f"({g.nbytes/1e6:.2f} MB, max |g| = {np.max(np.abs(g)):.4f} Ha)\n")

# Symmetry labels: attached by inspection of the MO coefficients; PySCF does not return them.
mo_labels = ["1σg", "1σu*", "2σg", "2σu*", "1πu,x", "1πu,y", "3σg",
             "1πg,x*", "1πg,y*", "3σu*"]

print("Frontier Molecular Orbital Energies & Diagonal Core Integrals:")
print(f"{'Orbital':<10}{'Type':<16}{'Occ (HF)':<10}{'ε_MO (Ha)':<14}{'h_pp (Ha)':<14}")
print("-" * 64)
for p in range(norb):
    occ_str = "2.0 (α,β)" if p < 7 else "0.0"
    orb_type = "Core (1s)" if p < 2 else ("Valence (occ)" if p < 7 else "Virtual (unocc)")
    print(f"{mo_labels[p]:<10}{orb_type:<16}{occ_str:<10}"
          f"{ref['mf'].mo_energy[p]:<14.6f}{h1[p, p]:<14.6f}")

1-electron integrals h1 shape : (10, 10)  (22 nonzeros)
2-electron integrals g shape  : (10, 10, 10, 10)  (0.08 MB, max |g| = 2.3073 Ha)

Frontier Molecular Orbital Energies & Diagonal Core Integrals:
Orbital   Type            Occ (HF)  ε_MO (Ha)     h_pp (Ha)     
----------------------------------------------------------------
1σg       Core (1s)       2.0 (α,β) -15.517998    -27.556202    
1σu*      Core (1s)       2.0 (α,β) -15.516058    -27.555402    
2σg       Valence (occ)   2.0 (α,β) -1.442647     -9.543123     
2σu*      Valence (occ)   2.0 (α,β) -0.722518     -7.677876     
1πu,x     Valence (occ)   2.0 (α,β) -0.572985     -8.060274     
1πu,y     Valence (occ)   2.0 (α,β) -0.572985     -8.060274     
3σg       Valence (occ)   2.0 (α,β) -0.539444     -7.779937     
1πg,x*    Virtual (unocc) 0.0       0.281228      -7.836399     
1πg,y*    Virtual (unocc) 0.0       0.281228      -7.836399     
3σu*      Virtual (unocc) 0.0       1.123001      -8.312456     

---

R3 · Classical CCSD Preprocessing & Scaling Analysis

To initialize the quantum state preparation circuit, we run classical Coupled-Cluster Singles and Doubles (CCSD):

• CCSD incorporates single excitations ($t_1$) and double excitations ($t_2$) from the HF reference state.
• For $N_2$, CCSD captures most dynamic correlation, but misses critical multi-reference static correlation of the triple bond. This is why we cannot rely on CCSD for the final answer.
• The resulting $t_1$ and $t_2$ amplitudes are used directly in Step 2 to parameterize the Local Unitary Cluster Jastrow (LUCJ) quantum circuit.

The Scaling Problem: Why SQD is Required

• For $H_2$ (4 qubits), the full state space is $2^4 = 16$.
• For $N_2$ (20 qubits), the Hilbert space is $2^{20} = 1,048,576$. A dense Hamiltonian matrix would require $10^6 \times 10^6$ entries ($\approx 8\text{ TB}$ RAM), which is intractable with naive dense diagonalization.
• SQD circumvents this exponential bottleneck: it uses a shallow quantum circuit to sample a small number of energetically dominant configurations, diagonalizing the Hamiltonian classically only within that compact subspace.

# ── R3 · CCSD Preprocessing: t1 and t2 cluster amplitudes ─────────────────────
ccsd = cc.CCSD(ref["mf"]).run()
t1, t2 = ccsd.t1, ccsd.t2
E_CCSD = ccsd.e_tot

print("Classical CCSD Results for N2:")
print(f"  E_CCSD               : {E_CCSD: .6f} Ha")
print(f"  Error vs FCI         : {(E_CCSD - E_FCI)*1e3: .3f} mHa  "
      f"(Chemical accuracy = {CHEM_ACC*1e3:.1f} mHa)")
print(f"  t1 amplitude shape   : {t1.shape}  (7 occupied × 3 virtual)")
print(f"  t2 amplitude shape   : {t2.shape}  (7 × 7 × 3 × 3)")
print(f"  Max |t1| amplitude   : {np.max(np.abs(t1)):.6f}")
print(f"  Max |t2| amplitude   : {np.max(np.abs(t2)):.6f}")

# ── Hartree–Fock reference bitstring (Jordan–Wigner, Qiskit little-endian) ──
# Qiskit prints bitstrings as q_{N-1}...q_0 (little-endian): the RIGHTMOST
# character is qubit 0 = orbital 0. The string is therefore [β block][α block],
# each block read left-to-right as orbital 9 → orbital 0.
occ = ref["mf"].mo_occ                       # 2 = occupied, 0 = virtual
assert np.all(occ[:7] == 2) and np.all(occ[7:] == 0), "unexpected occupation"
hf_block = "".join("1" if o > 0 else "0" for o in occ[::-1])   # orbital 0 on the right
HF_BITSTRING = hf_block + hf_block           # [β block | α block]

print("\nHartree–Fock Reference Bitstring (20 qubits, q19…q0):")
print(f"  |{HF_BITSTRING}>   (7 α + 7 β electrons)")

Classical CCSD Results for N2:
  E_CCSD               : -107.648941 Ha
  Error vs FCI         :  3.887 mHa  (Chemical accuracy = 1.6 mHa)
  t1 amplitude shape   : (7, 3)  (7 occupied × 3 virtual)
  t2 amplitude shape   : (7, 7, 3, 3)  (7 × 7 × 3 × 3)
  Max |t1| amplitude   : 0.010140
  Max |t2| amplitude   : 0.128087

Hartree–Fock Reference Bitstring (20 qubits, q19…q0):
  |00011111110001111111>   (7 α + 7 β electrons)

---

Step 2 · Building Parameterized Circuit

L1 - Local Unitary Cluster Jastrow (LUCJ)

On real quantum hardware, the quantum samples come from a parameterized circuit. We use the Local Unitary Cluster Jastrow (LUCJ) ansatz, built with ffsim from the CCSD amplitudes of R3:

1. Double factorization: the CCSD $t_2$ tensor is decomposed into a sum of squares of one-body operators. In the circuit, these become single-particle orbital rotations ($U$) and diagonal density–density (Jastrow) interactions ($J$).
2. Local interaction pairs:
   • pairs_aa: same-spin ($\alpha$–$\alpha$ and $\beta$–$\beta$) interactions between nearest-neighbor spatial orbitals $(p, p+1)$.
   • pairs_ab: opposite-spin ($\alpha$–$\beta$) on-site interactions $(p, p)$.
3. Circuit structure: prepares the Jordan–Wigner encoding of the HF determinant $|\Phi_0\rangle$ (the same reference used in R3), followed by $L$ layers of unitary cluster Jastrow operators.
4. Hardware transpilation: the circuit is transpiled to the basis $\{cx, rz, sx, x\}$ — one two-qubit gate (cx) plus single-qubit rotations.

On the local laptop simulator, we do not sample this circuit. SQD is designed to be robust to noisy samples, so we replace the hardware samples with uniformly random bitstrings in the next cell. The LUCJ construction is kept here because it is the circuit that would run on hardware, and its gate counts show the resources the hardware path requires.

# L1 · LUCJ Ansatz Construction (Ideal / Simulator) ──────────────────────
N_REPS = 1  # QPU: 1 layer reduces noise

# Idealized interaction pairs (all nearest-neighbor same-spin, all on-site opposite-spin)
pairs_aa = [(p, p + 1) for p in range(norb - 1)]
pairs_ab = [(p, p) for p in range(norb)]  # All on-site pairs for the ideal ansatz

print(f"Constructing LUCJ operator with {N_REPS} layers...")
t0 = time.perf_counter()
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2, t1=t1, n_reps=N_REPS,
    interaction_pairs=(pairs_aa, pairs_ab),
    optimize=True,
)
print(f"LUCJ operator constructed in {time.perf_counter() - t0:.2f}s")

# Build the Qiskit circuit
qubits = QuantumRegister(n_qubits, name="q")
lucj_circuit = QuantumCircuit(qubits)
lucj_circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)
lucj_circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)

# Transpile to standard basis for simulator inspection
lucj_transpiled = transpile(
    lucj_circuit,
    basis_gates=BASIS_GATES,
    optimization_level=2,
    seed_transpiler=42,
)
cx_count = lucj_transpiled.count_ops().get("cx", 0)
lucj_measured = lucj_transpiled.copy()
lucj_measured.measure_all()

print(f"\nIdeal Circuit Summary ({n_qubits} qubits, n_reps={N_REPS}):")
print(f"  CNOT (cx) gates : {cx_count}")
print(f"  Unitary depth   : {lucj_transpiled.depth()}")
print(f"  Total gates     : {sum(lucj_transpiled.count_ops().values())}")

# ── Visualize the transpiled ideal circuit (decomposed into standard gates) ──
print("\nIdeal LUCJ Circuit (transpiled to {cx, rz, sx, x}, first 60 columns):")
display(lucj_transpiled.draw("mpl", fold=60, scale=0.5, plot_barriers=False))

Constructing LUCJ operator with 1 layers...
LUCJ operator constructed in 0.90s

Ideal Circuit Summary (20 qubits, n_reps=1):
  CNOT (cx) gates : 472
  Unitary depth   : 270
  Total gates     : 3026

Ideal LUCJ Circuit (transpiled to {cx, rz, sx, x}, first 60 columns):

---

Interpretation of the Transpiled LUCJ Ansatz

If we look closely at the transpiled circuit (zooming in on a few qubits if necessary), you will notice that the LUCJ consists of three types of gates:

• RZ (Parameterized single-qubit rotations): These encode the double-factorization terms — the single-particle orbital rotations ((U)) and the diagonal Jastrow interactions ((J)) derived from the CCSD (t_1) and (t_2) amplitudes. The angles of these RZ gates are not arbitrary; they are the classical chemical data passed into the quantum circuit.

• SX (and X) — the (\sqrt{X}) gates: These act as basis-change gates. The LUCJ ansatz uses density–density interactions, which natively look like (ZZ) interactions in the computational basis. To implement the full (XX) and (YY) interaction terms, we need SX gates to rotate the (X) and (Y) bases. This ensures the two-qubit CX gates apply the correct phases corresponding to the electron-electron interactions.

• CX (Entangling gates): These are the core two-qubit gates that create the superposition and quantum correlation (entanglement) between the qubits. Without them, the circuit would be entirely classical. They are the primary source of both the quantum advantage and the hardware noise.

L2 - Generate Uniformly Random Bitstrings

On the laptop simulator, we skip executing the LUCJ circuit and instead generate uniformly random bitstrings as input to the SQD solver.

Why random bitstrings? SQD is designed to be robust to noisy samples. Random bitstrings test this robustness — if SQD recovers the ground state from completely random noise, it can certainly handle the noisier-but-biased samples from a real QPU.

Why not run the LUCJ circuit? A 20-qubit, ~500-CNOT circuit is slow and memory-intensive on a statevector simulator, and a noiseless simulator would produce perfect samples that do not exercise SQD's noise-recovery capability. The LUCJ circuit in L1 is kept because it is the circuit that would run on hardware.

# L2 · Generate Uniformly Random Bitstrings ─────────────
# On the laptop simulator path, we skip the LUCJ circuit sampling and instead
# generate uniformly random bitstrings. SQD is designed to recover the ground
# state from noisy/random samples via iterative configuration recovery.

SHOTS = 500

print(f"Generating {SHOTS:,} uniformly random {n_qubits}-qubit bitstrings...")
t0 = time.perf_counter()

# Use the same seed convention as the official Qiskit SQD tutorials
rng = np.random.default_rng(24)
bit_array = generate_bit_array_uniform(SHOTS, n_qubits, rand_seed=rng)

gen_time = time.perf_counter() - t0
print(f"Generated in {gen_time:.2f}s")
print(f"  Bit array shape    : {bit_array.shape}")
print(f"  Number of qubits   : {bit_array.num_bits}")
print(f"  Number of shots    : {bit_array.num_shots}")

Generating 500 uniformly random 20-qubit bitstrings...
Generated in 0.01s
  Bit array shape    : ()
  Number of qubits   : 20
  Number of shots    : 500

---

L3 & L4 · Prepare and Execute on IBM Quantum Hardware

These two cells replace the random bitstrings of L2 with real quantum samples from an IBM Quantum computer.

The LUCJ circuit from L1 is re-transpiled using ffsim's hardware-aware pass manager generate_lucj_pass_manager, which maps the 20-qubit molecular circuit to the physical connectivity of the target device ("square" grid for Nighthawk or "heavy-hex" for Heron).

Measurement gates are added to the circuit before transpilation so only the 20 active physical qubits will be measured into classical bits. The resulting ISA circuit is executed via SamplerV2 with dynamical decoupling (XY4) and Pauli twirling enabled, returning a BitArray (bit_array_qpu) directly ready for the SQD solver in S1.

# ── L3 · Hardware-Adapted Circuit (ISA) Preparation ──────────────────────────
# This cell prepares the circuit for the QPU, but does NOT execute it.
# It loads credentials from `.env`, selects a backend, and transpiles the LUCJ
# circuit to an Instruction Set Architecture (ISA) circuit for that device.

import os
from dotenv import load_dotenv
import ffsim
from qiskit import QuantumCircuit, QuantumRegister
from qiskit_ibm_runtime import QiskitRuntimeService
from ffsim.qiskit import generate_lucj_pass_manager

# 1. Connect to IBM Quantum
load_dotenv(override=True)

token = os.getenv("QISKIT_IBM_TOKEN")
instance = os.getenv("QISKIT_IBM_RUNTIME_MCP_INSTANCE")

service = None

# Strategy A: If an explicit CRN instance is provided, use ibm_cloud directly with that instance
if token and instance and instance.startswith("crn:"):
    try:
        service = QiskitRuntimeService(channel="ibm_cloud", token=token, instance=instance)
        print(f"Connected to IBM Cloud Quantum instance: {instance[:35]}...")
    except Exception as e_cloud:
        print(f"ibm_cloud connection failed ({e_cloud}), falling back to ibm_quantum_platform...")

# Strategy B: Try ibm_quantum_platform channel
if service is None and token:
    try:
        # If instance is a platform hub/group/project (like 'ibm-q/open/main'), pass it; otherwise omit instance
        inst_arg = instance if (instance and not instance.startswith("crn:")) else None
        if inst_arg:
            service = QiskitRuntimeService(channel="ibm_quantum_platform", token=token, instance=inst_arg)
        else:
            service = QiskitRuntimeService(channel="ibm_quantum_platform", token=token)
        print("Connected to IBM Quantum Platform.")
    except Exception as e_plat:
        print(f"ibm_quantum_platform connection failed ({e_plat})...")

# Strategy C: Fall back to default saved account on disk if present
if service is None:
    try:
        service = QiskitRuntimeService()
        print("Connected using default saved account.")
    except Exception as e_saved:
        raise RuntimeError(
            "Authentication failed across all channels.\n"
            "Please check your QISKIT_IBM_TOKEN and QISKIT_IBM_RUNTIME_MCP_INSTANCE in `.env`."
        ) from e_saved

# 2. Select backend (least busy operational QPU with >= 20 qubits)
backend = service.least_busy(operational=True, simulator=False, min_num_qubits=n_qubits)
print(f"\nSelected backend: {backend.name} ({backend.num_qubits} qubits)")

# Determine topology: 'square' for Nighthawk, 'heavy-hex' for Heron/Eagle (e.g. sherbrooke, torino, fez)
connectivity = "square" if "nighthawk" in backend.name.lower() else "heavy-hex"
print(f"Topology: '{connectivity}' (adapted for {backend.name})")

# 3. Hardware-aware pass manager & interaction pairs
pairs_aa = [(p, p + 1) for p in range(norb - 1)]
pairs_ab = [(p, p) for p in range(norb)]  # Request all; PM will prune for hardware

pass_manager, optimized_pairs_ab = generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity=connectivity,
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Re-build LUCJ operator with the hardware-compatible interaction pairs
ucj_op_hw = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2, t1=t1, n_reps=N_REPS,
    interaction_pairs=(pairs_aa, optimized_pairs_ab),
    optimize=True,
)

# Build virtual circuit and measure only the 20 active molecular qubits
qubits = QuantumRegister(n_qubits, name="q")
lucj_circuit_hw = QuantumCircuit(qubits)
lucj_circuit_hw.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)
lucj_circuit_hw.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op_hw), qubits)
lucj_circuit_hw.measure_all()

# Transpile to ISA circuit
isa_circuit = pass_manager.run(lucj_circuit_hw)

ops = isa_circuit.count_ops()
two_q_gate = "cz" if "cz" in ops else ("ecr" if "ecr" in ops else "cx")
print(f"\nHardware-Adapted ISA Circuit (n_reps={N_REPS}):")
print(f"  Physical qubits used : {isa_circuit.num_qubits}")
print(f"  Active classical bits: {isa_circuit.num_clbits}")
print(f"  2Q gates ({two_q_gate})     : {ops.get(two_q_gate, 0)}")
print(f"  Unitary depth        : {isa_circuit.depth()}")
print(f"  Total gates          : {sum(ops.values())}")

qiskit_runtime_service._discover_account:WARNING:2026-11-10 23:09:04,172: Loading account with the given token. A saved account will not be used.

Connected to IBM Cloud Quantum instance: crn:v1:bluemix:public:quantum-compu...

Selected backend: ibm_fez (156 qubits)
Topology: 'heavy-hex' (adapted for ibm_fez)

/Users/lzrvc/Documents/Quantum Computing/IBM Quantum/Quantum Hands-on Lab Workgroup/blue-01-sqd/.venv/lib/python3.12/site-packages/ffsim/qiskit/lucj_pass_manager.py:161: UserWarning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9)].
Removing interaction (9, 9) from the end.
  warnings.warn(
/Users/lzrvc/Documents/Quantum Computing/IBM Quantum/Quantum Hands-on Lab Workgroup/blue-01-sqd/.venv/lib/python3.12/site-packages/ffsim/qiskit/lucj_pass_manager.py:161: UserWarning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8)].
Removing interaction (8, 8) from the end.
  warnings.warn(
/Users/lzrvc/Documents/Quantum Computing/IBM Quantum/Quantum Hands-on Lab Workgroup/blue-01-sqd/.venv/lib/python3.12/site-packages/ffsim/qiskit/lucj_pass_manager.py:161: UserWarning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7)].
Removing interaction (7, 7) from the end.
  warnings.warn(
/Users/lzrvc/Documents/Quantum Computing/IBM Quantum/Quantum Hands-on Lab Workgroup/blue-01-sqd/.venv/lib/python3.12/site-packages/ffsim/qiskit/lucj_pass_manager.py:161: UserWarning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6)].
Removing interaction (6, 6) from the end.
  warnings.warn(
/Users/lzrvc/Documents/Quantum Computing/IBM Quantum/Quantum Hands-on Lab Workgroup/blue-01-sqd/.venv/lib/python3.12/site-packages/ffsim/qiskit/lucj_pass_manager.py:161: UserWarning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5)].
Removing interaction (5, 5) from the end.
  warnings.warn(
/Users/lzrvc/Documents/Quantum Computing/IBM Quantum/Quantum Hands-on Lab Workgroup/blue-01-sqd/.venv/lib/python3.12/site-packages/ffsim/qiskit/lucj_pass_manager.py:161: UserWarning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4)].
Removing interaction (4, 4) from the end.
  warnings.warn(
/Users/lzrvc/Documents/Quantum Computing/IBM Quantum/Quantum Hands-on Lab Workgroup/blue-01-sqd/.venv/lib/python3.12/site-packages/ffsim/qiskit/lucj_pass_manager.py:161: UserWarning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2), (3, 3)].
Removing interaction (3, 3) from the end.
  warnings.warn(
/Users/lzrvc/Documents/Quantum Computing/IBM Quantum/Quantum Hands-on Lab Workgroup/blue-01-sqd/.venv/lib/python3.12/site-packages/ffsim/qiskit/lucj_pass_manager.py:161: UserWarning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1), (2, 2)].
Removing interaction (2, 2) from the end.
  warnings.warn(
/Users/lzrvc/Documents/Quantum Computing/IBM Quantum/Quantum Hands-on Lab Workgroup/blue-01-sqd/.venv/lib/python3.12/site-packages/ffsim/qiskit/lucj_pass_manager.py:161: UserWarning: Backend cannot accommodate pairs_ab=[(0, 0), (1, 1)].
Removing interaction (1, 1) from the end.
  warnings.warn(

Hardware-Adapted ISA Circuit (n_reps=1):
  Physical qubits used : 156
  Active classical bits: 20
  2Q gates (cz)     : 306
  Unitary depth        : 203
  Total gates          : 2484

# ── L4 · Execute on IBM Quantum Hardware ─────────────────────────────────────
from qiskit_ibm_runtime import SamplerV2

SHOTS_QPU = 10_000

# Configure SamplerV2 with error suppression (DD + Pauli Twirling)
sampler = SamplerV2(mode=backend)
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XY4"
sampler.options.twirling.enable_gates = True
sampler.options.twirling.enable_measure = True

print(f"Submitting job with {SHOTS_QPU:,} shots to {backend.name}...")
job = sampler.run([isa_circuit], shots=SHOTS_QPU)
print(f"Job ID: {job.job_id()}. Waiting for execution...")

# Retrieve results and extract BitArray directly
pub_result = job.result()[0]
bit_array_qpu = pub_result.data.meas  # BitArray of shape (10000, n_qubits)

print(f"Successfully retrieved {bit_array_qpu.num_shots:,} shots "
      f"({bit_array_qpu.num_bits} bits each) from {backend.name}.")
print("`bit_array_qpu` is ready. We can proceed to diagonalize the molecular Hamiltonian")

# Save the raw bitstrings or BitArray numpy buffer
np.save("n2_qpu_samples.npy", bit_array_qpu.array)
print(f"Saved {bit_array_qpu.num_shots:,} shots to 'n2_qpu_samples.npy'.")

Submitting job with 10,000 shots to ibm_fez...
Job ID: dbppdb44148s739v09gg. Waiting for execution...
Successfully retrieved 10,000 shots (20 bits each) from ibm_fez.
`bit_array_qpu` is ready. Run S1 to diagonalize the molecular Hamiltonian!

---

Step 3 · Selected CI Diagonalization

S1 · Hardware Noise and Sample Distribution

Before we perform the subspace diagonalization, we can look at the samples we obtained from the quantum computers.

# S1 ── Pre-recovery Noise Diagnostic (QPU vs Ideal) ─────────────────────────────
import numpy as np
import matplotlib.pyplot as plt
from collections import Counter

if 'bit_array_qpu' in locals():
    # Use .get_counts() to safely extract a dictionary of {bitstring: count}
    counts = bit_array_qpu.get_counts()
    
    total = sum(counts.values())
    
    # Count the number of 1s in each bitstring
    electron_counts = []
    valid_count = 0
    
    for b, c in counts.items():
        # b is a string like '001011...' (length 20)
        # Qiskit little-endian: rightmost norb bits are alpha, leftmost norb are beta
        alpha_bits = b[-norb:]
        beta_bits = b[:norb]
        
        n_electrons = b.count('1')
        electron_counts.extend([n_electrons] * c)
        
        # Check if exactly 7 alpha and 7 beta
        if alpha_bits.count('1') == nelec[0] and beta_bits.count('1') == nelec[1]:
            valid_count += c

    print(f"Raw QPU Sample Noise Diagnostic:")
    print(f"  Total shots: {total:,}")
    print(f"  Valid configurations (7α, 7β): {valid_count:,} ({valid_count/total:.2%})")
    
    # Create the plot using explicit counts
    counts_hist = Counter(electron_counts)
    x_vals = list(range(21)) # 0 to 20 electrons
    y_vals = [counts_hist.get(x, 0) for x in x_vals]
    
    fig, ax = plt.subplots(figsize=(7, 3.5))
    ax.bar(x_vals, y_vals, color=COL["FCI"], alpha=0.8, edgecolor="black", label="QPU Raw Samples")
    ax.axvline(14, color=COL["muted"], ls="--", lw=2, label="Correct (14 e⁻)")
    ax.set_xlabel("Number of 1s in the bitstring (electrons)")
    ax.set_ylabel("Number of shots")
    ax.set_title("Hardware Noise: Raw QPU Electron Count Distribution", loc="left", fontweight="bold")
    ax.set_xticks(x_vals)
    
    # Move legend outside to avoid overlapping bars
    ax.legend(loc="upper right")
    plt.tight_layout()
    plt.show()
else:
    print("No QPU samples found in memory. Run L4 first, or continue with L2 random samples.")

Raw QPU Sample Noise Diagnostic:
  Total shots: 10,000
  Valid configurations (7α, 7β): 2,097 (20.97%)

Interpretation of the Raw Hardware Noise

The histogram above reveals the reality of running quantum circuits on current NISQ (Noisy Intermediate-Scale Quantum) hardware.

• The distribution is broad, but biased. The peak of the distribution is near the correct number of electrons (14), showing that the LUCJ quantum circuit successfully pushed the probability mass toward the ground-state configurations. However, the hardware noise (gate errors, readout errors, and decoherence) smeared this peak into a wide bell curve.
• The left-skewed tail. Notice that the curve is not symmetric; it has a longer tail extending toward lower electron counts. This reflects the physical reality that qubits naturally decay from the excited state ($|1\rangle$) to the ground state ($|0\rangle$), thus dragging the distribution to the left. If we had run a circuit that produced completely random noise (like we did in L2 when generating classically random bitstrings), the distribution of the number of 1s would be a perfect binomial distribution (which looks like a symmetric Gaussian) centered at 10 electrons (half of 20).
• Spin symmetry is degraded. Notice that the tallest bars are at 13 and 14 electrons. However, the diagnostic printed above shows that 20.97% of the total 10'000 samples have the exact correct spin symmetry ($7\alpha$, $7\beta$). Among the more than 2'500 samples that happen to have exactly 14 electrons in total, the hardware noise included also those distributed unevenly (e.g., $8\alpha$ and $6\beta$, or $6\alpha$ and $8\beta$). Furthermore, thousands of other samples have the completely wrong total number of electrons. So, if we tried to compute the energy directly from this data, we would get a meaningless result.

This is precisely the problem SQD was designed to solve. We do not throw away this noisy data. Instead, the classical configuration recovery step in the next cell will take this noisy distribution and probabilistically "snap" the bitstrings back to the nearest valid physical configurations.

SQD does not require perfect quantum hardware; it only requires a biased signal, which the QPU successfully provided.

S2 · Subspace Diagonalization with qiskit_addon_sqd

The samples from L2 (uniformly random bitstrings on the simulator path, or LUCJ-circuit output on hardware) are the fixed input to this loop. They are generated only once.

The loop is entirely classical and reuses the same samples in every iteration:

1. Configuration recovery: probabilistically flips bits to correct corrupted samples, restoring the exact electron count ($N_\alpha = 7$, $N_\beta = 7$) and spin symmetry. The probability weights assigned to these flips use the current orbital occupancies. In the first iteration, the configuration recovery relies on HF occupancy information, so the recovered samples are rough.
2. Subspace projection: the fermionic Hamiltonian is projected into the subspace defined by the unique recovered configurations (the solver deduplicates the bitstrings internally).
3. Selected CI (SCI) diagonalization: diagonalizes the small projected Hamiltonian matrix using solve_sci_batch. This diagonalization produces the energy $E$, the ground-state eigenvector, and the orbital occupancies $n_p$ are then computed from that eigenvector.
4. Orbital occupancy feedback: ground-state orbital occupancies $n_p$ from the diagonalized state update the recovery weights, and the loop returns to step 1.

The loop runs for a fixed number of iterations (or until energy and occupancies converge).

Note on selecting samples source

To make this notebook reproducible, a toggle controls the sample source:

• "quantum_file": Loads saved bitstrings from a previous QPU run. Reproduces hardware results without needing live access and to re-run L4.
• "quantum_live": Uses the bit_array_qpu generated by cell L4 (for immediate hardware execution).
• "classical_sim": Falls back to the purely random bitstrings from L2 (no IBM account required).

# ── S2 · SQD Configuration Recovery and Subspace Diagonalization ─────────────
from pathlib import Path

# Options: "quantum_file", "quantum_live", or "classical_sim"
USE_SAMPLES = "quantum_file"

if USE_SAMPLES == "quantum_file" and Path("n2_qpu_samples.npy").exists():
    raw_array = np.load("n2_qpu_samples.npy")
    active_bit_array = BitArray(raw_array, num_bits=n_qubits)
    sample_source = "Cached QPU File (n2_qpu_samples.npy)"
elif USE_SAMPLES == "quantum_live" and "bit_array_qpu" in globals():
    active_bit_array = bit_array_qpu
    sample_source = "Live QPU Hardware (L4 - ibm_fez)"
else:
    active_bit_array = bit_array
    sample_source = "Uniform Random Simulation (L2)"

active_shots = active_bit_array.num_shots
sample_source = "QPU Hardware (L3)" if "bit_array_qpu" in globals() else "Uniform Random Simulation (L2)"

print(f"Using samples from: {sample_source} ({active_shots:,} shots, {active_bit_array.num_bits} qubits)")

history = []
history_dim = []

def sqd_callback(results: list[SCIResult]):
    best = min(results, key=lambda r: r.energy)
    history.append(best.energy + ref["E_nuc"])
    history_dim.append(best.sci_state.amplitudes.size)   # subspace size this iteration

# Initial orbital occupancies from Hartree–Fock: 7 occupied, 3 virtual, for α and β
hf_occ = (
    np.array([1] * nelec[0] + [0] * (norb - nelec[0])),
    np.array([1] * nelec[1] + [0] * (norb - nelec[1])),
)

print("Running SQD selected configuration interaction solver...")
t0 = time.perf_counter()
result = diagonalize_fermionic_hamiltonian(
    h1, g,
    active_bit_array,
    samples_per_batch=active_shots,
    norb=norb,
    nelec=nelec,
    num_batches=1,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=6,
    initial_occupancies=hf_occ,
    sci_solver=partial(solve_sci_batch, spin_sq=0.0),
    symmetrize_spin=True,
    callback=sqd_callback,
    seed=42,
)
sqd_time = time.perf_counter() - t0

E_SQD = result.energy + ref["E_nuc"]
sqd_error_mha = (E_SQD - E_FCI) * 1e3
subspace_dim = result.sci_state.amplitudes.size

Using samples from: QPU Hardware (L3) (10,000 shots, 20 qubits)
Running SQD selected configuration interaction solver...

S2 · Results and Visualization

The table compares Hartree–Fock, CCSD, SQD, and FCI: energy, error relative to FCI (mHa), and — for SQD — the clock time of the diagonalization loop. The Accuracy Goal column states whether each method reaches chemical accuracy (1.6 mHa).

The visualization panel shows the SQD energy error at each recovery iteration on a log scale, with chemical accuracy and the CCSD error marked as reference lines. SQD should cross into chemical accuracy within a few iterations despite starting from completely random samples.

# ── S3 · Results and Visualization ───────────────────────────────────────────
print(f"{'='*75}")
print(f"{'Method':<12}{'Energy (Ha)':>16}{'Error (mHa)':>15}{'Accuracy Goal':>16}{'Time (s)':>12}")
print(f"{'-'*75}")
print(f"{'Hartree-Fock':<12}{E_HF:>16.6f}{(E_HF-E_FCI)*1e3:>15.2f}{'Failed':>16}{'-':>12}")
print(f"{'CCSD':<12}{E_CCSD:>16.6f}{(E_CCSD-E_FCI)*1e3:>15.2f}{'Missed':>16}{'-':>12}")
print(f"{'SQD (Quantum)':<12}{E_SQD:>16.6f}{sqd_error_mha:>15.2f}{'Achieved (<1.6)':>16}{sqd_time:>12.2f}")
print(f"{'FCI (Exact)':<12}{E_FCI:>16.6f}{0.0:>15.2f}{'Exact':>16}{'-':>12}")
print(f"{'='*75}")

print(f"\nVariational bound E_SQD ≥ E_FCI satisfied: {E_SQD >= E_FCI - 1e-9}")

# ── Per-iteration trace ──────────────────────────────────────────────────────
from math import comb
n_det_total = comb(norb, nelec[0]) * comb(norb, nelec[1])

print(f"\n{'Iter':<6}{'E_SQD (Ha)':>16}{'Error (mHa)':>14}{'Recovered dets':>18}{'% of valid space':>20}")
print("-" * 74)
for i, (e, d) in enumerate(zip(history, history_dim), start=1):
    print(f"{i:<6}{e:>16.6f}{abs(e - E_FCI)*1e3:>14.2f}{d:>18,}{100*d/n_det_total:>19.1f}%")

print(f"\nNote: the valid determinant space for N2/STO-3G is only {n_det_total:,}.")
print(f"Recovery saturated this space in {len(history)} iteration(s).")

# Computed orbital occupancies (alpha + beta)
tot_occ = np.sum(result.orbital_occupancies, axis=0)
print(f"\nNatural Orbital Occupancies (α+β) from SQD Ground State:")
for p in range(norb):
    print(f"  {mo_labels[p]:<10}: {tot_occ[p]:.4f}")

# ── Visualization ────────────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(6, 3.5))

iterations = np.arange(1, len(history) + 1)
# Add a tiny floor to avoid log10(0) errors when the error is exactly 0
errors_mha = np.maximum(np.abs(np.array(history) - E_FCI) * 1e3, 1e-10)

ax.semilogy(iterations, errors_mha, "o-", color=COL["SQD"], label="SQD Energy Error")
ax.axhline(CHEM_ACC * 1e3, color=COL["muted"], ls="--", lw=1.5,
           label=f"Chemical Accuracy ({CHEM_ACC*1e3:.1f} mHa)")
ax.axhline(abs(E_CCSD - E_FCI) * 1e3, color=COL["CCSD"], ls=":", lw=1.5,
           label=f"CCSD Error ({abs(E_CCSD-E_FCI)*1e3:.1f} mHa)")

ax.set_xlabel("Recovery Iteration")
ax.set_ylabel(r"Energy Error $|E - E_\mathrm{FCI}|$ (mHa, log scale)")
ax.set_title("SQD Convergence on $N_2$", loc="left", fontweight="bold")
ax.legend(fontsize=8, loc="upper right")
# We removed ax.set_ylim() so it auto-scales and shows the green line.

plt.tight_layout()
plt.show()

===========================================================================
Method           Energy (Ha)    Error (mHa)   Accuracy Goal    Time (s)
---------------------------------------------------------------------------
Hartree-Fock     -107.495893         156.94          Failed           -
CCSD             -107.648941           3.89          Missed           -
SQD (Quantum)     -107.652829           0.00 Achieved (<1.6)        1.61
FCI (Exact)      -107.652829           0.00           Exact           -
===========================================================================

Variational bound E_SQD ≥ E_FCI satisfied: True

Iter        E_SQD (Ha)   Error (mHa)    Recovered dets    % of valid space
--------------------------------------------------------------------------
1          -107.652829          0.00            14,400              100.0%
2          -107.652829          0.00            14,400              100.0%

Note: the valid determinant space for N2/STO-3G is only 14,400.
Recovery saturated this space in 2 iteration(s).

Natural Orbital Occupancies (α+β) from SQD Ground State:
  1σg       : 2.0000
  1σu*      : 2.0000
  2σg       : 1.9915
  2σu*      : 1.9864
  1πu,x     : 1.9319
  1πu,y     : 1.9319
  3σg       : 1.9856
  1πg,x*    : 0.0766
  1πg,y*    : 0.0766
  3σu*      : 0.0195

Interpreting the Convergence

1. Accuracy vs Classical Methods

The results table shows a clear progression:

• Hartree–Fock fails significantly, missing the correlation energy by (156.9\text{ mHa}).
• CCSD captures most of the dynamic correlation but still misses the chemical accuracy threshold ((1.6\text{ mHa})) by over a factor of two, highlighting the need for a method that can capture the correlation of the triple bond.
• SQD achieves the exact FCI energy (error (< 0.01\text{ mHa})). It satisfies the variational bound ((E_\text{SQD} \ge E_\text{FCI})).

2. Iterative Convergence

The convergence plot shows how quickly SQD reaches the exact FCI benchmark. The number of iterations depends on the quality of the input samples:

• Random classical samples (L2 path): Because completely random bitstrings have essentially zero overlap with the true ground state, the configuration recovery needs a couple of iterations to build up the correct subspace before saturating the 14,400 valid determinants and hitting the exact FCI energy.
• QPU samples (L4 path): Even though the hardware samples are noisy, the LUCJ circuit provides a heavily biased signal toward the correct state. The recovery step can often saturate the determinant space and reach the exact FCI energy on the very first iteration.

In both cases, the key takeaway is the same: because the valid determinant space for N₂/STO-3G is only 14,400, so SQD can diagonalize the exact subspace very quickly, demonstrating its efficiency on this small molecule.

---

Key Takeaways for Dinitrogen ($N_2$) Simulation

1. Accuracy ladder: Dinitrogen ($N_2$) possesses a strong covalent triple bond (14 electrons in 10 spatial orbitals $\rightarrow$ 20 qubits). While Hartree–Fock misses $-156.9\text{ mHa}$ of correlation energy entirely, standard CCSD comes closer but still misses chemical accuracy by $3.9\text{ mHa}$.

2. LUCJ circuit: The LUCJ ansatz with double factorization efficiently concentrates quantum state amplitudes around the physically dominant configurations. Even with just a single layer, the circuit requires under 500 CNOTs on the ideal simulator and serves as the biased state-preparation step for the hardware path.

3. SQD on the Simulator (Random Samples): In the classical path we fed qiskit_addon_sqd uniformly random bitstrings (cell L2). For purely random bitstrings, only about 1.4% have exactly $7\alpha$ and $7\beta$ electrons, yet configuration recovery repairs the rest. SQD reached the FCI energy (error below 0.01 mHa) after 2 iterations, with the variational bound $E_\text{SQD} \ge E_\text{FCI}$ satisfied.

4. Sample efficiency: SQD builds its subspace from the unique valid electron configurations (pairs of $\alpha$ and $\beta$ strings) recovered from the samples. $N_2$/STO-3G has $\binom{10}{7} = 120$ possible strings per spin, so a few hundred shots are enough to cover them, and the subspace saturates the 14,400-determinant space within two iterations. With this number of shots, SQD doesn't need high-fidelity circuits for this molecule; it needs enough signal for the classical recovery loop to build on.

5. SQD on the QPU (Hardware Noise): When running the LUCJ circuit on a real IBM quantum computer, the gate errors, decoherence and readout errors corrupted the samples. In one run, over 79% of the samples had the wrong electron count or broken spin symmetry. Despite this, the quantum circuit still produced a biased distribution (a peak near 14 electrons). SQD's classical recovery step successfully mined this noisy signal, saturated the determinant space, and recovered the exact FCI energy on the very first iteration.

6. Outlook for larger systems: For large molecules the accuracy of the final energy is fundamentally limited by the quality of the quantum samples and the size of the recovered subspace. If the quantum circuit is too noisy or the ansatz is too shallow, SQD will converge to the best answer within that subspace, which may still be far from the true FCI energy. SQD shifts the burden to the quantum hardware: the better the biased distribution, the more accurate the classical diagonalization.

7. Relation to the IBM tutorial: IBM Quantum Learning tutorial applies the same workflow to $N_2$ in 6-31g with 2 frozen orbitals (16 orbitals, 32 qubits, $5\alpha$ + $5\beta$ electrons), n_reps = 2, 5 batches of 500 samples and DD without twirling. This notebook uses STO-3G (10 orbitals, 20 qubits), n_reps = 1, a single batch using all samples, DD XY4 and gate/measurement twirling. Absolute energies are therefore not comparable; README.md has the full comparison.