Quantendiagonalisierung eines Chemie-Hamiltonians auf der Grundlage von Stichproben

Geschätzte Nutzungsdauer: weniger als eine Minute auf einem Heron r2 Prozessor (HINWEIS: Dies ist nur eine Schätzung. Ihre Laufzeit kann variieren.)

Lernziele

• Wie man das SQD-Qiskit-Add-on verwendet, um die Grundzustandsenergie eines molekularen Systems anhand von Bitstrings zu approximieren, die von einer Quantenverarbeitungseinheit (QPU) abgetastet wurden.
• Wie man mit ffsim eine LUCJ-Schaltung (Local Unitary Cluster Jastrow) für die quantenchemische Simulation erstellt.

Voraussetzungen

• Quantenchemie und zweite Quantisierung
• Verwendung der Sampler-Primitive zum Abtasten von Quantenschaltungen

Hintergrund

In diesem Tutorial zeigen wir, wie man verrauschte Quantenstichproben nachbearbeitet, um den Grundzustand des Stickstoffmoleküls $\text{N}_2$ bei Gleichgewichtsbindungslänge zu approximieren, indem man das SQD-Qiskit-Add-on zur Implementierung des Algorithmus der stichprobenbasierten Quantendiagonalisierung (SQD) verwendet. Weitere Informationen zur Software finden Sie in der entsprechenden Dokumentation, einschließlich eines einfachen Beispiels für den Einstieg.

Die Stickstoffchemie verbindet dieses Beispiel mit dem Thema Nachhaltigkeit durch die Herstellung von Ammoniak für Düngemittel. Das Verständnis, wie Katalysatoren Stickstoff aktivieren, ist Teil der Forschung zu weniger energieintensiven Produktionsverfahren. Die Berechnung des Grundzustands eines isolierten Stickstoffmoleküls führt zu Methoden der Elektronenstruktur, die für diese umfassendere Herausforderung von Bedeutung sind.

Dieses Tutorial wird für Nutzer empfohlen, die mit Quantenchemie vertraut sind, insbesondere mit der Bestimmung der Grundzustandsenergien eines Moleküls. Eine ausführliche Anleitung zum Arbeitsablauf finden Sie im Kurs zum Quanten-Diagonalisierungsalgorithmus.

SQD ist eine Methode zur Bestimmung von Eigenwerten und Eigenvektoren von Quantenoperatoren, wie beispielsweise dem Hamilton-Operator eines Quantensystems, durch die kombinierte Nutzung von Quanten- und verteilter klassischer Datenverarbeitung. Im klassischen verteilten Rechnen werden die von einem Quantenprozessor gewonnenen Samples verarbeitet, um einen Ziel-Hamilton-Operator in einen von ihnen aufgespannten Unterraum zu projizieren und zu diagonalisieren. Ein auf SQD basierender Arbeitsablauf umfasst die folgenden Schritte:

1. Wählen Sie einen Schaltungsansatz und wenden Sie ihn auf einem Quantencomputer auf einen Referenzzustand (in diesem Fall den Hartree-Fock-Zustand ) an.
2. Beispiel-Bitstrings aus dem resultierenden Quantenzustand.
3. Führen Sie die selbstkonsistente Konfigurationswiederherstellung für die Bitstrings durch, um die Grundzustandsnäherung zu erhalten.

Es ist bekannt, dass SQD gut funktioniert, wenn der Ziel-Eigenzustand spärlich ist: Die Wellenfunktion ist in einer Menge von Basiszuständen $\mathcal{S} = \{|x\rangle \}$ enthalten, deren Größe nicht exponentiell mit der Größe des Problems zunimmt.

Quantenchemie

Der Hamiltonian eines molekularen Systems kann geschrieben werden als

wobei $h_{pr}$ und $h_{prqs}$ komplexe Zahlen sind, die als molekulare Integrale bezeichnet werden und mit Hilfe eines Computerprogramms aus der Spezifikation des Moleküls berechnet werden können. In diesem Tutorium berechnen wir die Integrale mit Hilfe des PySCF software-Paket.

Wie der molekulare Hamiltonian abgeleitet wird, ist in einem Lehrbuch der Quantenchemie nachzulesen (z. B. Modern Quantum Chemistry von Szabo und Ostlund). Eine ausführliche Erklärung, wie Probleme der Quantenchemie auf Quantencomputer abgebildet werden, finden Sie in der Vorlesung Mapping Problems to Qubits von der Qiskit Global Summer School 2024.

Lokaler unitäre Cluster Jastrow (LUCJ) Ansatz

SQD benötigt einen Quantenkreis-Ansatz, aus dem Stichproben entnommen werden können. In diesem Tutorial verwenden wir den Ansatz des lokalen unitären Clusters Jastrow (LUCJ), da er sowohl physikalisch fundiert als auch hardwarefreundlich ist. Wir werden ffsim verwenden, um die Ansatzschaltung zu erstellen.

Der LUCJ-Ansatz lässt sich an QPUs mit eingeschränkter Qubit-Konnektivität anpassen. Die Spinorbitale werden so auf Qubits abgebildet, dass der Ansatz keine Weiterleitung mit SWAP-Gattern erfordert. IBM® Die Hardware weist entweder eine Qubit-Topologie mit quadratischem Gitter oder mit dichtem hexagonalen Gitter auf. Für das Heavy-Hex-Gitter können wir das unten dargestellte „Zick-Zack“-Muster verwenden. In diesem Modell werden Orbitalen mit gleichem Spin auf Qubits mit einer Linientopologie abgebildet (rote und blaue Kreise), und alle vier Orbitalen besteht eine Verbindung zwischen Orbitalen mit unterschiedlichem Spin, wobei diese Verbindung durch ein Ancilla-Qubit (violette Kreise) ermöglicht wird.

Selbstkonsistente Konfigurationswiederherstellung

Das Verfahren zur Wiederherstellung der selbstkonsistenten Konfiguration ist so konzipiert, dass möglichst viele Signale aus verrauschten Quantenproben extrahiert werden. Da der molekulare Hamiltonian die Teilchenzahl und den Spin Z beibehält, ist es sinnvoll, einen Schaltkreis-Ansatz zu wählen, der diese Symmetrien ebenfalls beibehält. Bei Anwendung auf den Hartree-Fock-Zustand hat der resultierende Zustand eine feste Teilchenzahl und einen festen Spin Z in der geräuschlosen Einstellung. Daher sollten die Spin- $\alpha$ und Spin- $\beta$ Hälften jeder Bitkette, die aus diesem Zustand abgetastet wird, das gleiche Hamming-Gewicht wie im Hartree-Fock-Zustand haben. Aufgrund des Rauschens in aktuellen Quantenprozessoren werden einige gemessene Bitstrings diese Eigenschaft verletzen. Eine einfache Form der Nachselektion würde diese Bitstrings verwerfen, aber das ist verschwenderisch, weil diese Bitstrings noch ein Signal enthalten könnten. Das Verfahren der selbstkonsistenten Wiederherstellung versucht, einen Teil dieses Signals in der Nachbearbeitung wiederherzustellen. Das Verfahren ist iterativ und erfordert als Eingabe eine Schätzung der durchschnittlichen Besetzungen der einzelnen Orbitale im Grundzustand, die zunächst aus den Rohproben berechnet wird. Das Verfahren wird in einer Schleife ausgeführt, wobei jede Iteration die folgenden Schritte umfasst:

1. Für jeden Bitstring, der gegen die angegebenen Symmetrien verstößt, werden die Bits mit einem probabilistischen Verfahren umgedreht, um den Bitstring näher an die aktuelle Schätzung der durchschnittlichen Orbitalbesetzungen heranzuführen und so einen neuen Bitstring zu erhalten.
2. Sammeln Sie alle alten und neuen Bitstrings, die die Symmetrien erfüllen, und nehmen Sie eine Teilstichprobe mit einer festen Größe, die im Voraus festgelegt wird.
3. Projizieren Sie für jede Teilmenge von Bitstrings den Hamiltonian in den Unterraum, der von den entsprechenden Basisvektoren aufgespannt wird (eine Beschreibung dieser Basisvektoren finden Sie im vorherigen Abschnitt ), und berechnen Sie eine Grundzustandsschätzung des projizierten Hamiltonian auf einem klassischen Computer.
4. Aktualisieren Sie die Schätzung der durchschnittlichen Orbitalbesetzungen mit der Grundzustandsschätzung mit der niedrigsten Energie.

SQD-Workflow-Diagramm

Der SQD-Arbeitsablauf ist im folgenden Diagramm dargestellt:

Anforderungen

Vergewissern Sie sich vor Beginn dieses Tutorials, dass Sie die folgenden Komponenten installiert haben:

• Qiskit SDK v1.0 oder später, mit Visual isierungsunterstützung
• Qiskit Runtime v0.22 oder später (pip install qiskit-ibm-runtime)
• SQD-Qiskit-Add-on v0.12 oder höher (pip install qiskit-addon-sqd)
• ffsim v0.0.75 oder höher (pip install ffsim)

Einrichtung

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

Beispiel für einen Simulator im kleinen Maßstab

In diesem Tutorial werden wir eine Näherung für den Grundzustand eines Stickstoffmoleküls in der Nähe seines Gleichgewichtsbindungsabstands finden. Zunächst verwenden wir einen kleinen Basis-Satz von „ STO-6G “, um das Experiment zu simulieren und sicherzustellen, dass es funktioniert.

Schritt 1: Klassische Eingaben einem Quantenproblem zuordnen

Zunächst legen wir das Molekül und seine Eigenschaften fest.

# Specify molecule properties
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, -0.5)], ["N", (0, 0, 0.5)]],
    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)

Bevor wir die LUCJ-Ansatzschaltung konstruieren, führen wir zunächst eine CCSD-Berechnung in der folgenden Codezelle durch. Die $t_1$ und $t_2$ Amplituden aus dieser Berechnung werden zur Initialisierung der Parameter des Ansatzes verwendet.

# 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.1283731437052352

Nun verwenden wir ffsim, um das Ansatz-Schaltbild zu erstellen. Da unser Molekül einen Hartree-Fock-Zustand mit geschlossener Schale aufweist, verwenden wir die spinausgeglichene Variante des UCJ-Ansatzes, UCJOpSpinBalanced. Wir setzen optimize=True in der from_t_amplitudes Methode, um die „komprimierte“ Doppelfaktorisierung der Amplituden von $t_2$ s zu aktivieren (Einzelheiten finden Sie im Abschnitt „The local unitary cluster Jastrow (LUCJ) ansatz“ in der Dokumentation von ffsim).

Da sich der LUCJ-Ansatz an die verfügbare Konnektivität der QPU anpasst, müssen wir das QPU-Backend initialisieren, bevor wir den Ansatz erstellen. Zunächst erstellen wir ein generisches Backend mit einer quadratischen Gitterkopplungsabbildung und einem Gatter-Set, in das sich der LUCJ-Ansatz auf natürliche Weise zerlegen lässt. Anschließend verwenden wir, ffsim.qiskit.generate_lucj_pass_manager um einen Pass-Manager zu erstellen, der darauf spezialisiert ist, den LUCJ-Ansatz gemäß einem „Ladder“-Layout (siehe die Erläuterung in der Dokumentation von ffsim) in das angegebene Backend zu transpilieren. Diese Funktion nutzt eine Bewertungsheuristik, um die mit dem ausgewählten Layout verbundenen Fehler zu minimieren. Dies ist wichtig, wenn es sich bei Ihrem Backend um eine echte QPU oder um einen Simulator mit einem Rauschmodell handelt. Diese Funktion gibt nicht nur den Pass-Manager zurück, sondern auch die Alpha-Beta-Kopplungspaare, die auf der Hardware implementiert werden können. Wenn nicht alle Paare implementiert werden können, wird eine Warnung ausgegeben.

import warnings

from qiskit.transpiler import CouplingMap

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

# Set ansatz properties
n_reps = 4
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
coupling_map = CouplingMap.from_grid(norb, 2)
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="square",
    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,
)

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

Schritt 2: Optimierung für die Ausführung auf Quantenhardware

Als Nächstes optimieren wir die Schaltung für eine bestimmte Hardware. In der Regel umfasst dieser Schritt die Initialisierung des Hardware-Backends und eines Pass-Managers für dieses Backend. Da der LUCJ-Ansatz jedoch an die Hardware-Konnektivität angepasst ist, haben wir diese Schritte bereits im vorherigen Schritt durchgeführt. Nun muss nur noch der Pass-Manager auf der Schaltung ausgeführt werden, um sie in eine ISA-Schaltung zu transpilieren, die direkt auf der QPU ausgeführt werden kann.

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

Output:

Gate counts: OrderedDict({'xx_plus_yy': 254, 'cp': 88, 'p': 64, 'measure': 16, 'x': 10, 'barrier': 1})

Schritt 3: Mit „ Qiskit primitives “ ausführen

Nachdem wir die Schaltung für die Hardware-Ausführung optimiert haben, können wir sie nun auf der Zielhardware ausführen und Messwerte für die Schätzung der Grundzustandsenergie erfassen. Da wir nur einen Circuit haben, werden wir den Ausführungsmodus „ IBM Quantum “ des Compute Service verwenden und unseren Circuit ausführen.

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

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

Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format

Eine nützliche Kennzahl zur Beurteilung der Qualität der QPU-Ausgabe ist die Anzahl der zurückgegebenen gültigen Konfigurationen. Eine gültige Konfiguration hat die richtige Partikelanzahl und den richtigen Spin Z, was bedeutet, dass die rechte Hälfte der Bitfolge ein Hamming-Gewicht hat, das der Anzahl der Spin-up-Elektronen entspricht, und die linke Hälfte ein Hamming-Gewicht hat, das der Anzahl der Spin-down-Elektronen entspricht. Die folgende Zelle berechnet den Anteil der gültigen Konfigurationen aus der Stichprobe.

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

Alle Bitstrings sind gültig, da wir die Schaltung auf einem rauschfreien Simulator abtasten. Bei Ausführung auf einer verrauschten QPU wird dieser Anteil kleiner als eins sein, aber hoffentlich größer als der Anteil, den man erwarten würde, wenn die Bitstrings gleichmäßig zufällig ausgewählt würden – dieser Wert wird in der folgenden Zelle berechnet.

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: "
    f"{expected_fraction_random}"
)

Output:

Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625

Nun schätzen wir die Grundzustandsenergie des Hamiltonianers mit Hilfe der Funktion diagonalize_fermionic_hamiltonian . Diese Funktion führt das Verfahren der selbstkonsistenten Konfigurationswiederherstellung durch, um die verrauschten Quantenproben iterativ zu verfeinern und die Energieschätzung zu verbessern. Wir übergeben eine Callback-Funktion, damit wir die Zwischenergebnisse für eine spätere Analyse speichern können. In der API-Dokumentation finden Sie Erläuterungen zu den Argumenten für diagonalize_fermionic_hamiltonian.

Hier verwenden wir das initial_occupancies Argument, diagonalize_fermionic_hamiltonian um die Hartree-Fock-Konfiguration als anfängliche Schätzung für die Orbitalbelegungen im Grundzustand anzugeben. Dieser Ansatz ist sinnvoll für Systeme, bei denen der Grundzustand eine signifikante Unterstützung für die Hartree-Fock-Konfiguration aufweist, aber in anderen Situationen möglicherweise nicht geeignet ist, obwohl fortgeschrittenere Berechnungsmethoden in diesen Fällen zu besseren Anfangsschätzungen führen könnten. Durch die Angabe von kann die initial_occupancies Konfigurationswiederherstellung auch dann ausgeführt werden, wenn keine gültigen Konfigurationen erfasst wurden, wie dies beispielsweise bei der Erfassung eines großen Schaltkreises auf einer verrauschten QPU der Fall sein kann. Ohne dieses Argument würde die Konfigurationswiederherstellung fehlschlagen und einen Fehler auslösen, wenn keine gültigen Konfigurationen angegeben wurden.

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 = 1000
max_dim = 50
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,
    max_dim=max_dim,
    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.59589540253083
		Subspace dimension: 2500
	Subsample 1
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 2
		Energy: -108.59589540253083
		Subspace dimension: 2500
Iteration 2
	Subsample 0
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 1
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 2
		Energy: -108.59589540253083
		Subspace dimension: 2500
Final energy: -108.59589540253083
Final energy error: 9.19484552071026e-05

Visualisieren Sie die Ergebnisse

Die erste Grafik zeigt, dass wir in dieser Simulation bereits nach der ersten Iteration nahe 1 mH am exakten Ergebnis liegen (die chemische Genauigkeit wird üblicherweise mit 1 kcal/mol$\approx$ angegeben 1.6 mH). Es handelt sich hierbei jedoch um ein kleines System, und da die Stichproben rauschfrei sind, ist eine Konfigurationswiederherstellung nicht erforderlich. Bei einem größeren System, das auf einer rauschbehafteten QPU läuft, sind möglicherweise mehrere Iterationen zur Konfigurationswiederherstellung erforderlich, und die endgültige Genauigkeit könnte geringer ausfallen. Im Allgemeinen lässt sich die Energie verbessern, indem man mehr Iterationen zur Konfigurationswiederherstellung zulässt oder die Anzahl der Stichproben pro Stapel erhöht.

Das zweite Diagramm zeigt die durchschnittliche Belegung der einzelnen räumlichen Orbitale nach der letzten Iteration. Wir können sehen, dass sowohl die Spin-up- als auch die Spin-down-Elektronen die ersten fünf Orbitale mit hoher Wahrscheinlichkeit in unseren Lösungen besetzen.

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

# Chemical accuracy (+/- 1.6 milli-Hartree)
chem_accuracy = 0.0016

# 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_yscale("log")
axs[0].set_ylim(1e-5, 1e-2)
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:

Beispiel für groß angelegte Hardware

Nun führen wir ein umfangreicheres Beispiel auf echter Quantenhardware durch. Hier leiten wir einen aktiven Raum für das Stickstoffmolekül aus dem Basissatz „ cc-pVDZ “ ab. Außerdem werden wir folgende Änderungen vornehmen:

• Verringern Sie den Wert n_reps , 1 um Störgeräusche durch Gate-Fehler zu reduzieren
• Stellen Sie die Ansatz-Konnektivität entsprechend der Backend-Prozessorfamilie ein („heavy-hex“ für Heron, „square“ für Nighthawk)
• Erhöhen Sie auf max_dim , 1000 um größere Unterraumdimensionen zu ermöglichen.

Schritte 1–4

Hier fassen wir alle Schritte zu einem einzigen, groß angelegten Arbeitsablauf zusammen, der dann auf echter Quantenhardware ausgeführt wird.

# ------------------------------ Step 1 ------------------------------
# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, -0.5)], ["N", (0, 0, 0.5)]],
    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)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

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

# Create pass manager
connectivity_map = {"Heron": "heavy-hex", "Nighthawk": "square"}
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity=connectivity_map[
        backend.configuration().processor_type["family"]
    ],
    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,
)

# 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: "
    f"{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 = 1000
max_dim = 1000
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,
    max_dim=max_dim,
    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]

# Chemical accuracy (+/- 1.6 milli-Hartree)
chem_accuracy = 0.0016

# 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_yscale("log")
axs[0].set_ylim(1e-3, 1e-1)
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.2879500329450041
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': 7036, 'rz': 6965, 'cz': 1858, 'measure': 52, 'x': 50, 'barrier': 1})
Fraction of sampled configurations that are valid: 0.01677
Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07
Iteration 1
	Subsample 0
		Energy: -109.17922180754668
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.17713541195594
		Subspace dimension: 962361
	Subsample 2
		Energy: -109.17407767345108
		Subspace dimension: 968256
Iteration 2
	Subsample 0
		Energy: -109.1873108379026
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.1876970140745
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.18255107152576
		Subspace dimension: 1000000
Iteration 3
	Subsample 0
		Energy: -109.19052148187927
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19316375313544
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19087348337145
		Subspace dimension: 1000000
Iteration 4
	Subsample 0
		Energy: -109.19374847699348
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19517039328333
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19422416032815
		Subspace dimension: 1000000
Iteration 5
	Subsample 0
		Energy: -109.19580768221108
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19608553581007
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19548131833729
		Subspace dimension: 1000000
Final energy: -109.19608553581007
Final energy error: 0.03194368084709254

Nächste Schritte