SqDRIFT-Algorithmus zur Grundzustandsabschätzung

Nutzungsschätzung: 180 Sekunden auf einem Heron-r3-Prozessor (HINWEIS: Dies ist nur eine Schätzung. Deine Laufzeit kann abweichen.)

Lernziele

• Lerne, wie du Circuits mit geringerer Tiefe im Vergleich zur Trotterisierung erstellst

• Durchlaufe einen End-to-End-Workflow zur Grundzustandsabschätzung mit qDRIFT und SQD

• Lerne, wie du qiskit-fermions zusammen mit anderen Qiskit-Addons einsetzt, um einen solchen Workflow umzusetzen

Dieses Tutorial wird zu Lehrzwecken als Python-Notebook präsentiert.

Voraussetzungen

• Lies die Übersicht zur probenbasierten Quantendiagonalisierung (SQD)

• Lies die Lektion zur probenbasierten Krylov-Quantendiagonalisierung (SKQD)

Hintergrund

SqDRIFT ist eine Variante von SKQD, die die Notwendigkeit, einen Ansatz zum Sampeln von Bitstrings zu wählen, durch ein Ensemble von Zeitentwicklungs-Circuits ersetzt, die direkt aus dem Ziel-Hamiltonoperator konstruiert werden. Dies wird erreicht, indem kleinere Zeitentwicklungsoperatoren anhand ihrer Koeffizienten aus dem Hamiltonoperator gesampelt werden, was als qDRIFT-Trotterisierungsmethode bekannt ist.

Dieses Tutorial verwendet Qiskit Fermions, um die natürlicheren fermionischen Circuits für den qDRIFT-Algorithmus zu erstellen, gefolgt von fermionischen Layout- und Synthese-Passes, bevor die Circuits in die traditionelle Qiskit-Pipeline zur Hardware-Ausführung eingespeist werden.

Der Hamiltonoperator habe die Form:

$$H = \sum_{i=1}^{N} c_i h_i$$

wobei wir ohne Beschränkung der Allgemeinheit $c_i > 0$ verlangen und dass der größte Eigenwert von $h_i$ betragsmäßig gleich $1$ ist. Jeder vorzeichenbehaftete oder komplexe Vorfaktor wird in $h_i$ absorbiert, sodass die Koeffizienten $c_i$ streng positive Gewichte sind, während die $h_i$ die Richtung jedes Terms tragen. Hier ist $N$ die Anzahl der Terme (oder, nach der Gruppierung, die Anzahl der Gruppen) im Hamiltonoperator; sie ist eine Eigenschaft des Hamiltonoperators und von der Anzahl der in einen einzelnen Circuit gesampelten Operatoren zu unterscheiden, die unten mit $n$ bezeichnet wird.

Der qDRIFT-Algorithmus realisiert dann für die Zielzeit $t$ einen Operator $V_k$, wobei $k$ von $1 \cdots K$ läuft und den $k_{th}$ SqDRIFT-Circuit bezeichnet, definiert als:

$$V_k = \prod_{j=1}^{n} e^{-i h_{k_j} \lambda t / n }$$

Hier ist $n$ die Anzahl der gesampelten Operatoren pro Circuit und $K$ die Anzahl der Circuits im Ensemble. Das Produkt läuft über die $n$ Ziehungen, nicht über alle $N$ Terme des Hamiltonoperators, und da die Terme mit Zurücklegen gezogen werden, kann dasselbe $h_i$ in einem einzelnen $V_k$ mehr als einmal vorkommen.

Die Größe:

$$\lambda = \sum_{i=1}^{N} c_i$$

ist die $L_1$-Norm der Koeffizienten, sodass jeder der $n$ Schritte unabhängig davon, welcher Term gezogen wurde, über dieselbe Dauer $\lambda t / n$ entwickelt. Die Einheitlichkeit des Schrittwinkels ist das charakteristische Merkmal von qDRIFT: Ein Koeffizient beeinflusst das Ergebnis dadurch, wie oft sein Term gezogen wird, nicht dadurch, wie weit dieser Term rotiert wird. Die Indizes werden aus der folgenden Verteilung gesampelt:

$$P[k_i] = \frac{c_i}{\lambda}$$

sodass die Folge $(k_1, \ldots, k_n)$ eine zufällige Folge von Termindizes ist, die aus dieser Verteilung gezogen wird. Da die $c_i$ positiv sind und sich zu $\lambda$ summieren, ist dies eine normierte Wahrscheinlichkeitsverteilung, und der Erwartungswert des resultierenden Kanals über die zufälligen Ziehungen approximiert die Entwicklung unter $H$, mit einem Fehler, der mit wachsendem $n$ abnimmt. Beachte, dass der Approximationsfehler von $\lambda$ abhängt und nicht von der Anzahl der Terme $N$.

(Das SqDRIFT-Paper schreibt die Anzahl der Terme als $\mathcal{N}$ und die Sequenzlänge als $N$; wir verwenden hier $N$ und $n$, um die beiden klar auseinanderzuhalten.)

Dieses Tutorial zeigt, wie man ein Ensemble solcher randomisierten Circuits erzeugt. Nachdem wir diese Circuits erstellt haben, sampeln wir ähnlich wie beim Aufbau eines Krylov-Unterraums für verschiedene Operatoren Bitstrings aus mehreren solchen Operatoren mit unterschiedlichen Zeitparametern. Das sorgt für eine höhere Überlappung zwischen den Grundzustandsvektoren und den gesampelten Bitstrings.

Anforderungen

Stelle vor dem Start dieses Tutorials sicher, dass du Folgendes installiert hast

• Eine virtuelle Python-Umgebung (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (Beachte, dass der Name im Plural steht)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Du kannst alle erforderlichen Pakete installieren mit:

pip install "qiskit~=2.5" "qiskit-fermions==0.1.0" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy

Setup

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "qiskit_aer": "qiskit-aer", "qiskit_fermions": "qiskit-fermions", "qiskit_ibm_runtime": "qiskit-ibm-runtime"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

# Qiskit core
from qiskit import transpile
from qiskit.primitives import BitArray

# Qiskit Aer
from qiskit_aer import AerSimulator

# IBM Quantum Compute Service
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler

# Qiskit Fermions
from qiskit_fermions.operators.library import FCIDump
from qiskit_fermions.operators import FermionOperator
from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms
from qiskit_fermions.operators.terms.grouping import (
    group_terms_by_electronic_structure,
)
from qiskit_fermions.operators.terms.ordering import canonical_order
from qiskit_fermions.circuit import FermionicCircuit
from qiskit_fermions.circuit.library import Evolution
from qiskit_fermions.transpiler import FermionicPassManager
from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
from qiskit_fermions.transpiler.passes import QDriftTrotterization
from qiskit_fermions.circuit.library import InitializeModes

# Qiskit addon SQD
from qiskit_addon_sqd.fermion import (
    diagonalize_fermionic_hamiltonian,
    SCIResult,
)

Simulatorbeispiel

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Den FCIDump lesen und vorbereiten

Für dieses Tutorial laden wir den Elektronenstruktur-Hamiltonoperator für Stickstoff (N2). Es gibt auch andere Möglichkeiten, fermionische Operatoren zu erstellen. Siehe die Dokumentation unter qiskit_fermions.operators.library.

Über diesen FCIDump. Die Datei N2_sto_3g beschreibt ein Stickstoffmolekül ($N_2$) in der minimalen STO-3G-Basis bei einem interatomaren Abstand von 1.09 $\AA$, der experimentellen Gleichgewichtsbindungslänge. Ihr Header deklariert NORB=10, NELEC=14 und MS2=0: 10 räumliche Orbitale (also 20 Spinorbitale und 20 Qubits unter Jordan-Wigner), 14 Elektronen in einem Spin-Singulett, also sieben $\alpha$- und sieben $\beta$-Elektronen. Allen Orbitalen ist das Symmetrielabel 1 zugewiesen, das heißt, es wird keine Punktgruppensymmetrie ausgenutzt. Da es sich um einen STO-3G-Dump des vollen Raums handelt, sind keine Orbitale eingefroren und der Korrelationsraum ist klein genug, dass eine exakte FCI-Referenzenergie zum Vergleich klassisch berechnet werden kann, wie in der nächsten Zelle gezeigt.

Eine äquivalente Datei kann mit PySCF neu erzeugt werden:

from pyscf import gto, scf, tools

mol = gto.M(atom="N 0 0 0; N 0 0 1.09", basis="sto-3g", symmetry=False)
mf = scf.RHF(mol).run()
tools.fcidump.from_scf(mf, "N2_sto_3g")

Da die Integrale von den konvergierten SCF-Orbitalen abhängen, kann sich eine neu erzeugte Datei von der mitgelieferten in Orbitalphase oder -reihenfolge unterscheiden; die Gesamtenergien sind davon nicht betroffen.

Die Datei beschaffen. Den FCIDump findest du in diesem GitHub-Repository. Du kannst die Zelle unten ausführen, um ihn an den Ort zu laden, den der Rest des Tutorials erwartet.

Zuerst verwenden wir den von pyscf bereitgestellten cisolver, um die Referenzenergie zu erhalten. Das ist die wahre Grundzustandsenergie des Moleküls, mit dem wir arbeiten. Dazu deklarieren wir zunächst norb und nelec, die Anzahl der Orbitale bzw. die Anzahl der Elektronen. Dann deklarieren wir h1e und h2e, die Ein- bzw. Zweielektronenintegrale. Alle diese werden später auch für SQD verwendet.

import os
from urllib.request import urlopen

# The FCIDump is stored with this tutorial in the Qiskit documentation repository.
FCIDUMP_URL = "https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g"
FCIDUMP_PATH = "fcidump_files/N2_sto_3g"

if not os.path.exists(FCIDUMP_PATH):
    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)
    with urlopen(FCIDUMP_URL) as response:
        contents = response.read()
    with open(FCIDUMP_PATH, "wb") as f:
        f.write(contents)
    print(f"Downloaded FCIDump to {FCIDUMP_PATH}")
else:
    print(f"Using existing FCIDump at {FCIDUMP_PATH}")

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

Den Hamiltonoperator laden

Wenn die nötigen Daten bereitstehen, lesen wir den Hamiltonoperator aus der FCI-Datei in einem Format ein, das mit qiskit-fermions kompatibel ist

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

Fermionische Workflows mit qiskit-fermions

Wir bilden den Hamiltonoperator zunächst mit qiskit-fermions auf ein fermionisches Circuit-Modell ab, das Transpiler-Passes und Gates speziell für fermionische Circuits bereitstellt. Diese werden später in diesem Workflow vor den traditionellen Transpiler-Passes von Qiskit verwendet.

Termgruppierung

Um die Reproduzierbarkeit der Ergebnisse sicherzustellen, verwenden wir zuerst canonical_order, um die Terme allein nach ihrer Struktur zu sortieren. Die Reihenfolge der Operatoren in der Liste canon ist damit festgelegt. Das gewährleistet die Reproduzierbarkeit der erstellten Operatoren, weil der Pass QDriftTrotterization, den wir später verwenden werden, zufällige Indizes zur Erstellung der qDRIFT-Operatoren sampelt.

In diesem Schritt nutzen wir die vielen Symmetrien aus, die im Elektronenstruktur-Hamiltonoperator vorhanden sind, indem wir verwandte Terme mit identischen Koeffizienten gruppieren. Dies verändert zwar die Verteilung der Operatorkoeffizienten, aus der das qDRIFT-Protokoll sampelt, beeinträchtigt aber nicht seine Konvergenzgarantien. Entscheidend ist, dass das Gruppieren symmetrieverwandter Terme zu einer günstigen Auslöschung von Pauli-Termen und insgesamt zu einer geringeren Circuit-Tiefe führt, wenn ein Zustand unter ihrer Wirkung zeitentwickelt wird.

qiskit-fermions stellt die Funktion group_terms_by_electronic_structure bereit, die diese Gruppierung für uns übernimmt.

Beachte, dass group_terms_by_electronic_structure normal geordnete Terme voraussetzt.

Diagonalterme filtern

Wir entfernen die Diagonalterme aus dem Hamiltonoperator, der zur Erzeugung der Circuits verwendet wird, damit die $n$ qDRIFT-Sampling-Slots für Terme aufgewendet werden, die Population zwischen Konfigurationen verschieben. Solche Terme filtert man am besten an dieser Stelle aus dem Hamiltonoperator heraus, bevor im nächsten Schritt das Evolution-Gate konstruiert wird.

Die fraglichen Terme sind diejenigen, die in der Besetzungszahlbasis diagonal sind, das heißt die Produkte von Anzahloperatoren $a^\dagger_i a_i$. Drei Arten von Termen fallen unter diese Beschreibung:

• der konstante Energieoffset, ein Produkt aus null Anzahloperatoren, dessen Zeitentwicklung nur eine globale Phase beiträgt;

• die einzelnen Anzahloperatoren $n_i$, deren Zeitentwicklung sich auf $Z$-Rotationen einzelner Qubits reduziert;

• die Produkte höherer Ordnung wie $n_i n_j$.

Für sich genommen verschiebt keiner von ihnen Population zwischen Besetzungszahlkonfigurationen; sie wirken nur auf die Phasen der bereits vorhandenen Konfigurationen. Sie sind jedoch nicht wirkungslos: Diese relativen Phasen fließen in die Interferenz ein, die von den Anregungstermen später im Circuit erzeugt wird, sodass ihr Herausfiltern die tatsächlich erzeugte Entwicklung verändert und die Sampling-Verteilung verändern kann. Dies ist eine bewusste Näherung im Schritt der Circuit-Erzeugung, die das Sampling auf Anregungsterme fokussieren soll, und kein Schritt, der die gesampelte Verteilung unberührt lässt. Anders als die obige Symmetriegruppierung, die die Konvergenzgarantien von qDRIFT intakt lässt, verändert dieser Filter den entwickelten Operator. Die Circuits approximieren daher nicht mehr die Entwicklung unter dem vollen Hamiltonoperator, und die qDRIFT-Fehlerschranken gelten für den gefilterten Operator statt für den ursprünglichen. Das ist hier akzeptabel, weil die Circuits nur eine Sampling-Heuristik zum Vorschlagen von Konfigurationen sind: Kein Term geht in der Energieabschätzung selbst verloren, da der Filter nur auf den Hamiltonoperator angewendet wird, der zum Bau der Circuits dient, während die spätere klassische Diagonalisierung den vollen Hamiltonoperator einschließlich der Diagonalterme verwendet. Die Genauigkeit von SQD hängt von diesem klassischen Schritt ab, der unabhängig davon, wie die Konfigurationen vorgeschlagen wurden, im gesampelten Unterraum variationell bleibt.

Die Funktion filter_diagonal_terms() entfernt solche Terme in place aus einem Operator. Sie erkennt sie an ihrer normal geordneten Struktur — die Multimenge der Erzeugungsmoden stimmt mit der Multimenge der Vernichtungsmoden überein — und ist daher nur für einen Operator gültig, der bereits normal geordnet ist. Diese Annahme wird zur Laufzeit nicht geprüft.

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

5060

Nachdem wir die Terme im Hamiltonoperator gruppiert haben, legen wir die folgenden Parameter zur Erzeugung des Circuit-Ensembles fest:

• Die Anzahl der zu erzeugenden Circuits: num_circuits
• Die Länge jedes Circuits in Anregungsgruppen: num_exc
• Der Faktor für die verschiedenen Evolutionszeiten: times

Fermionische Circuits erstellen

Wir erstellen nun fermionische Circuits für jeden der Zeitschritte. Jeder Circuit besteht aus einem einzelnen Evolutions-Gate mit der Evolutionszeit, die wir zuvor deklariert haben. Der Evolutionsoperator ist der Hamiltonoperator. Später führen wir Transpiler-Passes auf diesen Circuits aus, um qDRIFT-Circuits zu erzeugen.

Ansatz-Vorbereitung

Wir präparieren den Hartree-Fock-Zustand mit der Klasse InitializeModes. Für Stickstoff besteht der Vorgang einfach darin, X-Gates auf die ersten num_elec_a Qubits und dann auf die num_elec_b Qubits anzuwenden, die beide für Stickstoff gleich sieben sind. Dieser Zustand repräsentiert die sieben $\alpha$- und sieben $\beta$-Elektronen von Stickstoff.

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []

hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

Schritt 2: Das Problem für die Ausführung auf Quantenhardware optimieren

Nachdem wir nun unsere Circuits haben, verwenden wir zunächst die in qiskit-fermions verfügbaren Passes, um Optimierungen auf fermionischer Ebene durchzuführen, und transpilieren unseren Circuit anschließend für das gewählte Backend. Da dies ein Simulator-Experiment ist, tun wir dies zunächst für den AerSimulator. Gewichtsberechnung für jede Gruppe

In diesem Schritt führen wir das qDRIFT-Sampling der Terme stochastisch durch, mit Wahrscheinlichkeiten proportional zu ihren Koeffizienten im Hamiltonian. Der qDRIFT-Transpiler-Pass übernimmt das für uns. Wir können nun flachere Circuits erstellen, die sich trotz begrenzter Qubit-Konnektivität effizienter auf der Hardware ausführen lassen, selbst wenn der Hamiltonian langreichweitige Kopplungen und Terme höherer als quadratischer Ordnung enthält. Nach der Termgruppierung werden die Operatoren anhand ihrer Gewichte gesampelt. Für jeden Operator $h_i$ ist das Gewicht $W_{h_i}$ wie folgt definiert:

$$W_{h_i} = |c_i| / \lambda$$

Da die Terme in Schritt 1 gruppiert wurden, ist jedes $h_i$ hier eine ganze Gruppe: $c_i$ ist der mittlere Absolutwert der Koeffizienten der Terme in Gruppe $i$, und jeder Term der Gruppe wird mit seinem auf das Vorzeichen reduzierten Koeffizienten zeitentwickelt.

Fermionische und hardwarenative Optimierungen

Die Funktion generate_preset_jw_pass_manager() gibt einen MultiStagePassManager zurück, der einen FermionicCircuit entgegennimmt und einen optimierten finalen Circuit erzeugt, den wir für die Ausführung auf unserer Hardware transpilieren können. Wir ersetzen seine standardmäßige Optimierungsstufe durch einen FermionicPassManager, der unseren QDriftTrotterization-Pass enthält:

• Der QDriftTrotterization-Pass nutzt intern die Gewichtsberechnung und das Sampling, um die Circuits zu erzeugen, die wir für das Sampling verwenden werden

• Der RelabelModes-Pass ist ein weiterer Optimierungs-Pass, mit dem sich die fermionischen Moden permutieren lassen, um die Konnektivität über die Qubits hinweg zu optimieren und die Gate-Tiefe zu verringern; mehr dazu in der API-Referenz

Die übrigen Stufen des MultiStagePassManager laufen automatisch ab und übernehmen die vollständige Fermion-zu-Qubit-Abbildung:

• F2QLayout: Der voreingestellte Pass Manager wendet den TrivialF2QLayout-Pass an, der $n$ fermionische Bits trivial auf $n$ Qubits abbildet.

• F2QSynth: Ein Transpilations-Pass, der fermionbasierte Circuit-Anweisungen auf qubitbasierte abbildet.

qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

400

Nachdem die Optimierungen auf fermionischer Ebene abgeschlossen sind, können wir die Circuits für die Ausführung auf dem Simulator transpilieren.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Schritt 3: Mit Qiskit-Primitiven ausführen

Nachdem wir nun unsere Circuits haben, können wir sie mit Qiskit-Primitiven auf dem AerSimulator ausführen. Wir fassen alle Counts aus den verschiedenen Circuits zusammen. Wir wandeln sie in boolesche Vektoren um, bevor wir sie abschließend mit SQD nachverarbeiten.

print(
    f"Executing {len(transpiled_circuits)} circuits with {shots} shots each..."
)

job = simulator.run(transpiled_circuits, shots=shots)
result = job.result()

all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]

print(len(all_counts), "length before post processing")

Executing 400 circuits with 100 shots each...
400 length before post processing

Schritt 4: Nachverarbeiten und Ergebnis im gewünschten klassischen Format zurückgeben

Bitstrings für SQD verwenden

Wir können nun das Diagonalisierungsverfahren auf die ausgewählten Bitstrings anwenden, um den niedrigsten Eigenwert zu finden, der der Grundzustandsenergie des Moleküls entspricht. Wir erstellen eine Callback-Funktion, legen die anfänglichen Besetzungen fest und setzen die Parameter, bevor wir das Diagonalisierungsverfahren schließlich ausführen. Die Callback-Funktion wird verwendet, um bei jeder Iteration die aktuelle Iteration und die aktuelle Eigenwertschätzung auszugeben.

Um schließlich die Schätzung des Grundzustands zu erhalten, addieren wir die nuclear_repulsion_energy zur resultierenden Energie.

Hinweis: Die Dimension des Unterraums ist über die Iterationen hinweg nicht fest, selbst auf dem rauschfreien Simulator — jedes Subsample zieht eine andere Menge an Konfigurationen, und der Recovery-Schritt formt den Pool zwischen den Iterationen um, sodass die gemeldete Dimension von einem Subsample zum nächsten variiert. Rauschfreies Sampling legt die Dimension des ausgewählten Unterraums also nicht von selbst fest. Der Hardwarelauf liefert dagegen tendenziell systematisch größere Unterräume, weil verrauschte Shots die Teilchenzahlsymmetrie brechen und die Konfigurationswiederherstellung sie in zusätzliche Basisvektoren verwandelt. Deshalb führen wir im Hardware-Abschnitt außerdem einen weiteren Schritt zum Beschneiden von Bitstrings ein.

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

40000
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64767025226178
		Subspace dimension: 5538
	Subsample 1
		Energy: -107.64772799119115
		Subspace dimension: 5670
	Subsample 2
		Energy: -107.64765512281548
		Subspace dimension: 5767
Iteration 2
	Subsample 0
		Energy: -107.64795948524682
		Subspace dimension: 6080
	Subsample 1
		Energy: -107.64806617355072
		Subspace dimension: 6300
	Subsample 2
		Energy: -107.64802260640258
		Subspace dimension: 6308
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712
 0.99301927 0.0373282  0.0373266  0.00944508]
Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737
 0.99302145 0.03733536 0.03733399 0.0094585 ]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6480661736 Ha
Error:            1.1811817564e-04 Ha

Hardwarebeispiel

Dieses Beispiel verwendet 20 Qubits (10 räumliche Orbitale). Diese Wahl ist eine Vereinfachung für ein Tutorial, das schnell laufen soll, und keine harte Obergrenze der Methode.

Die Kosten des klassischen Schritts werden nicht direkt von der Qubit-Anzahl bestimmt. SQD diagonalisiert den Hamiltonian, projiziert auf den Unterraum, der von den gesampelten Konfigurationen aufgespannt wird. Die klassischen Kosten werden also von der Dimension dieses ausgewählten Unterraums getrieben — hier bestimmt durch samples_per_batch, num_batches und die Anzahl der verschiedenen Konfigurationen, die die Circuits tatsächlich erzeugen — zusammen mit der dünnbesetzten linearen Algebra, die nötig ist, um den projizierten Hamiltonian anzuwenden. Der vollständige CI-Raum wächst kombinatorisch mit Orbitalen und Elektronen, aber der ausgewählte Unterraum ist ein kleiner, einstellbarer Ausschnitt davon, dessen Größe wir direkt steuern. Folglich lassen sich die Anzahl der Qubits und die klassische Schwierigkeit einigermaßen unabhängig voneinander variieren: Ein breiterer Orbitalraum, der in einen bescheidenen Unterraum gesampelt wird, kann günstiger sein als ein kleineres System, das über einen sehr großen Unterraum diagonalisiert wird.

In der Praxis hängt die machbare Systemgröße also von der Unterraumdimension ab, die du für die gewünschte Genauigkeit benötigst, sowie vom Speicher und den Kernen, die dem Eigenwertlöser zur Verfügung stehen. Größere Orbitalräume erfordern typischerweise tatsächlich einen größeren Unterraum, um chemische Genauigkeit zu erreichen, und genau das motiviert letztlich verteilte Ressourcen — siehe qiskit-addon-sqd-hpc, um diesen Schritt zu skalieren. Statt eine feste Grenze anzunehmen, besteht der praktische Ansatz darin, die gemeldete Unterraumdimension und die Energiekonvergenz über die Iterationen hinweg zu beobachten und die Unterraumgröße zu erhöhen, bis sich die Energie nicht mehr verbessert oder der verfügbare Speicher erschöpft ist.

Hinweis: Aufgrund von Sampling-Fehlern durch das Hardwarerauschen ist der für die Diagonalisierung erzeugte Unterraum im Hardwarelauf größer als der, den wir mit dem Simulator erhalten. Das vergrößert zwar die Dimension des Unterraums, den wir diagonalisieren wollen, doch der Workflow liefert dank der Robustheit von SQD gegenüber Rauschen weiterhin eine genaue Antwort.

Beschneiden fehlerhafter Strings

Hier können wir einen zusätzlichen Schritt durchführen. Wenn wir alle Bitstrings aus den Circuit-Ausführungen haben, können wir entweder die ungültigen Bitstrings herausfiltern, bevor wir SQD ausführen, oder ohne Beschneiden fortfahren. Auf das Beschneiden zu verzichten ist bei Hardwarelaufen im Allgemeinen vorzuziehen, weil die symmetriegebrochenen Shots dann für die Konfigurationswiederherstellung erhalten bleiben, die sie in gültige Konfigurationen überführen und so den Unterraum erweitern kann, statt diese Shots einfach zu verwerfen.

Da Stickstoff nur sieben $\alpha$- und sieben $\beta$-Elektronen haben kann, können alle Bitstrings verworfen werden, die in der ersten und zweiten Hälfte der Ausgabe mehr oder weniger als sieben Einsen enthalten. Wir definieren eine Funktion, die prüft, ob die Bitstrings gültig sind, und sie andernfalls verwirft. Sobald wir die fehlerhaften Bitstrings herausgefiltert haben, werden die übrigen an das Diagonalisierungsverfahren übergeben. Verwende das Flag PRUNE unten, um zwischen den beiden Verhaltensweisen umzuschalten.

Beachte, dass das Beschneiden nur eine von mehreren Entscheidungen ist, die den endgültigen Unterraum prägen, neben der Anzahl der Circuits, der Menge der Evolutionszeiten und dem Filtern diagonaler Terme. Ein beschnittener Lauf lässt sich nur dann sinnvoll mit einem unbeschnittenen vergleichen, wenn alles andere konstant gehalten wird; die C++-Version von diesem Tutorial behandelt das ausführlicher, da sie postselektiert statt wiederherzustellen und sich auch in diesen anderen Parametern unterscheidet.

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []
hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

# Calculate weights for sampling (one per group)
qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService(channel="ibm_quantum_platform")

# Select backend (choose based on qubit requirements)
backend = service.least_busy(
    operational=True,
    simulator=False,
    min_num_qubits=2 * norb,
)

print(f"Selected backend: {backend.name} ({backend.num_qubits} qubits)")

# Transpile for hardware
transpiled_circuits = transpile(
    sqdrift_circuits,
    backend=backend,
    optimization_level=3,
    seed_transpiler=42,
)

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-SqDRIFT"]

job = sampler.run(transpiled_circuits, shots=shots)
result = job.result()

# Extract counts from SamplerV2 results
all_counts = [pub_result.data.meas.get_counts() for pub_result in result]

# Set to True to filter out bitstrings that violate electron-number conservation
PRUNE = False

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
    )

if PRUNE:
    all_counts_filtered = []
    for counts in all_counts:
        filtered_count = {}
        for key in counts:
            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):
                continue
            elif key not in filtered_count.keys():
                filtered_count[key] = counts[key]
            else:
                filtered_count[key] += counts[key]
        all_counts_filtered.append(filtered_count)
    all_counts = all_counts_filtered

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print("Electron configuration:")
print(f"  Total electrons: {nelec}")
print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

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

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha
5060
400

Selected backend: ibm_aachen (156 qubits)
40000
Electron configuration:
  Total electrons: 14
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64593072647523
		Subspace dimension: 7221
	Subsample 1
		Energy: -107.6458270048177
		Subspace dimension: 7209
	Subsample 2
		Energy: -107.64007673117075
		Subspace dimension: 7138
Iteration 2
	Subsample 0
		Energy: -107.64757372124944
		Subspace dimension: 9009
	Subsample 1
		Energy: -107.64674060104392
		Subspace dimension: 8245
	Subsample 2
		Energy: -107.64731360491942
		Subspace dimension: 8178
Iteration 3
	Subsample 0
		Energy: -107.64765518770588
		Subspace dimension: 8835
	Subsample 1
		Energy: -107.64767975712016
		Subspace dimension: 8649
	Subsample 2
		Energy: -107.64761634415606
		Subspace dimension: 8648
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295
 0.99298797 0.03728154 0.03728186 0.00938359]
Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174
 0.99298424 0.03728121 0.03728169 0.00939159]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6476797571 Ha
Error:            5.0453460619e-04 Ha

Nächste Schritte

Empfehlungen

Wenn du diese Arbeit interessant fandest, interessieren dich vielleicht die folgenden Materialien:

• Sample-based Krylov quantum diagonalization of a fermionic lattice model – ein verwandtes Tutorial, das Zeitentwicklungs-Circuits statt eines variationellen Ansatzes verwendet.
• Sample-based quantum diagonalization of a chemistry Hamiltonian – ein Tutorial dazu, wie man einen Local-Unitary-Cluster-Jastrow-Circuit (LUCJ) für die Quantenchemie-Simulation konstruiert.
• Das Paper SqDRIFT – die Literatur, auf der dieses Tutorial basiert. (Beachte, dass einige der in diesem Paper besprochenen Optimierungen derzeit in Arbeit sind und dieses Tutorial sich künftig je nach Weiterentwicklung der verwendeten Bibliotheken ändern kann.)