---
title: Scale SQD chemistry workflows with the SBD eigensolver
description: Scale SQD chemistry workflows with the SBD eigensolver for the latest version of Sample-based quantum diagonalization (SQD)
source: https://quantum.cloud.ibm.com/docs/en/addons/qiskit-addon-sqd/guides/integrate-sbd
---

# Scale SQD chemistry workflows with the SBD eigensolver

This guide shows how to use the [SBD eigensolver](https://github.com/Qiskit/sbd-eigensolver-python) as an alternative selected configuration interaction (SCI) solver to diagonalize larger fermionic problems beyond what the default PySCF solver can support. SBD (Selected Basis Diagonalization) runs in-process, is parallelized with MPI, and offers both CPU and CUDA backends. For information on how to install and use it, see the [SBD documentation](https://github.com/Qiskit/sbd-eigensolver-python).

SBD compiles a C++ extension against the MPI and BLAS libraries available on the machine, so you must install both before you install the `sbd` extra:

```bash
pip install qiskit-addon-sqd[sbd]
```

Because `solve_sci_batch` already matches the solver interface that `diagonalize_fermionic_hamiltonian` expects, you can pass it directly to `sci_solver` after you bind its configuration options with `functools.partial`. To run on a CUDA-capable GPU instead of the CPU, pass `device_config=DeviceConfig.gpu()` from `sbd.device_config`.

For more details on the SQD code used in this example, see the [chemistry Hamiltonian tutorial](/docs/tutorials/sample-based-quantum-diagonalization).

```python
from functools import partial

import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
from sbd.sbd_solver import solve_sci_batch

# Specify molecule properties
num_orbitals = 16
num_elec_a = num_elec_b = 5
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)

# Compute exact energy
exact_energy = cas.run().e_tot

# Create a seed to control randomness throughout this workflow
rng = np.random.default_rng(24)


# Generate random samples
bit_array = generate_bit_array_uniform(
    10_000, num_orbitals * 2, rand_seed=rng
)

# Bind the SBD solver options. method=0 selects the Davidson eigensolver.
sbd_solver = partial(
    solve_sci_batch,
    sbd_config={"method": 0, "eps": 1e-5, "max_it": 10, "max_nb": 10},
)

# Run SQD
result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=300,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=5,
    max_iterations=5,
    sci_solver=sbd_solver,
    symmetrize_spin=True,
    seed=rng,
)
```

```python
print(f"Exact energy: {exact_energy}")
print(f"Estimated energy: {result.energy + nuclear_repulsion_energy}")
```