Sample-based quantum diagonalization of a chemistry Hamiltonian

Usage estimate: under one minute on a Heron r2 processor (NOTE: This is an estimate only. Your runtime might vary.)

Learning outcomes

After going through this tutorial, users should understand:

How to use the SQD Qiskit addon to approximate the ground state energy of a molecular system using bitstrings sampled from a quantum processing unit (QPU).
How to use ffsim to construct a local unitary cluster Jastrow (LUCJ) circuit for quantum chemistry simulation.

Prerequisites

We recommend that users familiarize themselves with the following topics before going through this tutorial:

Quantum chemistry and second quantization
Using the Sampler primitive to sample from quantum circuits

Background

In this tutorial, we show how to post-process noisy quantum samples to approximate the ground state of the nitrogen molecule $\text{N}_2$ at equilibrium bond length, by using the SQD Qiskit addon to implement the sample-based quantum diagonalization (SQD) algorithm. More details on the software can be found in the corresponding documentation, including a simple example to get started.

This tutorial is recommended for users who are familiar with quantum chemistry: specifically, familiarity with finding ground state energies of a molecule. For a detailed walkthrough on the workflow, refer to the quantum diagonalization algorithm course.

SQD is a technique for finding eigenvalues and eigenvectors of quantum operators, such as a quantum system Hamiltonian, by using quantum and distributed classical computing together. Classical distributed computing is used to process samples obtained from a quantum processor, and to project and diagonalize a target Hamiltonian in a subspace that they span. An SQD-based workflow has the following steps:

Choose a circuit ansatz and apply it on a quantum computer to a reference state (in this case, the Hartree-Fock state).
Sample bitstrings from the resulting quantum state.
Run the self-consistent configuration recovery procedure on the bitstrings to obtain the ground state approximation.

SQD is known to work well when the target eigenstate is sparse: the wave function is supported in a set of basis states $\mathcal{S} = \{|x\rangle \}$ whose size does not increase exponentially with the size of the problem.

Quantum chemistry

The Hamiltonian of a molecular system can be written as

\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \frac12 \sum_{ \substack{prqs\\\sigma\tau} } h_{prqs} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma},

where the $h_{pr}$ and $h_{prqs}$ are complex numbers called molecular integrals that can be calculated from the specification of the molecule using a computer program. In this tutorial, we compute the integrals using the PySCF software package.

For details about how the molecular Hamiltonian is derived, consult a textbook on quantum chemistry (for example, Modern Quantum Chemistry by Szabo and Ostlund). For a high-level explanation of how quantum chemistry problems are mapped onto quantum computers, check out the lecture Mapping Problems to Qubits from Qiskit Global Summer School 2024.

Local unitary cluster Jastrow (LUCJ) ansatz

SQD requires a quantum circuit ansatz to draw samples from. In this tutorial, we'll use the local unitary cluster Jastrow (LUCJ) ansatz due to its combination of physical motivation and hardware-friendliness. We'll use ffsim to construct the ansatz circuit.

The LUCJ ansatz adapts to QPUs with restricted qubit connectivity. The spin orbitals are mapped to qubits such that the ansatz does not require routing with SWAP gates. IBM® hardware has a heavy-hex lattice qubit topology, in which case we can adopt a "zig-zag" pattern, depicted below. In this pattern, orbitals with the same spin are mapped to qubits with a line topology (red and blue circles), and a connection between orbitals of different spin is present at every 4th spatial orbital, with the connection being facilitated by an ancilla qubit (purple circles).

Qubit mapping diagram for the LUCJ ansatz on a heavy-hex lattice

Self-consistent configuration recovery

The self-consistent configuration recovery procedure is designed to extract as much signal as possible from noisy quantum samples. Because the molecular Hamiltonian conserves particle number and spin Z, it makes sense to choose a circuit ansatz that also conserves these symmetries. When applied to the Hartree-Fock state, the resulting state has a fixed particle number and spin Z in the noiseless setting. Therefore, the spin- $\alpha$ and spin- $\beta$ halves of any bitstring sampled from this state should have the same Hamming weight as in the Hartree-Fock state. Due to the presence of noise in current quantum processors, some measured bitstrings will violate this property. A simple form of postselection would discard these bitstrings, but this is wasteful because those bitstrings might still contain some signal. The self-consistent recovery procedure attempts to recover some of that signal in post-processing. The procedure is iterative and requires as input an estimate of the average occupancies of each orbital in the ground state, which is first computed from the raw samples. The procedure is run in a loop, and each iteration has the following steps:

For each bitstring that violates the specified symmetries, flip its bits with a probabilistic procedure designed to bring the bitstring closer to the current estimate of the average orbital occupancies, to obtain a new bitstring.
Collect all of the old and new bitstrings that satisfy the symmetries, and subsample subsets of a fixed size, chosen in advance.
For each subset of bitstrings, project the Hamiltonian into the subspace spanned by the corresponding basis vectors (see the previous section for a description of these basis vectors), and compute a ground state estimate of the projected Hamiltonian on a classical computer.
Update the estimate of the average orbital occupancies with the ground state estimate with the lowest energy.

SQD workflow diagram

The SQD workflow is depicted in the following diagram:

Requirements

Before starting this tutorial, ensure that you have the following installed:

Qiskit SDK v1.0 or later, with visualization support
Qiskit Runtime v0.22 or later (pip install qiskit-ibm-runtime)
SQD Qiskit addon v0.11 or later (pip install qiskit-addon-sqd)
ffsim v0.0.75 or later (pip install ffsim)

Setup

import math

import ffsim
import matplotlib.pyplot as plt
import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.primitives import StatevectorSampler
from qiskit.providers.fake_provider import GenericBackendV2
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

Small-scale simulator example

In this tutorial, we will find an approximation to the ground state of a nitrogen molecule at equilibrium. We first use a small STO-6G basis set so that we can simulate the experiment and make sure it's working.

Step 1: Map classical inputs to a quantum problem

First, we specify the molecule and its properties.

# Specify molecule properties
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="sto-6g",
    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()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
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), norb)

# Compute exact energy using FCI
reference_energy = cas.run().e_tot

print(f"norb = {norb}")
print(f"nelec = {nelec}")

Output:

converged SCF energy = -108.464957764796
CASCI E = -108.595987350986  E(CI) = -32.4115475088426  S^2 = 0.0000000
norb = 8
nelec = (5, 5)

Before constructing the LUCJ ansatz circuit, we first perform a CCSD calculation in the following code cell. The $t_1$ and $t_2$ amplitudes from this calculation will be used to initialize the parameters of the ansatz.

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

Output:

E(CCSD) = -108.5933309085008  E_corr = -0.1283731437052354

Now, we use ffsim to create the ansatz circuit. Since our molecule has a closed-shell Hartree-Fock state, we use the spin-balanced variant of the UCJ ansatz, UCJOpSpinBalanced. We set optimize=True in the from_t_amplitudes method to enable the "compressed" double-factorization of the $t_2$ amplitudes (see The local unitary cluster Jastrow (LUCJ) ansatz at ffsim's documentation for details).

Because the LUCJ ansatz adapts to the available connectivity of the QPU, we have to initialize the QPU backend before creating the ansatz. For now, we'll create a generic backend with a heavy-hex coupling map and a gate set that the LUCJ ansatz naturally decomposes to. Then, we'll use ffsim.qiskit.generate_lucj_pass_manager to create a pass manager that is specialized for transpiling the LUCJ ansatz to the given backend according to the "zig-zag" layout described in the background section on the LUCJ ansatz. This function uses a scoring heuristic to minimize the errors associated with the selected layout, which is important if your backend is a real QPU, or a simulator with a noise model. In addition to returning the pass manager, this function also returns the alpha-beta coupling pairs that can be implemented on the hardware. If not all pairs can be implemented, it emits a warning.

import warnings

from qiskit.transpiler import CouplingMap

warnings.formatwarning = lambda msg, *args, **kwargs: f"Warning: {msg}\n"

# Set ansatz properties
n_reps = 1
pairs_aa = [(p, p + 1) for p in range(norb - 1)]
pairs_ab = None  # Let generate_lucj_pass_manager determine the alpha-beta interactions

# Initialize backend
coupling_map = CouplingMap.from_heavy_hex(3)
backend = GenericBackendV2(
    coupling_map.size(),
    coupling_map=coupling_map,
    basis_gates=["cp", "xx_plus_yy", "p", "x", "swap"],
)

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="heavy-hex",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
    # Limit the number of optimization iterations to prevent the code cell from running
    # too long. Removing this line may improve results.
    options=dict(maxiter=1000),
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

Step 2: Optimize for quantum hardware execution

Next, we optimize the circuit for a target hardware. Typically, this step involves initializing the hardware backend and a pass manager for that backend. However, since the LUCJ ansatz is adapted to the hardware connectivity, we already performed these actions in the previous step. All that is left to do is run the pass manager on the circuit to transpile it to an ISA circuit that can be directly executed on the QPU.

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

Output:

Gate counts: OrderedDict({'xx_plus_yy': 86, 'p': 16, 'measure': 16, 'cp': 15, 'x': 10, 'swap': 2, 'barrier': 1})

Step 3: Execute using Qiskit primitives

After optimizing the circuit for hardware execution, we are ready to run it on the target hardware and collect samples for ground state energy estimation. As we only have one circuit, we will use Qiskit Runtime's Job execution mode and execute our circuit.

rng = np.random.default_rng()
sampler = StatevectorSampler(seed=rng)
job = sampler.run([isa_circuit], shots=100_000)

Output:

Warning: Trying to add QuantumRegister to a QuantumCircuit having a layout

primitive_result = job.result()
pub_result = primitive_result[0]

Step 4: Post-process and return result in desired classical format

A useful metric to judge the quality of the QPU output is the number of valid configurations returned. A valid configuration has the correct particle number and spin Z, which means that the right half of the bitstring has Hamming weight equal to the number of spin-up electrons, and the left half has Hamming weight equal to the number of spin-down electrons. The following cell computes the fraction of sampled configurations that are valid.

def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )


bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(f"Fraction of sampled configurations that are valid: {valid_fraction}")

Output:

Fraction of sampled configurations that are valid: 1.0

All of the bitstrings are valid because we are sampling the circuit on a noiseless simulator. When running on a noisy QPU, the fraction will be less than one, but hopefully it will be larger than the fraction you would expect if the bitstrings were sampled uniformly at random, which is computed in the following cell.

expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: {expected_fraction_random}"
)

Output:

Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625

Now, we estimate the ground state energy of the Hamiltonian using the diagonalize_fermionic_hamiltonian function. This function performs the self-consistent configuration recovery procedure to iteratively refine the noisy quantum samples to improve the energy estimate. We pass a callback function so that we can save the intermediate results for later analysis. See the API documentation for explanations of the arguments to diagonalize_fermionic_hamiltonian.

Here, we use the initial_occupancies argument to diagonalize_fermionic_hamiltonian to specify the Hartree-Fock configuration as the initial guess for the orbital occupancies in the ground state. This approach is sensible for systems where the ground state has significant support on the Hartree-Fock configuration, but it might not be appropriate in other situations, though more advanced computational methods might yield better initial guesses in those cases. Specifying initial_occupancies also allows configuration recovery to run even if no valid configurations were sampled, as may be the case when sampling a large circuit on a noisy QPU. Without this argument, the configuration recovery would fail and raise an error if no valid configurations were provided.

from functools import partial

from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
    solve_sci_batch,
)

# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 300
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the sci_solver argument
# in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)

# List to capture intermediate results
result_history = []


def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )


result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {energy_error}")

Output:

Iteration 1
	Subsample 0
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 1
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 2
		Energy: -108.59275573641656
		Subspace dimension: 900
Iteration 2
	Subsample 0
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 1
		Energy: -108.59275573641656
		Subspace dimension: 900
	Subsample 2
		Energy: -108.59275573641656
		Subspace dimension: 900
Final energy: -108.59275573641656
Final energy error: 0.0032316145694579745

Visualize the results

The first plot shows that in this simulation we are already within 1 mH of the exact answer after the first iteration (chemical accuracy is typically accepted to be 1 kcal/mol $\approx$ 1.6 mH). This is a small system, though, and because the samples are noiseless, configuration recovery is not needed. On a larger system run on a noisy QPU, multiple configuration recovery iterations might be needed, and the final accuracy might be worse. Generally, the energy can be improved by allowing more configuration recovery iterations or by increasing the number of samples per batch.

The second plot shows the average occupancy of each spatial orbital after the final iteration. We can see that both the spin-up and spin-down electrons occupy the first five orbitals with high probability in our solutions.

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-4)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

plt.tight_layout()
plt.show()

Output:

Large-scale hardware example

Now, we run a larger example on real quantum hardware. Here, we'll derive an active space for the nitrogen molecule from the cc-pVDZ basis set.

Steps 1-4

Here we put all the steps together into a single workflow at a larger scale, which is then run on real quantum hardware.

# ------------------------------ Step 1 ------------------------------
# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="cc-pvdz",
    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()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
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), norb)

# Store reference energy from SCI calculation performed separately
reference_energy = -109.22802921665716

print(f"norb = {norb}")
print(f"nelec = {nelec}")

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

# Set ansatz properties
n_reps = 1
pairs_aa = [(p, p + 1) for p in range(norb - 1)]
pairs_ab = None  # Let generate_lucj_pass_manager determine the alpha-beta interactions

# Initialize backend
service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=133
)
print(f"Using backend {backend.name}")

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="heavy-hex",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
    # Limit the number of optimization iterations to prevent the code cell from running
    # too long. Removing this line may improve results.
    options=dict(maxiter=1000),
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()


# ------------------------------ Step 2 ------------------------------

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")


# ------------------------------ Step 3 ------------------------------
sampler = Sampler(mode=backend)
sampler.options.environment.job_tags = ["TUT_SQD"]
job = sampler.run([isa_circuit], shots=100_000)
primitive_result = job.result()
pub_result = primitive_result[0]


# ------------------------------ Step 4 ------------------------------

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(f"Fraction of sampled configurations that are valid: {valid_fraction}")
expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: {expected_fraction_random}"
)
# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 300
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the sci_solver argument
# in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)

# List to capture intermediate results
result_history = []


result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    initial_occupancies=initial_occupancies,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {energy_error}")

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-4)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

plt.tight_layout()
plt.show()

Output:

converged SCF energy = -108.929838385609
norb = 26
nelec = (5, 5)
E(CCSD) = -109.2177884185544  E_corr = -0.2879500329450045
Using backend ibm_boston

Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].
Removing interaction (24, 24) from the end.
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].
Removing interaction (20, 20) from the end.

Gate counts: OrderedDict({'sx': 7039, 'rz': 6990, 'cz': 1858, 'x': 61, 'measure': 52, 'barrier': 1})
Fraction of sampled configurations that are valid: 0.02124
Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07
Iteration 1
	Subsample 0
		Energy: -109.13889134249762
		Subspace dimension: 120409
	Subsample 1
		Energy: -109.11785470455858
		Subspace dimension: 110889
	Subsample 2
		Energy: -109.13234360554011
		Subspace dimension: 130321
Iteration 2
	Subsample 0
		Energy: -109.16392179579177
		Subspace dimension: 223729
	Subsample 1
		Energy: -109.16281938332986
		Subspace dimension: 223729
	Subsample 2
		Energy: -109.16955816711932
		Subspace dimension: 233289
Iteration 3
	Subsample 0
		Energy: -109.17905772999075
		Subspace dimension: 324900
	Subsample 1
		Energy: -109.17532445048462
		Subspace dimension: 357604
	Subsample 2
		Energy: -109.1733168689756
		Subspace dimension: 348100
Iteration 4
	Subsample 0
		Energy: -109.18437778820451
		Subspace dimension: 474721
	Subsample 1
		Energy: -109.18450164209159
		Subspace dimension: 476100
	Subsample 2
		Energy: -109.18493571190754
		Subspace dimension: 487204
Iteration 5
	Subsample 0
		Energy: -109.18616522497996
		Subspace dimension: 622521
	Subsample 1
		Energy: -109.18652868888333
		Subspace dimension: 644809
	Subsample 2
		Energy: -109.18753326484406
		Subspace dimension: 585225
Final energy: -109.18753326484406
Final energy error: 0.040495951813099396

Next steps

Recommendations

If you found this work interesting, you might be interested in the following material:

Sample-based Krylov quantum diagonalization of a fermionic lattice model - a related tutorial using time evolution circuits instead of a variational ansatz
Scale SQD chemistry workflows with Dice solver - a page showing how to use the more efficient Dice software for diagonalization
SQD addon API documentation - reference for the diagonalize_fermionic_hamiltonian function
Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer - the paper this tutorial is based on

Was this page helpful?

Report a bug, typo, or request content on GitHub.