Scale SQD chemistry workflows with the SBD eigensolver
This guide shows how to use the SBD eigensolver 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.
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:
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.
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,
)print(f"Exact energy: {exact_energy}")
print(f"Estimated energy: {result.energy + nuclear_repulsion_energy}")