Scale SQD chemistry workflows with Fulqrum

Fulqrum is a software package for enabling the computation of eigenvalues and eigenvectors of Hamiltonians projected into a subspace. It achieves this by constructing a portable representation of the projected Hamiltonian that can then be passed 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 used for diagonalization in the SQD addon. For information on how to install and use Fulqrum, visit the docs.

For more details on the SQD code used in this example, check out tutorial 1.

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

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}")

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 <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.LinearOperator.html>`__, you are free to diagonalize it with any eigensolver of your choosing rather than being tied to scipy.sparse.linalg.eigsh. This makes it easy to swap in a solver that is better suited to your problem size or hardware.

As an example, we demonstrate using PRIMME, a high-performance eigensolver library. Its Python interface, primme.eigsh, follows a similar calling convention as scipy.sparse.linalg.eigsh, so it can be used as a drop-in replacement.

The only change from the previous solver is that we call 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}")

Exact energy: -109.04667177808028
Estimated energy: -109.03937429196165