SqDRIFT-Algorithmus zur Grundzustandsabschätzung

Python

Python-Notebook herunterladen

Lade das Python-Notebook herunter, um dieses Tutorial selbst auszuführen. Das Notebook lädt die FCIDump-Eingabe herunter, falls sie noch nicht vorhanden ist.

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

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

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-Hamiltonoperator konstruiert werden. Dies wird erreicht, indem kleinere Zeitentwicklungsoperatoren anhand seiner Koeffizienten aus dem Hamiltonoperator gezogen werden, was als qDRIFT-Trotterisierungsmethode bekannt ist.

Beide Implementierungen verwenden Qiskit Fermions, um den molekularen Hamiltonian darzustellen und Circuits für den qDRIFT-Algorithmus zu konstruieren. Die Python-Implementierung nutzt fermionische Transpiler-Passes, während die C++-Implementierung Gruppen explizit samplet und auf Qubit-Operatoren abbildet.

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 Hamiltonoperator; es ist eine Eigenschaft des Hamiltonoperators 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$ Hamiltonoperator-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

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

Implementierung

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Die folgenden Python-Schritte verwenden den Aer-Simulator. Das Hardware-Beispiel folgt danach.

Lesen und Vorbereiten des FCIDump

Für dieses Tutorial laden wir den Hamiltonoperator 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 im nächsten Code-Block 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 den untenstehenden Code 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 = "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}")

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

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, das Transpiler-Passes und Gates speziell für fermionische Circuits bereitstellt, auf ein fermionisches Circuit-Modell ab. 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 Circuit-Tiefe 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 Circuits 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 Circuit 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 Circuit-Erzeugung, 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 Circuits 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 Circuits 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 Circuits 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 Circuits zu erzeugen:

• Die Anzahl der zu erzeugenden Circuits: num_circuits
• Die Länge jedes Circuits in Anregungsgruppen: num_exc
• Der Faktor für die unterschiedlichen Evolutionszeiten: 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 Circuits haben, verwenden wir zunächst die in qiskit-fermions verfügbaren Passes, um fermionische Optimierungen durchzuführen, gefolgt vom Transpilieren unseres Circuits 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 Circuits 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$$

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 auf sein Vorzeichen reduziertem Koeffizienten evolviert.

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 Circuit 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 Circuits 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 Circuits 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 Wiederherstellungsschritt 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 Konfigurationswiederherstellung 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 Circuits 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 Circuit-Ausfü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 Konfigurationswiederherstellung 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 Circuits, der Menge der Evolutionszeiten 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++-Tab geht darauf genauer ein, da er Postselektion statt Konfigurationswiederherstellung verwendet 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")

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

C++

C++-Code herunterladen

Der C++-Quellcode, das CMake-Projekt und die Build-Anweisungen befinden sich im Verzeichnis docs/tutorials/assets/sqdrift/cpp des Qiskit-Dokumentations-Repositorys. Die FCIDump-Eingabe liegt im benachbarten Verzeichnis fcidump_files. Klone das Repository oder lade beide Verzeichnisse herunter, um dieses Tutorial selbst zu bauen und auszuführen.

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

Lernziele (C++)

• 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

Voraussetzungen (C++)

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

• Lies die Lektion zu Sample-based Krylov Quantum Diagonalization (SKQD)

Hintergrund (C++)

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 seiner Koeffizienten aus dem Hamiltonoperator gezogen werden, was als qDRIFT-Trotterisierungsmethode bekannt ist.

Beide Implementierungen verwenden Qiskit Fermions, um den molekularen Hamiltonian darzustellen und Circuits für den qDRIFT-Algorithmus zu konstruieren. Die Python-Implementierung nutzt fermionische Transpiler-Passes, während die C++-Implementierung Gruppen explizit samplet und auf Qubit-Operatoren abbildet.

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 Hamiltonoperator; es ist eine Eigenschaft des Hamiltonoperators 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$ Hamiltonoperator-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 (C++)

Diese Implementierung läuft auf IBM Quantum®-Hardware und verwendet für die klassische Diagonalisierung eine externe SBD-Programmdatei. Das Python-Simulatorbeispiel ist darin nicht enthalten.

Folge den C++-Build-Anweisungen, um den Compiler, CMake, Rust, Boost und die Qiskit-Abhängigkeiten zu installieren. Installiere Rust mit rustup: Die Qiskit-Abhängigkeiten legen ihre eigene Rust-Version fest, die rustup automatisch herunterlädt, aber ein über Homebrew oder einen Linux-Paketmanager installiertes Rust ist möglicherweise zu alt. Der abschließende SBD-Schritt erfordert zusätzlich MPI, OpenMP und BLAS/LAPACK. Windows wird von dieser Implementierung derzeit nicht unterstützt.

Verwende den vollständigen C++-Quellcode und das CMake-Projekt aus dem Repository. Behalte beim Herunterladen des Beispiels das danebenliegende Verzeichnis fcidump_files bei. Baue aus einem Checkout des Dokumentations-Repositorys mit den folgenden Befehlen. Der erste Build lädt alle Abhängigkeiten herunter, kompiliert sie und dauert mehrere Minuten.

cd docs/tutorials/assets/sqdrift/cpp
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --config Release

Der Build schreibt die ausführbare Datei SqDRIFT nach docs/tutorials/assets/sqdrift/cpp (das Projektstammverzeichnis), nicht nach build/. Führe ./SqDRIFT im Projektstammverzeichnis aus, damit es ../fcidump_files/N2_sto_3g findet. Alle anderen Befehle in diesem Tab werden ebenfalls im Projektstammverzeichnis ausgeführt, sofern nichts anderes angegeben ist. Wenn ./SqDRIFT mit einem Fehler wegen einer fehlenden Shared Library fehlschlägt, folge den plattformspezifischen Anweisungen zum Bibliothekspfad in der Build-Anleitung; unter macOS läuft die ausführbare Datei normalerweise auch ohne sie.

Einrichtung (C++)

Die folgenden Snippets sind Auszüge aus dem vollständigen SqDRIFT.cpp-Programm. Baue und führe diese Datei als ein einzelnes Programm aus.

Stelle vor dem Ausführen von SqDRIFT sicher, dass dein IBM Quantum API-Token für die IBM Cloud Platform konfiguriert ist.

Wenn du Python mit installiertem qiskit-ibm-runtime hast, erstellst du am einfachsten eine korrekt formatierte Konfigurationsdatei so:

from qiskit_ibm_runtime import QiskitRuntimeService
QiskitRuntimeService.save_account(channel="ibm_quantum_platform", token="YOUR_IBM_QUANTUM_TOKEN_HERE")

Das schreibt ~/.qiskit/qiskit-ibm.json (Linux/macOS) bzw. %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) in der richtigen Form. Wenn du die Datei lieber manuell erstellst:

{
  "default-ibm-quantum-platform": {
    "channel": "ibm_quantum_platform",
    "token": "YOUR_IBM_QUANTUM_TOKEN_HERE"
  }
}

Hole dir deinen API-Schlüssel von IBM Quantum Platform.

Hinweis: Eine fehlende oder fehlerhafte Konfigurationsdatei führt beim Start zu einer wenig hilfreichen Rust-unwrap()-Panik statt zu einer sauberen Fehlermeldung. Prüfe den Schlüsselnamen (default-ibm-quantum-platform) und den Channel-Wert (ibm_quantum_platform) sorgfältig, wenn ein Problem auftritt.

Implementierung (C++)

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden (C++)

Den molekularen Hamiltonoperator laden

char filename[] = "../fcidump_files/N2_sto_3g";
QfFCIDump* fcidump = qf_fcidump_from_file(filename);
if (fcidump == nullptr) {
    std::cerr << "Failed to load FCIDump!" << std::endl;
    return 1;
}

uint32_t norb = qf_fcidump_norb(fcidump);
uint32_t nelec = qf_fcidump_nelec(fcidump);
uint32_t n_alpha = nelec / 2;
uint32_t n_beta  = nelec - n_alpha;
uint32_t num_modes = 2 * norb;
uint32_t num_qubits = num_modes;
std::cout << "✓ Loaded N2 molecule (" << norb << " orbitals, "
          << num_qubits << " qubits, " << nelec << " electrons)" << std::endl;

Output:

✓ Loaded N2 molecule (10 orbitals, 20 qubits, 14 electrons)

Lädt den Hamiltonian des $N_2$-Moleküls aus einer FCIDump-Datei. Die Anzahl der Qubits entspricht der doppelten Anzahl der Orbitale (Spin-up und Spin-down).

Den Hamiltonian erstellen und normalordnen

// Create and normal order Hamiltonian
QfFermionOperator* hamiltonian = qf_ferm_op_from_fcidump(fcidump);
QfFermionOperator* normal_ordered = qf_ferm_op_normal_ordered(hamiltonian, nullptr);

Wandelt die FCIDump-Daten in einen fermionischen Operator um und wendet Normalordnung an, um die Operatorstruktur zu vereinfachen.

Terme nach elektronischer Struktur gruppieren

// Group terms by electronic structure
QfExitCode exit_code = qf_ferm_op_group_terms_by_electronic_structure(
    normal_ordered, num_modes, false
);
if (exit_code != QfExitCode_Success) {
    std::cerr << "Failed to group terms!" << std::endl;
    return 1;
}

uint32_t num_groups = qf_ferm_op_num_groups(normal_ordered);
std::cout << "Grouped into " << num_groups << " groups" << std::endl;

// Split into group operators
QfFermionOperator** group_ops = new QfFermionOperator*[num_groups];
qf_ferm_op_split_out_groups(normal_ordered, nullptr, 0, group_ops);

Output:

Grouped into 1902 groups

Gruppiert zusammengehörige Hamiltonian-Terme nach elektronischer Struktur, bevor ihre Evolutions-Circuits konstruiert werden.

Sampling-Gewichte berechnen, dann normalisieren und auf Qubit-Operatoren abbilden

Die Reihenfolge ist hier wichtig. Die Normalisierung skaliert jeden Koeffizienten an Ort und Stelle auf Betrag eins, daher müssen die Sampling-Gewichte zuerst berechnet werden. Läuft die Normalisierung zuerst, liest sich jedes $|c_j|$ als exakt 1.0 zurück, und das Gewicht jeder Gruppe schrumpft auf die bloße Anzahl ihrer Terme zusammen — die Sampling-Verteilung würde dann von der Gruppengröße statt von den physikalischen Koeffizienten bestimmt.

Das Gewicht jeder Gruppe ist der Mittelwert der Absolutwerte der Koeffizienten ihrer Terme. Der QDriftTrotterization-Pass des Python-Tabs verwendet dasselbe Gewicht für einen gruppierten Hamiltonian. Würde man stattdessen die Summe verwenden, skaliert $\lambda$ mit der Gruppengröße, was wiederum die Evolutionszeit $t_k = \lambda\tau/n$ verzerrt.

// Calculate sampling weights BEFORE normalization.
std::vector<double> weights(num_groups);
double total_weight = 0.0;

for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    double abs_sum = 0.0;
    for (uint64_t j = 0; j < num_terms; j++) {
        abs_sum += std::sqrt(coeffs[j].re * coeffs[j].re +
                             coeffs[j].im * coeffs[j].im);
    }
    weights[i] = (num_terms > 0) ? (abs_sum / static_cast<double>(num_terms)) : 0.0;
    total_weight += weights[i];
}

std::cout << "Total weight (λ): " << total_weight << std::endl;

// Normalize the evolution operators term-by-term, now that the sampling
// weights above have already captured the original magnitudes.
for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    for (uint64_t j = 0; j < num_terms; j++) {
        const double magnitude = std::sqrt(coeffs[j].re * coeffs[j].re +
                                           coeffs[j].im * coeffs[j].im);
        if (magnitude > 0.0) {
            coeffs[j].re /= magnitude;
            coeffs[j].im /= magnitude;
        }
    }
}

QkObs** qubit_ops = new QkObs*[num_groups];
for (uint32_t i = 0; i < num_groups; i++) {
    QfExitCode jw_exit = qf_ferm_op_jordan_wigner(group_ops[i], num_qubits, &qubit_ops[i]);
    if (jw_exit != QfExitCode_Success) {
        std::cerr << "Failed to map group " << i << std::endl;
        return 1;
    }
}
std::cout << "Mapped all " << num_groups << " normalized groups to qubit operators" << std::endl;

Output:

Total weight (λ): 337.261
Mapped all 1902 normalized groups to qubit operators

Das Gewicht $\lambda_i = \frac{1}{|G_i|}\sum_{j \in G_i} |c_{ij}|$ ist der mittlere Absolutwert der Koeffizienten über die Terme der Gruppe $G_i$, gelesen aus den ursprünglichen Koeffizienten. Die Normalisierung skaliert dann vor der Jordan-Wigner-Transformation jeden fermionischen Term auf Betrag eins und bewahrt im evolvierten Operator nur seine Phase bzw. sein Vorzeichen. Das entspricht der Evolution normalisierter Terme in qDRIFT: Die Beträge der Koeffizienten bestimmen über $\lambda$ die Sampling-Wahrscheinlichkeiten, während die Circuit-Evolution den normalisierten Term verwendet, sodass große Koeffizienten nicht doppelt gezählt werden.

Operatormengen samplen

// SqDRIFT Sampling: Create operator sets
    const int num_circuits = 100;    // Number of circuits to create
    const int ops_per_circuit = 10;  // Operators per circuit
    const double time_step = 1;    // Time step for evolution
    
    std::cout << "\n Generating " << num_circuits << " circuits with " 
              << ops_per_circuit << " operators each..." << std::endl;
    
    std::mt19937 gen(42);
    std::discrete_distribution<> dist(weights.begin(), weights.end());
    
    // Store sampled operator indices
    std::vector<std::vector<int>> operator_sets(num_circuits);
    for (int i = 0; i < num_circuits; i++) {
        operator_sets[i].resize(ops_per_circuit);
        for (int j = 0; j < ops_per_circuit; j++) {
            operator_sets[i][j] = dist(gen);
        }
    }

Output:

 Generating 100 circuits with 10 operators each...

Zieht die Operatormengen für 100 zufällige Circuits, die jeweils 10 aus der gewichteten Verteilung gesampelte Gruppen enthalten. Der nächste Schritt macht aus diesen Mengen das Ensemble der qDRIFT-Circuits.

Schritt 2: Die qDRIFT-Circuits erstellen (C++)

In der C++-Version baut dieser Schritt die Circuits aus den gesampelten Operatormengen auf. Die Transpilation für das Ziel-Backend erfolgt in Schritt 3, unmittelbar bevor jeder Circuit eingereicht wird.

    // Create Suzuki-Trotter circuits by composing all operators
    std::vector<Qiskit::circuit::QuantumCircuit> circuits;
    circuits.reserve(num_circuits);

    const uint32_t trotter_order = 1;
    const uint32_t trotter_reps = 1;
    const bool preserve_order = false;
    const bool insert_barriers = false;

    // Gate name mapping is static — compute once before all loops
    auto name_map = Qiskit::circuit::get_standard_gate_name_mapping();

    for (int circuit_idx = 0; circuit_idx < num_circuits; circuit_idx++) {
        // Create initial circuit with state preparation using C++ API
        Qiskit::circuit::QuantumCircuit qc(num_qubits, num_qubits);

        // State preparation: apply X gates to spin-up qubits 0..n_alpha-1
        // and spin-down qubits norb..norb+n_beta-1 (derived from FCIDump NELEC)
        for (uint32_t q = 0; q < n_alpha && q < num_qubits; q++) {
            qc.x(q);
        }
        for (uint32_t q = norb; q < norb + n_beta && q < num_qubits; q++) {
            qc.x(q);
        }

        // For each operator in this circuit's set
        for (int op_idx = 0; op_idx < ops_per_circuit; op_idx++) {
            int group_idx = operator_sets[circuit_idx][op_idx];
            double evolution_time = (total_weight * time_step) / ops_per_circuit;
            
            // Create Suzuki-Trotter evolution circuit
            QkCircuit* evolution_circuit = qk_circuit_library_suzuki_trotter(
                qubit_ops[group_idx], trotter_order, trotter_reps,
                evolution_time, preserve_order, insert_barriers
            );
            
            if (evolution_circuit != nullptr) {
                // Get number of instructions
                uint32_t num_ops = qk_circuit_num_instructions(evolution_circuit);
                
                // Manually append each instruction
                for (uint32_t i = 0; i < num_ops; i++) {
                    QkCircuitInstruction inst;
                    qk_circuit_get_instruction(evolution_circuit, i, &inst);

                    // Prepare qubit vector
                    std::vector<uint32_t> qubits(inst.num_qubits);
                    for (uint32_t j = 0; j < inst.num_qubits; j++) {
                        qubits[j] = inst.qubits[j];
                    }

                    // Get operation kind
                    QkOperationKind kind = qk_circuit_instruction_kind(evolution_circuit, i);

                    // Get mutable pointer to circuit
                    QkCircuit* mutable_circuit = qc.get_rust_circuit().get();

                    // Append based on operation type
                    if (kind == QkOperationKind_Gate) {
                        // Convert gate name string to QkGate enum
                        std::string gate_name(inst.name);
                        auto gate_it = name_map.find(gate_name);
                        if (gate_it != name_map.end()) {
                            QkGate gate_enum = gate_it->second.gate_map();
                            qk_circuit_parameterized_gate(
                                mutable_circuit,
                                gate_enum,           // Use QkGate enum, not string
                                qubits.data(),
                                inst.params
                            );
                        } else {
                            std::cerr << "Unknown gate in evolution circuit: " << gate_name << std::endl;
                            qk_circuit_instruction_clear(&inst);
                            qk_circuit_free(evolution_circuit);
                            return 1;
                        }
                    } else if (kind == QkOperationKind_Barrier) {
                        qk_circuit_barrier(
                            mutable_circuit,
                            qubits.data(), 
                            inst.num_qubits
                        );
                    } else if (kind == QkOperationKind_Reset) {
                        qk_circuit_reset(
                            mutable_circuit,
                            qubits[0]
                        );
                    }

                    qk_circuit_instruction_clear(&inst);
                }

                qk_circuit_free(evolution_circuit);
            }
        }

        circuits.push_back(qc);

        // Display progress
        if ((circuit_idx + 1) % 10 == 0) {
            std::cout << "  Created " << (circuit_idx + 1) << "/" 
                    << num_circuits << " circuits" << std::endl;
        }
    }

    std::cout << "✓ Created all " << num_circuits << " Suzuki-Trotter circuits" << std::endl;

Output:

  Created 10/100 circuits
  Created 20/100 circuits
  Created 30/100 circuits
  Created 40/100 circuits
  Created 50/100 circuits
  Created 60/100 circuits
  Created 70/100 circuits
  Created 80/100 circuits
  Created 90/100 circuits
  Created 100/100 circuits
✓ Created all 100 Suzuki-Trotter circuits

Das vollständige Programm gibt anschließend Circuit-Statistiken und die Operatorindizes von Circuit 0 aus, die hier weggelassen sind.

Für jeden gesampelten Operator wird ein Suzuki-Trotter-Evolutions-Circuit mit der Zeit $t_k = \frac{\lambda \tau}{n}$ erstellt, wobei $n$ die Anzahl der pro Circuit gesampelten Operatoren und $\tau$ der Zeitschritt ist. Da jede gesampelte Gruppe vor der Abbildung term-weise normalisiert wurde, wendet diese Evolution bei der Rotation nur die Phase bzw. das Vorzeichen jedes Terms an; die Beträge gehen nur über die Sampling-Verteilung ein. Der Hartree-Fock-Anfangszustand wird durch Anwenden von X-Gates auf die Qubits 0-6 und 10-16 präpariert (sieben Elektronen in jedem Spinsektor). Die Anweisungen jedes Evolutions-Circuits werden mit der C++-API manuell an den Haupt-Circuit angehängt.

Schritt 3: Ausführung mit Qiskit-Primitives (C++)

std::cout << "\n Connecting to IBM Quantum Cloud..." << std::endl;

// Initialize IBM Quantum Compute Service using qiskit-cpp
Qiskit::service::QiskitRuntimeService service;
std::cout << "✓ Connected to IBM Quantum" << std::endl;

// Get backend
const std::string backend_name = "ibm_fez";
auto backend = service.backend(backend_name);

if (backend.name().empty()) {
    std::cerr << "Backend " << backend_name
              << " is not available to this account." << std::endl;
    return 1;
}

std::cout << "✓ Selected backend: " << backend.name() << std::endl;

// Create sampler primitive
const int32_t shots = 100;  // Number of shots per circuit
Qiskit::primitives::BackendSamplerV2 sampler(backend, shots);

std::cout << "\n Submitting " << num_circuits << " circuits ("
            << shots << " shots each)..." << std::endl;

// Submit one job per circuit.
std::vector<std::shared_ptr<Qiskit::primitives::BasePrimitiveJob>> jobs;

for (int i = 0; i < num_circuits; i++) {
    for (uint32_t q = 0; q < num_qubits; q++) {
        circuits[i].measure(q, q);
    }
    Qiskit::circuit::QuantumCircuit transpiled_qc =
        Qiskit::compiler::transpile(circuits[i], backend, 2, 1.0, 42);
    std::vector<Qiskit::primitives::SamplerPub> pubs;
    pubs.push_back(Qiskit::primitives::SamplerPub(transpiled_qc, shots));
    jobs.push_back(sampler.run(pubs));
}

std::cout << "\n⏳ Waiting for " << num_circuits << " jobs to complete..." << std::endl;

// Collect results — poll each job and extract bitstrings
std::vector<boost::dynamic_bitset<>> all_bitstrings;

for (int i = 0; i < num_circuits; i++) {
    while (!jobs[i]->in_final_state()) {
        sleep(5);
    }
    if (jobs[i]->done()) {
        auto result = jobs[i]->result();
        if (result.size() > 0) {
            auto& pub_result = result[0];
            auto& bit_array  = pub_result.data();
            for (size_t s = 0; s < bit_array.num_shots(); s++) {
                auto sample = bit_array[s].to_string();
                boost::dynamic_bitset<> bs(num_qubits);
                for (uint32_t bit = 0; bit < num_qubits; bit++) {
                    bs[bit] = (sample[num_qubits - 1 - bit] == '1');
                }
                all_bitstrings.push_back(bs);
            }
        }
    } else {
        std::cout << "  Circuit " << (i+1) << " did not complete successfully" << std::endl;
    }
}

std::cout << "\n✓ Collected " << all_bitstrings.size()
            << " total bitstrings from all circuits" << std::endl;

Verbindet sich mit IBM Quantum, wählt das Backend nach Namen aus, transpiliert jeden Circuit für dieses Backend und reicht einen Job pro Circuit ein: 100 einzelne Sampler-Jobs mit je 100 Shots auf einem 20-Qubit-Register. Anschließend fragt es jeden Job ab, bis er abgeschlossen ist, und sammelt die Bitstrings aus result[0]. Jeder Job verbraucht QPU-Zeit und wartet separat in der Warteschlange, daher dauert dieser Schritt deutlich länger als ein einzelner gebündelter Job.

Die auf dieser Seite genannten Hardware-Ergebnisse stammen aus einem Lauf auf ibm_pittsburgh (100 Circuits × 100 Shots), wobei backend_name gegenüber dem Standardwert ibm_fez geändert wurde. Die gespeicherten Samples wurden mit der in Schritt 4 gezeigten korrigierten Bitreihenfolge der Determinanten nachverarbeitet. Das Hardware-Rauschen variiert zwischen den Läufen, daher unterscheiden sich deine Counts und die finale Energie.

Output:

 Connecting to IBM Quantum Cloud...
✓ Connected to IBM Quantum
✓ Selected backend: ibm_pittsburgh

 Submitting 100 circuits (100 shots each)...

⏳ Waiting for 100 jobs to complete...

✓ Collected 10000 total bitstrings from all circuits

Das Backend ändern: Ändere backend_name in SqDRIFT.cpp, um ein anderes Gerät anzusteuern; service.backends() listet auf, was dein Konto erreichen kann.

Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format (C++)

Bitstrings postselektieren und CI-Strings schreiben

// Postselect bitstrings by Hamming weight (derived from FCIDump NELEC)
std::vector<double> bitstring_weights(all_bitstrings.size(), 1.0);
auto [filtered_bitstrings, filtered_weights] = Qiskit::addon::sqd::postselect_bitstrings(
    all_bitstrings,
    bitstring_weights,
    Qiskit::addon::sqd::MatchesRightLeftHamming<uint32_t>(n_alpha, n_beta)
);
std::ignore = filtered_weights;

std::cout << "✓ Postselected " << filtered_bitstrings.size()
          << " bitstrings with Hamming weight ("
          << n_alpha << "," << n_beta << ")" << std::endl;

// Convert bitstrings to CI strings using SQD addon
auto ci_strings = Qiskit::addon::sqd::bitstrings_to_ci_strings_symmetrize_spin(
    filtered_bitstrings,
    std::nullopt  // No dimension limit
);

std::cout << "✓ Generated " << ci_strings.size() << " CI strings" << std::endl;

// Write CI strings to file for SBD
std::ofstream alpha_file("alphadets_from_sqd.txt");
for (const auto& ci_string : ci_strings) {
    std::string bitstr;
    // SBD expects the highest orbital index first (bit zero on the right).
    boost::to_string(ci_string, bitstr);
    alpha_file << bitstr << "\n";
}
alpha_file.close();

std::cout << "✓ Wrote " << ci_strings.size()
          << " CI strings to alphadets_from_sqd.txt" << std::endl;
std::cout << " Ready for SBD diagonalization!" << std::endl;

Output:

✓ Postselected 1137 bitstrings with Hamming weight (7,7)
✓ Generated 49 CI strings
✓ Wrote 49 CI strings to alphadets_from_sqd.txt
 Ready for SBD diagonalization!

Postselektiert Bitstrings mit dem korrekten Hamming-Gewicht (n_alpha Spin-up-, n_beta Spin-down-Elektronen, abgeleitet aus dem FCIDump-Feld NELEC), wandelt sie mittels Spinsymmetrisierung in Configuration-Interaction-(CI-)Strings um und schreibt sie für die anschließende Selected Basis Diagonalization nach alphadets_from_sqd.txt.

Ressourcen freigeben

// Cleanup
// circuits vector will be automatically cleaned up (RAII)
// Free qubit operators
for (uint32_t i = 0; i < num_groups; i++) {
    qk_obs_free(qubit_ops[i]);
}
delete[] qubit_ops;
// Free fermionic operators
for (uint32_t i = 0; i < num_groups; i++) {
    qf_ferm_op_free(group_ops[i]);
}
delete[] group_ops;
qf_ferm_op_free(normal_ordered);
qf_ferm_op_free(hamiltonian);
qf_fcidump_free(fcidump);
return 0;

Gibt den gesamten allokierten Speicher und alle Ressourcen ordnungsgemäß frei. Der C++-std::vector<QuantumCircuit> wird per RAII automatisch aufgeräumt. Qubit-Operatoren, fermionische Operatoren und die FCIDump-Daten werden manuell freigegeben.

Nachdem wir nun die Basis erstellt haben, auf die wir unseren Hamiltonian projizieren, können wir mit der Diagonalisierung fortfahren, um die Grundzustandsabschätzung zu erhalten. Führe die folgenden Befehle im Projektstammverzeichnis docs/tutorials/assets/sqdrift/cpp aus.

Die SBD-diag-Binärdatei bauen

CMake lädt SBD nur als Quellcode nach deps/sbd/ (der Haupt-Build überspringt das Kompilieren bewusst). Die SBD-App verwendet ein einfaches Makefile, und SBD selbst benötigt MPI, OpenMP und BLAS/LAPACK — ein einfaches make genügt also erst, wenn diese installiert sind und die Configuration-Datei der App zu deiner Toolchain passt.

Das Makefile liest seine Compiler- und Linker-Flags aus der Configuration-Datei im selben Verzeichnis. Die im Upstream mitgelieferte Version ist auf ein bestimmtes macOS/Homebrew-Layout zugeschnitten und lässt sich auf einem Standardrechner nicht unverändert bauen — vor allem enthält -I/opt/homebrew/opt/llvm/include kein omp.h (es liegt in einem versionsspezifischen Clang-Ressourcenverzeichnis), und Apples clang++ lehnt -fopenmp ganz ab:

clang++: error: unsupported option '-fopenmp'

Das folgende Rezept ist auf macOS (Apple Silicon, Darwin 25.6, Open MPI 5.x + Homebrew LLVM) nachweislich funktionsfähig.

macOS (Apple Silicon) — getestet

# 1. Install MPI and an OpenMP-capable compiler.
#    Apple clang cannot do OpenMP, so Homebrew's LLVM provides both
#    clang++ and libomp. BLAS/LAPACK come from Apple's Accelerate
#    framework, which is already part of macOS — no openblas needed.
brew install open-mpi llvm

Ersetze dann deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization/Configuration durch:

# Path to the SBD library
SBD_PATH=../..

# MPI C++ compiler. OMPI_CXX points Open MPI's wrapper at Homebrew's clang++,
# which (unlike Apple clang) supports -fopenmp.
CCCOM=OMPI_CXX=/opt/homebrew/opt/llvm/bin/clang++ mpicxx

# Build flags. omp.h lives in a version-specific clang resource dir, so let
# Homebrew's clang++ find it itself rather than hardcoding an -I path.
CCFLAGS= -std=c++17 -fopenmp -O3

# Link flags: OpenMP runtime from Homebrew LLVM, BLAS/LAPACK from the
# Accelerate framework that ships with macOS (no openblas install needed).
SYSLIB= -L/opt/homebrew/opt/llvm/lib -lomp -framework Accelerate

# 2. Build
cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization
make
cd ../../../..

Ein erfolgreicher Build gibt die beiden Compile-/Link-Zeilen aus und hinterlässt neben dem Makefile eine diag-Binärdatei. Eine Meldung ld: warning: ignoring duplicate libraries: '-lomp' ist harmlos.

Intel-macOS: Homebrews Präfix ist /usr/local statt /opt/homebrew, ersetze es daher in beiden vorangehenden Pfaden.

Linux (Ubuntu/Debian) — ungetestet

Die folgenden Voraussetzungen decken MPI, OpenMP und BLAS/LAPACK ab. GCC unterstützt -fopenmp nativ, daher ist keine Überschreibung des Compilers nötig. Das wurde unter Linux nicht verifiziert:

sudo apt install -y libopenmpi-dev libomp-dev libblas-dev liblapack-dev

SBD_PATH=../..
CCCOM=mpicxx
CCFLAGS= -std=c++17 -fopenmp -O3
SYSLIB= -llapack -lblas

SqDRIFT ausführen, um CI-Strings zu erzeugen

Wenn du ./SqDRIFT bereits durch die Schritte 1–4 ausgeführt hast, existiert alphadets_from_sqd.txt schon im Projektstammverzeichnis und du kannst diesen Schritt überspringen. Führe es andernfalls im Projektstammverzeichnis aus:

./SqDRIFT
# Writes alphadets_from_sqd.txt in the project root

Die Diagonalisierung ausführen

cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization

# The FCIDump lives in the shared sqdrift/fcidump_files/ directory, one level
# above this tutorial's project root, hence five levels up from here.
ln -sf ../../../../../fcidump_files/N2_sto_3g fcidump.txt

# alphadets_from_sqd.txt is written by ./SqDRIFT into the project root itself,
# which is four levels up.
ln -sf ../../../../alphadets_from_sqd.txt alphadets.txt

./diag \
--fcidump fcidump.txt \
--adetfile alphadets.txt \
--method 0 \
--iteration 100 \
--block 10 \
--tolerance 1e-8 \
| tee ../../../../sbd_output.txt

cd ../../../..

Parameter:

• --fcidump: Pfad zur FCIDump-Datei, die den Hamiltonian enthält
• --adetfile: Pfad zur Datei mit den Alpha-Determinanten (CI-Strings aus SqDRIFT)
• --method 0: Davidson-Diagonalisierungsmethode
• --iteration 100: Maximale Anzahl an Davidson-Iterationen
• --block 10: Blockgröße für den Davidson-Algorithmus
• --tolerance 1e-8: Konvergenztoleranz für die Energie

 Elapsed time for helper construction 0.001512 (sec)
 Elapsed time for init 1e-06 (sec)
 Davidson iteration 0.0 (tol=0.1857493937033264): -107.493531425221
 Davidson iteration 0.1 (tol=0.0455699976495391): -107.5098150493379 -105.93368397386
 Davidson iteration 0.2 (tol=0.009291923190816917): -107.5108356028703 -106.1245977311231 -105.6823366883346
 Davidson iteration 0.3 (tol=0.001488674426419356): -107.5108711084406 -106.3219246016818 -105.7376912083761 -105.4392974334491
 Davidson iteration 0.4 (tol=0.0003188697654737439): -107.5108722237947 -106.5935152393201 -106.1108017566682 -105.4932621297472
 Davidson iteration 0.5 (tol=3.080273116687244e-05): -107.5108722587463 -106.7245930484005 -106.2313070085243 -105.6423237790271
 Davidson iteration 0.6 (tol=2.981303155550973e-06): -107.5108722590151 -106.8008116251279 -106.330411722324 -105.69489630447
 Davidson iteration 0.7 (tol=4.565075388900609e-07): -107.5108722590186 -106.8090569201597 -106.3345525819528 -106.0578325555555
 Davidson iteration 0.8 (tol=4.870689101320012e-08): -107.5108722590186 -106.8098939189322 -106.3973756614146 -106.1617813840531
 Davidson iteration 0.9 (tol=6.224330513865847e-09): -107.5108722590187 -106.8159869231416 -106.4315891597535 -106.2222021964013
 Elapsed time for davidson 0.031276 (sec)
 Elapsed time for diagonalization 0.031286 (sec)
 Elapsed time for mult 0.002511 (sec)
 Energy = -107.5108722590187
 Elapsed time for measurement 4.6e-05 (sec)
 Sample-based diagonalization: Energy = -107.5108722590187
 Sample-based diagonalization: density = [1.999989340109303,1.99999662423008,1.999199121034164,1.981735389885314,1.999825046204946,1.999847567448985,1.995347419300126,0.01136512984757364,0.009844780833945199,0.002849581105572152
 Sample-based diagonalization: carryover bitstrings = [], size = 0

Die Zeile Energy = -107.5108722590187 ist die Grundzustandsabschätzung über dem gesampelten Unterraum.

Wie sich das vom Python-Beispiel unterscheidet

Die beiden Beispiele unterscheiden sich gleichzeitig in mehreren Punkten, nicht nur im Umgang mit Bitstrings mit falscher Symmetrie, daher ist nicht zu erwarten, dass sie identische Energien liefern:

	C++-Tab	Python-Tab
Bitstrings mit falscher Symmetrie	Verworfen durch postselect_bitstrings anhand des Hamming-Gewichts	Repariert durch iterative Konfigurationswiederherstellung
Schleifenstruktur	Einzelner Durchlauf: samplen → postselektieren → einmal diagonalisieren	Äußere Schleife: diagonalisieren, Orbitalbesetzungen auslesen, Konfigurationen wiederherstellen, erneut diagonalisieren
Anzahl der Circuits	100 Circuits (num_circuits = 100)	400 Circuits (200 pro Evolutionszeit)
Evolutionszeiten	Ein einzelnes $\tau$ (time_step = 1)	Zwei Zeiten, times = [1.0, 10.0]
Diagonale Terme	Im gesampelten Operator beibehalten	Durch filter_diagonal_terms entfernt, was die Sampling-Verteilung und die Circuits verändert
Termordnung	Nur gruppiert	canonical_order vor dem Gruppieren angewendet, für Reproduzierbarkeit
Diagonalisierungs-Solver	SBD (diag, externe MPI/OpenMP-Binärdatei)	qiskit-addon-sqd-Solver im selben Prozess
Backend	Fest per Name (backend_name = "ibm_fez"; der Lauf auf dieser Seite verwendete ibm_pittsburgh)	least_busy betriebsbereites Backend mit genügend Qubits
Transpilation	Jeder Circuit separat, Optimierungsstufe 2	Alle Circuits zusammen, Optimierungsstufe 3
Job-Einreichung	100 Jobs, je ein Circuit	Ein Sampler-Job mit allen 400 Circuits
Simulatorlauf	Keiner; nur Hardware	Rauschfreies Simulatorbeispiel vor dem Hardware-Lauf

Beide Versionen verwenden denselben Hamiltonian, gruppieren seine Terme nach elektronischer Struktur, gewichten jede Gruppe mit dem mittleren Absolutwert der Koeffizienten ihrer Terme, starten vom Hartree-Fock-Zustand und samplen 10 Gruppen und 100 Shots pro Circuit.

Jeder Unterschied beeinflusst den gesampelten Unterraum und damit die Energie unabhängig von den anderen. Das größere Ensemble und die zweite, längere Evolutionszeit weiten beide die Konfigurationen, die das Python-Beispiel erkundet; das Sampling zu mehr als einer Zeit verleiht SqDRIFT seine Krylov-artige Abdeckung des niederenergetischen Raums. Das Filtern diagonaler Terme ändert, welche Operatoren der Sampler ziehen kann, und deren Gewichte, sodass die beiden Läufe unterschiedliche Verteilungen samplen.

Die Postselektion ist ein strikter Filter, und ein Shot mit der falschen Teilchenzahl wird gänzlich verworfen. Im obigen Lauf verwarf das den Großteil der Daten: 1137 von 10000 Shots überstanden die (7,7)-Hamming-Gewicht-Prüfung (etwa 11 %) und fielen auf nur 49 verschiedene $\alpha$-Determinanten zusammen. Die Konfigurationswiederherstellung nutzt stattdessen mittlere Orbitalbesetzungen aus einer vorherigen Diagonalisierung, um Bits umzukehren und solche Shots zu reparieren, und verwertet so, was die Postselektion verwerfen würde. Diese Rückkopplung macht sie von Natur aus iterativ: recover_configurations in qiskit-addon-sqd-hpc verlangt ein Argument avg_occupancies, das erst existiert, nachdem eine Diagonalisierung gelaufen ist.

Dieser C++-Lauf liefert einen Unterraum aus 49 $\alpha$-Determinanten und eine Energie von -107.5109 Ha. Der vollständige $\alpha$-Determinantenraum für sieben Elektronen in 10 Orbitalen umfasst nur $\binom{10}{7} = 120$ Strings, sodass die Übergabe aller 120 an dieselbe diag-Binärdatei die exakte FCI-Energie -107.6482 Ha ergibt. Der Lauf liegt somit etwa 0.1373 Ha über dem exakten Wert und deckt dabei nur 49 der 120 $\alpha$-Determinanten ab. Da SQD im ausgewählten Unterraum variationell ist, kann ein kleinerer Unterraum die Energie nur erhöhen, daher ist der unvollständige Unterraum die direkteste Erklärung für die Lücke.

Diese Prüfung gegen die exakte FCI-Energie ist nur möglich, weil das Beispiel klein ist; sie ist nicht Teil des SqDRIFT-Workflows, und für Systemgrößen, bei denen SQD tatsächlich benötigt wird, existiert keine solche Referenz. Sie dient hier nur dazu, zu quantifizieren, wie viel des Raums dieser Lauf erfasst hat.

Ein Vergleich mit dem Python-Beispiel ist nicht kontrolliert, daher sollte ein Unterschied nicht allein dem Fehlen der Konfigurationswiederherstellung zugeschrieben werden; beachte, dass sich Circuit-Anzahl, die einzelne Evolutionszeit, die Filterung diagonaler Terme und das Hardware-Rauschen alle gleichzeitig unterscheiden und alle den Unterraum prägen.

Um die Lücke zu schließen, findest du hier einige mögliche nächste Schritte nach Abschluss dieses Tutorials: erhöhe shots, num_circuits oder ops_per_circuit; samplee zu mehr als einem time_step; und füge auf Qiskit::addon::sqd::recover_configurations aufbauend eine Wiederherstellungsschleife hinzu.

Nächste Schritte

Empfehlungen

• Sample-basierte Krylov-Quantendiagonalisierung eines fermionischen Gittermodells — ein verwandtes Tutorial, das Zeitentwicklungs-Circuits statt eines variationellen Ansatzes verwendet.

• Sample-basierte Quantendiagonalisierung eines chemischen Hamiltonians — ein Tutorial dazu, wie man einen Local-Unitary-Cluster-Jastrow-(LUCJ-)Circuit für die Simulation in der Quantenchemie konstruiert.

• Die SqDRIFT-Veröffentlichung — die Literatur, auf der dieses Tutorial basiert. (Beachte, dass einige der in dieser Veröffentlichung diskutierten Optimierungen derzeit in Arbeit sind und dieses Tutorial sich in Zukunft aufgrund der Weiterentwicklung der verwendeten Bibliotheken ändern kann.)