Scale SQD chemistry workflows with Fulqrum

Fulqrum is a software package that computes eigenvalues and eigenvectors of Hamiltonians projected into a subspace. It constructs a portable representation of the projected Hamiltonian that you can pass to an eigenvalue solver of your choice, including standard off-the-shelf solvers. These solvers can provide a much faster alternative to the default PySCF solver that the SQD addon uses for diagonalization. For information on how to install and use Fulqrum, see the Fulqrum documentation.

For more details on the SQD code used in this example, see the chemistry Hamilton tutorial.

import fulqrum as fq
import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
import scipy.sparse.linalg
from fulqrum.convert.integrals import integrals_to_fq_fermionic_op
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import (
    SCIResult,
    SCIState,
    diagonalize_fermionic_hamiltonian,
)

# 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
)

# Convert Hamiltonian to Fulqrum format
fermionic_op = integrals_to_fq_fermionic_op(
    one_body_integrals=hcore, two_body_integrals=eri
)
fulqrum_operator = fermionic_op.extended_jw_transformation()

# Initialize a Hamiltonian with a trivial subspace to be updated later
Hsub = fq.SubspaceHamiltonian(
    fulqrum_operator,
    fq.Subspace([["0" * num_orbitals], ["0" * num_orbitals]]),
)


# Define fulqrum solver
def solve_sci_batch(ci_batches, _h1, _h2, norb, nelec, tol=1e-12):
    # Ignore the _h1 and _h2 inputs because we store the Hamiltonian separately
    results = []
    for ints_a, ints_b in ci_batches:
        strs_a = [format(x, f"0{norb}b") for x in ints_a]
        strs_b = [format(x, f"0{norb}b") for x in ints_b]
        subspace = fq.Subspace([strs_a, strs_b])
        Hsub.update_subspace(subspace)
        Hsub_csr_linop = Hsub.to_csr_linearoperator_fast()
        diag = Hsub.diagonal_vector()
        v0 = np.zeros(len(subspace), dtype=Hsub.dtype)
        v0[np.argmin(diag)] = 1.0
        eigs, vecs = scipy.sparse.linalg.eigsh(
            Hsub_csr_linop, k=1, which="SA", tol=tol, v0=v0
        )
        vec = vecs[:, 0]
        occs = subspace.get_orbital_occupancies(
            probs=np.abs(vec) ** 2, norb=norb
        )
        # SQD ordering and Fulqrum ordering are transpose of each other.
        # ascontiguousarray() is needed for matching C contiguous ordering inside PySCF
        state = SCIState(
            amplitudes=np.ascontiguousarray(
                vec.reshape(len(strs_b), len(strs_a)).T
            ),
            ci_strs_a=ints_a,
            ci_strs_b=ints_b,
            norb=norb,
            nelec=nelec,
        )
        results.append(
            SCIResult(
                energy=eigs[0], sci_state=state, orbital_occupancies=occs
            )
        )
    return results


# 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=solve_sci_batch,
    symmetrize_spin=True,
    seed=rng,
)

Output:

converged SCF energy = -108.835236570775
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

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

Output:

Exact energy: -109.04667177808028
Estimated energy: -109.0340266755937

Use a different eigensolver

Because Fulqrum exposes the projected Hamiltonian as a standard scipy.sparse.linalg.LinearOperator, you can diagonalize it with any eigensolver of your choosing rather than being tied to scipy.sparse.linalg.eigsh. This makes it straightforward to substitute a solver that is better suited to your problem size or hardware.

As an example, this guide demonstrates PRIMME, a high-performance eigensolver library. Its Python interface, primme.eigsh, follows a calling convention similar to scipy.sparse.linalg.eigsh, so you can use it as a direct replacement.

The only change from the previous solver is the call to primme.eigsh instead of scipy.sparse.linalg.eigsh.

import primme


# Define a fulqrum solver that uses PRIMME to diagonalize the projected Hamiltonian
def solve_sci_batch_primme(ci_batches, _h1, _h2, norb, nelec, tol=1e-12):
    # Ignore the _h1 and _h2 inputs because we store the Hamiltonian separately
    results = []
    for ints_a, ints_b in ci_batches:
        strs_a = [format(x, f"0{norb}b") for x in ints_a]
        strs_b = [format(x, f"0{norb}b") for x in ints_b]
        subspace = fq.Subspace([strs_a, strs_b])
        Hsub.update_subspace(subspace)
        Hsub_csr_linop = Hsub.to_csr_linearoperator_fast()
        diag = Hsub.diagonal_vector()
        v0 = np.zeros(len(subspace), dtype=Hsub.dtype)
        v0[np.argmin(diag)] = 1.0
        # primme.eigsh follows the same calling convention as scipy.sparse.linalg.eigsh,
        # except that it expects the initial guess v0 as a 2D array of column vectors.
        eigs, vecs = primme.eigsh(
            Hsub_csr_linop, k=1, which="SA", tol=tol, v0=v0.reshape(-1, 1)
        )
        vec = vecs[:, 0]
        occs = subspace.get_orbital_occupancies(
            probs=np.abs(vec) ** 2, norb=norb
        )
        # SQD ordering and Fulqrum ordering are transpose of each other.
        # ascontiguousarray() is needed for matching C contiguous ordering inside PySCF
        state = SCIState(
            amplitudes=np.ascontiguousarray(
                vec.reshape(len(strs_b), len(strs_a)).T
            ),
            ci_strs_a=ints_a,
            ci_strs_b=ints_b,
            norb=norb,
            nelec=nelec,
        )
        results.append(
            SCIResult(
                energy=eigs[0], sci_state=state, orbital_occupancies=occs
            )
        )
    return results


# Run SQD using the PRIMME-based solver
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=solve_sci_batch_primme,
    symmetrize_spin=True,
    seed=rng,
)

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

Output:

Exact energy: -109.04667177808028
Estimated energy: -109.03937429196165