SqDRIFT-Algorithmus zur Grundzustandsabschätzung

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

Suchst du die C++-Version?

Dieses Tutorial verwendet Python. Für die C++-Implementierung, einschließlich Quellcode und Build-Anweisungen, siehe das C++-SqDRIFT-Tutorial.

Lernziele

• Lerne, wie man Circuits mit geringerer Tiefe im Vergleich zur Trotterisierung erstellt

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

• Lerne, wie man qiskit-fermions zusammen mit anderen Qiskit-Addons verwendet, um einen solchen Workflow zu implementieren

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

Voraussetzungen

• Lies den Überblick zu Sample-based quantum diagonalization (SQD)

• Lies die Lektion zu Sample-based Krylov Quantum Diagonalization (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-Hamiltonian konstruiert werden. Dies wird erreicht, indem kleinere Zeitentwicklungsoperatoren basierend auf ihren Koeffizienten aus dem Hamiltonian subgesampelt werden, was als qDRIFT-Trotterisierungsmethode bekannt ist.

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

Sei der Hamiltonian von der Form:

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

wobei wir ohne Beschränkung der Allgemeinheit $c_i > 0$ fordern und dass der größte Eigenwert von $h_i$ im Betrag gleich $1$ ist. Jeder vorzeichenbehaftete oder komplexe Vorfaktor wird in $h_i$ absorbiert, sodass die Koeffizienten $c_i$ strikt positive Gewichte sind, während die $h_i$ die Richtung jedes Terms tragen. Hier ist $N$ die Anzahl der Terme (oder, nach Gruppierung, die Anzahl der Gruppen) im Hamiltonian; es ist eine Eigenschaft des Hamiltonians und unterscheidet sich von der Anzahl der Operatoren, die in einen einzelnen Circuit gesampelt werden, unten mit $n$ bezeichnet.

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$ Hamiltonian-Terme, und da die Terme mit Zurücklegen gezogen werden, kann dasselbe $h_i$ mehr als einmal in einem einzelnen $V_k$ auftreten.

Die Größe:

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

ist die $L_1$-Norm der Koeffizienten, sodass jeder der $n$ Schritte für dieselbe Dauer $\lambda t / n$ entwickelt wird, unabhängig davon, welcher Term gezogen wurde. Die Gleichförmigkeit 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 Verteilung gesampelt:

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

sodass die Reihe $(k_1, \ldots, k_n)$ eine zufällige Sequenz von Term-Indizes ist, die aus dieser Verteilung gezogen werden. Da die $c_i$ positiv sind und sich zu $\lambda$ summieren, handelt es sich um eine normalisierte Wahrscheinlichkeitsverteilung, und der Erwartungswert des resultierenden Kanals über die zufälligen Ziehungen approximiert die Entwicklung unter $H$, mit einem Fehler, der abnimmt, wenn $n$ wächst. 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 zu unterscheiden.)

Dieses Tutorial zeigt, wie man ein Ensemble solcher randomisierter Circuits erzeugt. Nachdem wir diese Circuits erstellt haben, ähnlich wie wir einen Krylov-Unterraum für verschiedene Operatoren erzeugen, sampeln wir Bitstrings von mehreren solchen Operatoren mit unterschiedlichen Zeitparametern. Dies stellt eine höhere Überlappung zwischen den Grundzustandsvektoren und den gesampelten Bitstrings sicher.

Anforderungen

Bevor du mit diesem Tutorial beginnst, stelle sicher, dass du Folgendes installiert hast

• Eine Python-(>=3.10)-virtuelle Umgebung
• 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

Einrichtung

# Added by doQumentation — required packages for this notebook
!pip install -q numpy pyscf qiskit qiskit-addon-sqd qiskit-aer qiskit-fermions qiskit-ibm-runtime

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

Simulator-Beispiel

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Lesen und Vorbereiten des FCIDump

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

Über dieses 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 ein Full-Space-STO-3G-Dump handelt, werden 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 wird.

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 in Orbitalphase oder -reihenfolge von der mitgelieferten unterscheiden; die Gesamtenergien sind davon nicht betroffen.

Datei beschaffen. Du findest das FCIDump in diesem GitHub-Repository. Du kannst die untenstehende Zelle ausführen, um sie an den Ort zu holen, den der Rest des Tutorials erwartet.

Zunächst verwenden wir den von pyscf bereitgestellten cisolver, um die Referenzenergie zu erhalten. Dies ist die wahre Grundzustandsenergie des Moleküls, mit dem wir arbeiten. Dazu deklarieren wir zunächst norb und nelec, also die Anzahl der Orbitale beziehungsweise der Elektronen. Anschließend deklarieren wir h1e und h2e, also die Ein- beziehungsweise Zwei-Elektronen-Integrale. All dies wird 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 = "assets/sqdrift/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 = "assets/sqdrift/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

Laden des Hamiltonoperators

Nachdem die notwendigen Daten bereitstehen, lesen wir den Hamiltonoperator aus der FCI-Datei in einem Format, 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

Zunächst bilden wir den Hamiltonoperator mit qiskit-fermions auf ein fermionisches Schaltkreismodell ab, das transpilerspezifische Passes und Gates für fermionische Schaltkreise bereitstellt. Diese werden später vor Qiskits traditionellen Transpiler-Passes für diesen Workflow verwendet.

Term-Gruppierung

Um die Reproduzierbarkeit der Ergebnisse sicherzustellen, verwenden wir zunächst canonical_order, um die Terme ausschließlich anhand ihrer Struktur zu sortieren. Die Reihenfolge der Operatoren in der Liste canon ist damit festgelegt. Dies gewährleistet die Reproduzierbarkeit der erzeugten Operatoren, da der QDriftTrotterization-Pass, den wir später verwenden werden, zufällige Indizes sampelt, um die qDRIFT-Operatoren zu erzeugen.

In diesem Schritt nutzen wir die vielen Symmetrien, die im Hamiltonoperator der elektronischen Struktur vorhanden sind, indem wir verwandte Terme mit identischen Koeffizienten gruppieren. Auch wenn dies die Verteilung der Operatorkoeffizienten verändert, aus der das qDRIFT-Protokoll sampelt, beeinträchtigt dies nicht dessen Konvergenzgarantien. Entscheidend ist, dass die Gruppierung symmetrieverwandter Terme zu einer günstigen Auslöschung von Pauli-Termen und insgesamt zu einer kürzeren Schaltkreistiefe bei der Zeitentwicklung eines Zustands unter ihrer Wirkung führt.

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 von normalgeordneten Termen ausgeht.

Filtern diagonaler Terme

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

Bei den betreffenden Termen handelt es sich um jene, die in der Besetzungszahlbasis diagonal sind, also die Produkte von Besetzungszahloperatoren $a^\dagger_i a_i$. Drei Arten von Termen fallen unter diese Beschreibung:

• der konstante Energie-Offset, ein Produkt aus null Besetzungszahloperatoren, dessen Zeitentwicklung nur eine globale Phase beiträgt;

• die einzelnen Besetzungszahloperatoren $n_i$, deren Zeitentwicklung sich auf Einzel-Qubit-$Z$-Rotationen reduziert;

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

Für sich genommen verschiebt keiner dieser Terme 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 später im Schaltkreis durch die Anregungsterme erzeugt wird, sodass ihr Herausfiltern die tatsächlich erzeugte Zeitentwicklung verändert und die Sampling-Verteilung ändern kann. Dies ist eine bewusste Näherung im Schritt der Schaltkreiserzeugung, die vorgenommen wird, um das Sampling auf Anregungsterme zu fokussieren, und kein Schritt, der die gesampelte Verteilung unverändert lässt. Anders als die oben beschriebene Symmetriegruppierung, die die qDRIFT-Konvergenzgarantien intakt lässt, verändert dieser Filter den Operator, der zeitentwickelt wird. Die Schaltkreise approximieren daher nicht mehr die Zeitentwicklung unter dem vollständigen Hamiltonoperator, und die qDRIFT-Fehlerschranken gelten für den gefilterten statt für den ursprünglichen Operator. Dies ist hier akzeptabel, weil die Schaltkreise nur eine Sampling-Heuristik sind, die zum Vorschlagen von Konfigurationen dient: Kein Term geht aus der Energieschätzung selbst verloren, da der Filter nur auf den Hamiltonoperator angewendet wird, der zur Erzeugung der Schaltkreise dient, während die spätere klassische Diagonalisierung den vollständigen Hamiltonoperator einschließlich der diagonalen Terme 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 identifiziert sie anhand ihrer normalgeordneten Struktur — der Multimenge der Erzeugungsmoden, die der Multimenge der Vernichtungsmoden entspricht —, weshalb sie nur für einen bereits normalgeordneten Operator gültig ist. Diese Annahme wird zur Laufzeit nicht überprü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 fest, um das Ensemble von Schaltkreisen zu erzeugen:

• Die Anzahl der zu erzeugenden Schaltkreise: num_circuits
• Die Länge jedes Schaltkreises in Anregungsgruppen: num_exc
• Der Faktor für die unterschiedlichen Zeitentwicklungszeiten: times

Erzeugung fermionischer Schaltkreise

Wir erzeugen nun fermionische Schaltkreise für jeden der Zeitschritte. Jeder Schaltkreis besteht aus einem einzelnen Evolution-Gate mit der zuvor deklarierten Zeitentwicklungszeit. Der Zeitentwicklungsoperator ist der Hamiltonoperator. Später führen wir Transpiler-Passes auf diesen Schaltkreisen aus, um qDRIFT-Schaltkreise zu erzeugen.

Ansatz-Vorbereitung

Wir bereiten den Hartree-Fock-Zustand mithilfe der Klasse InitializeModes vor. Für Stickstoff besteht der Prozess einfach darin, X-Gates zunächst auf die ersten num_elec_a-Qubits und dann auf die num_elec_b-Qubits anzuwenden, wobei 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: Problem für die Ausführung auf Quantenhardware optimieren

Nachdem wir unsere Schaltkreise haben, verwenden wir zunächst die in qiskit-fermions verfügbaren Passes, um fermionische Optimierungen durchzuführen, gefolgt vom Transpilieren unseres Schaltkreises für das gewählte Backend. Da es sich hier um ein Simulator-Experiment handelt, tun wir dies zunächst für den AerSimulator. Gewichtsberechnung für jede Gruppe

In diesem Schritt führen wir das stochastische qDRIFT-Sampling von Termen mit Wahrscheinlichkeiten proportional zu ihren Koeffizienten im Hamiltonoperator durch. Der qDRIFT-Transpiler-Pass übernimmt dies für uns. Wir können nun flachere Schaltkreise erzeugen, die trotz begrenzter Qubit-Konnektivität effizienter auf der Hardware ausgeführt werden können, selbst wenn der Hamiltonoperator langreichweitige Kopplungen und Terme höherer als quadratischer Ordnung enthält. Nach der Term-Gruppierung sampelt es die Operatoren basierend auf ihren Gewichten. Für jeden Operator $h_i$ ist das Gewicht $W_{h_i}$ wie folgt definiert:

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

Fermionische und hardware-native Optimierungen

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

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

• Der RelabelModes-Pass ist ein weiterer Optimierungspass, der verwendet werden kann, um die fermionischen Modi zu permutieren, um die Konnektivität über die Qubits zu optimieren und die Gate-Tiefe zu reduzieren; mehr dazu in der API-Referenz

Die verbleibenden 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, um fermionenbasierte Schaltkreisinstruktionen auf qubitbasierte abzubilden.

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 wir mit den Optimierungen auf fermionischer Ebene fertig sind, können wir die Schaltkreise für die Ausführung auf dem Simulator transpilieren.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Schritt 3: Ausführung mit Qiskit-Primitives

Nachdem wir unsere Schaltkreise haben, können wir sie mit Qiskit-Primitives auf dem AerSimulator ausführen. Wir kombinieren alle Counts der verschiedenen Schaltkreise. Wir wandeln sie in boolesche Vektoren um, bevor wir sie schließlich mit SQD nachbearbeiten.

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: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format

Verwendung von Bitstrings für SQD

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

Um schließlich die Grundzustandsschätzung 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 Teilsample zieht eine andere Menge von Konfigurationen, und der Recovery-Schritt formt den Pool zwischen den Iterationen um, sodass die gemeldete Dimension von einem Teilsample zum nächsten variiert. Rauschfreies Sampling legt die Dimension des ausgewählten Unterraums für sich genommen nicht fest. Der Hardware-Lauf liefert jedoch tendenziell systematisch größere Unterräume, da verrauschte Shots die Teilchenzahlsymmetrie brechen und die Konfigurations-Recovery sie in zusätzliche Basisvektoren verwandelt. Aus diesem Grund führen wir im Hardware-Abschnitt auch einen weiteren Schritt zum Bereinigen 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

Hardware-Beispiel

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

Die Kosten des klassischen Schritts werden nicht direkt durch die Anzahl der Qubits bestimmt. SQD diagonalisiert den Hamiltonoperator projiziert auf den Unterraum, der von den gesampelten Konfigurationen aufgespannt wird. Was die klassischen Kosten bestimmt, ist daher die Dimension dieses ausgewählten Unterraums — hier gesteuert durch samples_per_batch, num_batches und die Anzahl der unterschiedlichen Konfigurationen, die die Schaltkreise tatsächlich erzeugen — zusammen mit der dünnbesetzten linearen Algebra, die zur Anwendung des projizierten Hamiltonoperators benötigt wird. Der vollständige CI-Raum wächst kombinatorisch mit Orbitalen und Elektronen, aber der ausgewählte Unterraum ist ein kleiner, anpassbarer Ausschnitt davon, und wir steuern seine Größe direkt. Folglich können die Anzahl der Qubits und die klassische Schwierigkeit bis zu einem gewissen Grad unabhängig voneinander variiert werden: 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 daher von der Dimension des Unterraums ab, die du für die gewünschte Genauigkeit benötigst, sowie vom verfügbaren Speicher und den verfügbaren Kernen für den Eigenlöser. Größere Orbitalräume erfordern typischerweise 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. Anstatt von einer festen Obergrenze auszugehen, 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 des Sampling-Fehlers durch das Rauschen in der Hardware wird der für die Diagonalisierung im Hardware-Lauf erzeugte Unterraum größer sein als der, den wir bei Verwendung des Simulators erhalten. Auch wenn dies die Dimension des zu diagonalisierenden Unterraums erhöht, liefert der Workflow dank der Robustheit von SQD gegenüber Rauschen dennoch eine genaue Antwort.

Bereinigung fehlerhafter Strings

Hier können wir uns entscheiden, einen zusätzlichen Schritt durchzuführen. Sobald wir alle Bitstrings aus den Schaltkreisausführungen haben, können wir entweder die ungültigen Bitstrings herausfiltern, bevor wir SQD ausführen, oder ohne Bereinigung fortfahren. Das Überspringen der Bereinigung ist bei Hardware-Läufen im Allgemeinen vorzuziehen, da dadurch die symmetriegebrochenen Shots für die Konfigurations-Recovery verfügbar bleiben, die sie in gültige Konfigurationen umwandeln und dadurch den Unterraum erweitern kann, anstatt diese Shots einfach zu verwerfen.

Da Stickstoff nur sieben $\alpha$- und sieben $\beta$-Elektronen haben kann, können alle Bitstrings, die mehr oder weniger als sieben Einsen in der ersten und der zweiten Hälfte der Ausgabe enthalten, verworfen werden. 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 in das Diagonalisierungsschema eingespeist. Verwende das Flag PRUNE unten, um zwischen den beiden Verhaltensweisen zu wechseln.

Denke daran, dass Bereinigung nur eine von mehreren Entscheidungen ist, die den finalen Unterraum formen, neben der Anzahl der Schaltkreise, der Menge der Zeitentwicklungszeiten und der Filterung diagonaler Terme. Ein bereinigter Lauf im Vergleich zu einem unbereinigten ist nur aussagekräftig, wenn alles andere konstant gehalten wird; der C++-Begleitcode geht darauf genauer ein, da er postselektiert statt zu recovern und sich auch in diesen anderen Parametern unterscheidet.

name = "assets/sqdrift/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, könnten dich folgende Materialien interessieren:

• Sample-basierte Krylov-Quantendiagonalisierung eines fermionischen Gittermodells - ein verwandtes Tutorial, das Zeitentwicklungsschaltkreise statt eines variationellen Ansatzes verwendet.
• Sample-basierte Quantendiagonalisierung eines chemischen Hamiltonoperators - ein Tutorial darüber, wie man einen Local-Unitary-Cluster-Jastrow-(LUCJ)-Schaltkreis für die Simulation von Quantenchemie konstruiert.
• Das SqDRIFT-Paper - die Literatur, auf der dieses Tutorial basiert. (Beachte, dass einige der in diesem Paper diskutierten Optimierungen derzeit in Arbeit sind und sich dieses Tutorial in Zukunft basierend auf der Weiterentwicklung der verwendeten Bibliotheken ändern kann.)