SQD zur Energieschätzung eines Chemie-Hamiltonoperators

In dieser Lektion wenden wir SQD an, um die Grundzustandsenergie eines Moleküls zu schätzen.

Dabei werden wir die folgenden Themen mithilfe des $4$-stufigen Qiskit-Musteransatzes diskutieren:

1. Schritt 1: Problem auf Quantenschaltkreise und Operatoren abbilden
   • Den molekularen Hamiltonoperator für $N_2$ aufstellen.
   • Den chemisch inspirierten und hardwarefreundlichen Local Unitary Cluster Jastrow (LUCJ) [1] erklären
2. Schritt 2: Für die Zielhardware optimieren
   • Gate-Anzahl und Layout des Ansatzes für die Hardwareausführung optimieren
3. Schritt 3: Auf der Zielhardware ausführen
   • Den optimierten Schaltkreis auf einem echten QPU ausführen, um Stichproben des Unterraums zu erzeugen
4. Schritt 4: Ergebnisse nachverarbeiten
   • Die selbstkonsistente Konfigurationswiederherstellungsschleife [2] einführen
     • Den vollständigen Satz an Bitstring-Stichproben nachverarbeiten, unter Verwendung von Vorwissen über die Teilchenzahl und die mittlere Orbitalbesetzung aus der jüngsten Iteration
     • Probabilistisch Chargen von Teilstichproben aus wiederhergestellten Bitstrings erstellen
     • Den molekularen Hamiltonoperator über jeden abgetasteten Unterraum projizieren und diagonalisieren
     • Die niedrigste Grundzustandsenergie über alle Chargen hinweg speichern und die mittlere Orbitalbesetzung aktualisieren

Wir verwenden in dieser Lektion mehrere Softwarepakete.

• PySCF zum Definieren des Moleküls und Aufstellen des Hamiltonoperators
• ffsim-Paket zum Konstruieren des LUCJ-Ansatzes
• Qiskit zum Transpilieren des Ansatzes für die Hardwareausführung
• Qiskit IBM Runtime zum Ausführen des Schaltkreises auf einem QPU und Sammeln von Stichproben
• Qiskit addon SQD für Konfigurationswiederherstellung und Grundzustandsenergieschätzung mittels Unterraumprojektion und Matrixdiagonalisierung

1. Problem auf Quantenschaltkreise und Operatoren abbilden

Molekularer Hamiltonoperator

Ein molekularer Hamiltonoperator hat die allgemeine Form:

$$\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \sum_{ \substack{prqs\\\sigma\tau} } \frac{(pr|qs)}{2} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma}$$

$\hat{a}^\dagger_{p\sigma}$/$\hat{a}_{p\sigma}$ sind die fermionischen Erzeugungs-/Vernichtungsoperatoren, die dem $p$-ten Basissatzelement und dem Spin $\sigma$ zugeordnet sind. $h_{pr}$ und $(pr|qs)$ sind die Einzel- und Zweikörper-Elektronenintegrale. Mit PySCF definieren wir das Molekül und berechnen die Ein- und Zweikörperintegrale des Hamiltonoperators für den Basissatz 6-31g.

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

import warnings
import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# Specify molecule properties
open_shell = False
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],  # Two N atoms 1 angstrom apart
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)  # hcore: one-body integrals
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)  # eri: two-body integrals

# Compute exact energy for comparison
exact_energy = cas.run().e_tot

converged SCF energy = -108.835236570774
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

In dieser Lektion verwenden wir die Jordan-Wigner (JW)-Transformation, um eine fermionische Wellenfunktion auf eine Qubit-Wellenfunktion abzubilden, damit sie mit einem Quantenschaltkreis präpariert werden kann. Die JW-Transformation bildet den Fock-Raum von Fermionen in M räumlichen Orbitalen auf den Hilbertraum von 2M Qubits ab – das heißt, ein räumliches Orbital wird in zwei Spin-Orbitale aufgeteilt: eines für ein Spin-auf ($\alpha$)-Elektron und eines für Spin-ab ($\beta$). Ein Spin-Orbital kann besetzt oder unbesetzt sein. Wenn wir von der Anzahl der Orbitale sprechen, meinen wir in der Regel die Anzahl der räumlichen Orbitale; die Anzahl der Spin-Orbitale ist doppelt so groß. In Quantenschaltkreisen repräsentieren wir jedes Spin-Orbital durch ein Qubit. Somit stellt eine Gruppe von Qubits Spin-auf- oder $\alpha$-Orbitale dar, und eine andere Gruppe repräsentiert Spin-ab- oder $\beta$-Orbitale. Das $N_2$-Molekül hat zum Beispiel für den Basissatz 6-31g $16$ räumliche Orbitale (also $16$ $\alpha$ + $16$ $\beta$ = $32$ Spin-Orbitale). Daher benötigen wir einen $32$-Qubit-Quantenschaltkreis (eventuell werden zusätzliche Hilfs-Qubits benötigt, wie später besprochen). Die Qubits werden in der Rechenbasis gemessen, um Bitstrings zu erzeugen, die elektronische Konfigurationen oder (Slater-)Determinanten darstellen. In dieser Lektion verwenden wir die Begriffe Bitstrings, Konfigurationen und Determinanten synonym. Die Bitstrings geben die Elektronenbesetzung in Spin-Orbitalen an: Eine $1$ an einer Bitposition bedeutet, dass das entsprechende Spin-Orbital besetzt ist, während eine $0$ bedeutet, dass es leer ist. Da Elektronenstrukturprobleme teilchenzahlerhaltend sind, muss eine feste Anzahl von Spin-Orbitalen besetzt sein. Das $N_2$-Molekül hat $5$ Spin-auf ($\alpha$)- und $5$ Spin-ab ($\beta$)-Elektronen. Daher muss jeder Bitstring, der die $\alpha$- und $\beta$-Orbitale darstellt, für das $N_2$-Molekül jeweils fünf $1\text{en}$ haben.

1.1 Quantenschaltkreis zur Stichprobenerzeugung: Der LUCJ-Ansatz

In dieser Lektion verwenden wir den Local Unitary Coupled Cluster Jastrow (LUCJ) \[1\]-Ansatz für die Quantenzustandspräparation und anschließende Stichprobennahme. Zunächst erklären wir die verschiedenen Bausteine des vollständigen UCJ-Ansatzes und die beim lokalen Ansatz vorgenommenen Näherungen. Anschließend konstruieren wir mit dem ffsim-Paket den LUCJ-Ansatz und optimieren ihn mit dem Qiskit Transpiler für die Hardwareausführung.

Der UCJ-Ansatz hat folgende Form (für ein Produkt aus $L$ Schichten oder Wiederholungen des UCJ-Operators):

$$|\psi\rangle = \prod_{\mu=1}^{L}{(e^{K^{\mu}} \times {e^{iJ^{\mu}}} \times {e^{-K^{\mu}}})} |\Phi_{0}\rangle$$

wobei $\vert \Phi_{0} \rangle$ ein Referenzzustand ist, typischerweise der Hartree-Fock (HF)-Zustand. Da der Hartree-Fock-Zustand so definiert ist, dass die niedrigsten Orbitale besetzt sind, beinhaltet die HF-Zustandspräparation das Anwenden von X-Gates, um die Qubits der besetzten Orbitale auf eins zu setzen. Zum Beispiel könnte der HF-Zustandspräparationsblock für 4 räumliche Orbitale und je 2 Spin-auf- und Spin-ab-Elektronen wie folgt aussehen: Eine einzelne Wiederholung des UCJ-Operators ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ besteht aus einer diagonalen Coulomb-Entwicklung ($e^{iJ^{(\mu)}}$), eingebettet zwischen Orbitalrotationen ($e^{K^{(\mu)}}$ und $e^{-K^{(\mu)}}$). Orbitalrotationsblöcke wirken auf eine einzelne Spinspezies ($\alpha$ (Spin-auf)/$\beta$ (Spin-ab)). Für jede Elektronenspezies besteht die Orbitalrotation aus einer Schicht von Einzel-Qubit-$R_{z}$-Gates, gefolgt von einer Sequenz aus Zwei-Qubit-Givens-Rotationsgates ($XX + YY$-Gates).

Die Zwei-Qubit-Gates wirken auf benachbarte Spin-Orbitale (nächste Nachbar-Qubits) und sind daher auf IBM®-QPUs ohne SWAP-Gates implementierbar. Der $e^{iJ^{(\mu)}}$-Operator, auch bekannt als diagonaler Coulomb-Operator, besteht aus drei Blöcken. Zwei davon wirken auf gleiche Spinbereiche ($e^{iJ_{\alpha \alpha}^{(\mu)}}$ und $e^{iJ_{\beta \beta}^{(\mu)}}$), und einer wirkt zwischen zwei Spinbereichen ($e^{iJ_{\alpha \beta}^{(\mu)}}$).

Alle Blöcke in $e^{iJ^{(\mu)}}$ bestehen aus Zahl-Zahl-Gates $U_{nn}(\phi)$ [1]. Ein $U_{nn}(\phi)$-Gate lässt sich weiter in ein $R_{ZZ}(\frac{\phi}{2})$-Gate, gefolgt von zwei Einzel-Qubit-$Rz(-\frac{\phi}{2})$-Gates auf zwei separaten Qubits, zerlegen.

Gleichspinige Komponenten ($J_{\alpha \alpha}$ und $J_{\beta \beta}$) haben $U_{nn}$-Gates zwischen allen möglichen Qubit-Paaren. Da supraleitende QPUs jedoch eine eingeschränkte Konnektivität haben, müssen Qubits getauscht (SWAPped) werden, um Gates zwischen nicht benachbarten Qubits zu realisieren.

Betrachten wir zum Beispiel den folgenden $e^{iJ_{\alpha \alpha}^{(\mu)}}$- (bzw. $e^{iJ_{\beta \beta}^{(\mu)}}$-) Block für $N = 4$ räumliche Orbitale. Bei einer linearen Qubit-Konnektivität sind die letzten drei Gates nicht direkt implementierbar, da sie zwischen nicht benachbarten Qubits wirken (zum Beispiel sind Q0 und Q2 nicht direkt verbunden). Daher benötigen wir SWAP-Gates, um sie benachbart zu machen (die folgende Abbildung zeigt ein Beispiel mit $3$ SWAP-Gates). Als Nächstes implementiert der $J_{\alpha \beta}$-Block Gates zwischen gleichindizierten Orbitalen aus verschiedenen Spinbereichen (zum Beispiel zwischen $0\alpha$ und $0\beta$). Ebenso erfordern diese Gates SWAPs, wenn die Qubits auf einem QPU physisch nicht benachbart sind. Aus der obigen Diskussion ergibt sich, dass der UCJ-Ansatz bei der HW-Ausführung auf Hindernisse stößt, da er aufgrund nicht benachbarter Qubit-Wechselwirkungen SWAP-Gates benötigt. Die lokale Variante des UCJ-Ansatzes, LUCJ, begegnet dieser Herausforderung, indem einige $U_{nn}$-Gates aus dem diagonalen Coulomb-Operator entfernt werden.

In den gleichspinigen Blöcken ($J_{\alpha \alpha}$ und $J_{\beta \beta}$) behalten wir in der LUCJ-Version nur die $U_{nn}$-Gates, die mit Nächste-Nachbar-Konnektivität kompatibel sind, und entfernen Gates zwischen nicht benachbarten Qubits. Die folgende Abbildung zeigt den LUCJ-Block nach Entfernung der nicht benachbarten Gates. Als Nächstes kann die LUCJ-Version des $J_{\alpha \beta}$-Blocks, der zwischen verschiedenen Elektronenspezies wirkt, je nach Gerätetopologie unterschiedliche Formen annehmen.

Auch hier entfernt die LUCJ-Version inkompatible Gates. Die folgende Abbildung zeigt Varianten des $J_{\alpha \beta}$-Blocks für verschiedene Qubit-Topologien, darunter Gitter, hexagonal, Heavy-Hex und linear.

• Quadratisch: Wir können $U_{nn}$-Gates zwischen allen $\alpha$- und $\beta$-Orbitalen ohne SWAPs haben und müssen daher keine $U_{nn}$-Gates entfernen.
• Heavy-Hex: Die $\alpha$-$\beta$-Wechselwirkungen werden zwischen jedem $4$-ten indizierten Spin-Orbital (zum Beispiel dem 0., 4. und 8.) beibehalten und sind hilfsqubit-vermittelt, das heißt, wir benötigen Hilfs-Qubits zwischen den linearen Ketten, die $\alpha$- und $\beta$-Orbitale darstellen. Diese Anordnung benötigt eine begrenzte Anzahl von SWAPs.
• Hexagonal: Jedes zweite Orbital, zum Beispiel das 0., 2. und 4. indizierte Orbital, wird zu Nächsten Nachbarn, wenn $\alpha$ und $\beta$ in zwei benachbarten linearen Ketten angeordnet sind.
• Linear: Nur ein $\alpha$- und ein $\beta$-Orbital sind verbunden, was bedeutet, dass der $J_{\alpha \beta}$-Block nur ein Gate hat. Das Entfernen von Gates aus dem UCJ-Ansatz zur Erstellung der LUCJ-Version macht ihn zwar HW-kompatibler, aber der Ansatz verliert etwas Ausdrucksstärke. Daher können mehr Wiederholungen ($L$) des modifizierten UCJ-Operators beim LUCJ-Ansatz notwendig sein.

1.2 Initialisierung des LUCJ-Ansatzes

Der LUCJ ist ein parametrisierter Ansatz, und wir müssen die Parameter vor der Hardwareausführung initialisieren. Eine Möglichkeit zur Initialisierung des Ansatzes ist die Verwendung von t1- und t2-Amplituden aus der klassischen Coupled Cluster Singles and Doubles (CCSD)-Methode, wobei t1-Amplituden die Koeffizienten der Einzel-Anregungsoperatoren und t2-Amplituden die der Doppel-Anregungsoperatoren sind.

Beachte, dass die Initialisierung des LUCJ-Ansatzes mit t1- und t2-Amplituden zwar anständige Ergebnisse liefert, die Ansatzparameter jedoch möglicherweise weiterer Optimierung bedürfen.

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

t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 Konstruktion des LUCJ-Ansatzes mit ffsim

Wir verwenden das ffsim-Paket, um den Ansatz zu erstellen und mit den oben berechneten t1- und t2-Amplituden zu initialisieren. Da unser Molekül einen geschlossenschaligen Hartree-Fock-Zustand hat, verwenden wir die spinausgewogene Variante des UCJ-Ansatzes, UCJOpSpinBalanced.

Da IBM-Hardware eine Heavy-Hex-Topologie hat, übernehmen wir das in [1] verwendete und oben erläuterte Zickzack-Muster für Qubit-Wechselwirkungen. In diesem Muster sind Orbitale (Qubits) mit gleichem Spin mit einer Linientopologie verbunden (rote und blaue Kreise). Aufgrund der Heavy-Hex-Topologie haben Orbitale unterschiedlicher Spins Verbindungen zwischen jedem 4. Orbital, also dem 0., 4., 8. und so weiter (lila Kreise).

import ffsim
from qiskit import QuantumCircuit, QuantumRegister

n_reps = 2
alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]
alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]

ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),
)

nelec = (num_elec_a, num_elec_b)

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

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

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()
# circuit.decompose().draw("mpl", scale=0.5, fold=-1)

Der LUCJ-Ansatz mit wiederholten Schichten kann durch Zusammenführen benachbarter Blöcke optimiert werden. Betrachte den Fall n_reps=2: Die beiden Orbitalrotationsblöcke in der Mitte können zu einem einzigen Orbitalrotationsblock zusammengeführt werden. Das ffsim-Paket hat einen Pass Manager namens ffsim.qiskit.PRE_INIT, um den Schaltkreis durch Zusammenführen solcher benachbarten Blöcke zu optimieren.

2. Für die Zielhardware optimieren

Zunächst rufen wir ein Backend unserer Wahl ab. Wir optimieren unseren Schaltkreis für das Backend und führen dann den optimierten Schaltkreis auf demselben Backend aus, um Stichproben für den Unterraum zu erzeugen.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

Als Nächstes empfehlen wir die folgenden Schritte, um den Ansatz zu optimieren und hardwarekompatibel zu machen.

• Physikalische Qubits (initial_layout) von der Zielhardware auswählen, die dem oben beschriebenen Zickzack-Muster (zwei lineare Ketten mit Hilfs-Qubit dazwischen) entsprechen. Das Layout der Qubits in diesem Muster führt zu einem effizienten, hardwarekompatiblen Schaltkreis mit weniger Gates.
• Einen gestuften Pass Manager mit der Funktion generate_preset_pass_manager aus Qiskit mit dem gewählten backend und initial_layout erstellen.
• Die pre_init-Stufe des gestuften Pass Managers auf ffsim.qiskit.PRE_INIT setzen. ffsim.qiskit.PRE_INIT enthält Qiskit-Transpiler-Passes, die Gates in Orbitalrotationen zerlegen und diese dann zusammenführen, was zu weniger Gates im endgültigen Schaltkreis führt.
• Den Pass Manager auf deinen Schaltkreis anwenden.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]
spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]

initial_layout = spin_a_layout + spin_b_layout

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend, initial_layout=initial_layout
)

# without PRE_INIT passes
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}")

# with PRE_INIT passes
# We will use the circuit generated by this pass manager for hardware execution
pass_manager.pre_init = ffsim.qiskit.PRE_INIT
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}")

Gate counts (w/o pre-init passes): OrderedDict({'rz': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})

3. Auf der Zielhardware ausführen

Nach der Optimierung des Schaltkreises für die Hardwareausführung sind wir bereit, ihn auf der Zielhardware auszuführen und Stichproben für die Grundzustandsenergieschätzung zu sammeln. Da wir nur einen Schaltkreis haben, verwenden wir den Job-Ausführungsmodus von Qiskit Runtime und führen unseren Schaltkreis aus.

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
sampler.options.dynamical_decoupling.enable = True

job = sampler.run([isa_circuit], shots=10_000)  # Takes approximately 5sec of QPU time

# Run cell after IQX job completion
primitive_result = job.result()
pub_result = primitive_result[0]
counts = pub_result.data.meas.get_counts()

4. Ergebnisse nachverarbeiten

Der Nachverarbeitungsteil des SQD-Workflows lässt sich mit folgendem Diagramm zusammenfassen.

Das Abtasten des LUCJ-Ansatzes in der Rechenbasis erzeugt einen Pool verrauschter Konfigurationen $\tilde{\mathcal{\chi}}$, die in der Nachverarbeitungsroutine verwendet werden. Diese umfasst eine Methode namens Konfigurationswiederherstellung (Details folgen), um Konfigurationen mit falscher Elektronenzahl probabilistisch zu korrigieren. Konfigurationen mit korrekter Elektronenzahl $\tilde{\mathcal{\chi}}_{R}$ werden dann in mehrere Chargen aufgeteilt und auf Basis der Häufigkeit jeder eindeutigen Konfiguration verteilt. Jede Charge von Stichproben definiert einen Unterraum ($\mathcal{S^{(k)}}$). Anschließend wird der molekulare Hamiltonoperator $H$ auf die Unterräume projiziert:

$$H_{\mathcal{S}^{(k)}} = P_{\mathcal{S}^{(k)}} H _{\mathcal{S}^{(k)}} \text{ with } P_{\mathcal{S}^{(k)}} = \sum_{x \in \mathcal{S}^{(k)}} \vert x \rangle \langle x \vert$$

Jeder projizierte Hamiltonoperator $H_{\mathcal{S}^{(k)}}$ wird dann an einen Eigenwertlöser übergeben, der ihn diagonalisiert, um Eigenwerte und Eigenvektoren zur Rekonstruktion eines Eigenzustands zu berechnen. In dieser Lektion projizieren und diagonalisieren wir den Hamiltonoperator mit dem qiskit-addon-sqd-Paket, das die Davidson-Methode aus PySCF zur Diagonalisierung verwendet.

$$H_{\mathcal{S}^{(k)}} \vert \psi^{(k)} \rangle = E^{(k)} \vert \psi^{(k)} \rangle$$

Wir sammeln dann den niedrigsten Eigenwert (Energie) aus den Chargen und berechnen auch die mittlere Orbitalbesetzung $\text{n}$. Die Information über die mittlere Besetzung wird im Konfigurationswiederherstellungsschritt verwendet, um verrauschte Konfigurationen probabilistisch zu korrigieren.

Als Nächstes erklären wir die selbstkonsistente Konfigurationswiederherstellungsschleife im Detail und zeigen konkrete Codebeispiele zur Implementierung der oben genannten Schritte zur Schätzung der Grundzustandsenergie des $N_2$-Hamiltonoperators.

4.1 Konfigurationswiederherstellung: Überblick

Jedes Bit in einem Bitstring (Slater-Determinante) repräsentiert ein Spin-Orbital. Die rechte Hälfte eines Bitstrings repräsentiert Spin-auf-Orbitale, und die linke Hälfte repräsentiert Spin-ab-Orbitale. Eine 1 bedeutet, dass das Orbital von einem Elektron besetzt ist, und eine 0 bedeutet, dass das Orbital leer ist. Wir kennen die korrekte Teilchenzahl (sowohl Spin-auf- als auch Spin-ab-Elektronen) a priori. Angenommen, wir haben eine Determinante $x$ mit $N_x$ Elektronen (also $N_x$ Einsen im Bitstring). Die korrekte Teilchenzahl ist $N$. Wenn $N_x \neq N$, wissen wir, dass der Bitstring durch Rauschen verfälscht wurde. Die selbstkonsistente Konfigurationsroutine versucht, den Bitstring zu korrigieren, indem sie probabilistisch $|N_x - N|$ Bits unter Nutzung von Informationen zur mittleren Orbitalbesetzung umkippt. Die mittlere Orbitalbesetzung ($n$) gibt an, wie wahrscheinlich es ist, dass ein Orbital von einem Elektron besetzt wird. Wenn $N_x < N$ ist, haben wir zu wenige Elektronen und müssen einige $0$en in $1$en umkippen, und umgekehrt.

Die Flip-Wahrscheinlichkeit kann $|x[i] - avg\_occupancy[i]|$ für das $i$-te Spin-Orbital sein. In [2] verwendeten die Autoren eine gewichtete Flip-Wahrscheinlichkeit mit einer modifizierten ReLU-Funktion.

$$\begin{align} w(y) = \begin{cases} \delta \frac{y}{h} & \text{if } y \leq h\\ \nonumber \delta + (1 - \delta) \frac{y - h}{1 - h} & \text{if } y > h \end{cases} \end{align}$$

Hierbei definiert $h$ die Position der „Ecke" der ReLU-Funktion, und der Parameter $\delta$ definiert den Wert der ReLU-Funktion an der Ecke. Für $\delta = 0$ wird $w$ zur echten ReLU-Funktion, und für $\delta > 0$ zur modifizierten ReLU. In der Publikation verwendeten die Autoren $\delta = 0.01$ und $h =$ Anzahl der Alpha- (oder Beta-)Teilchen / Anzahl der Alpha- (oder Beta-)Spin-Orbitale $= N/M$ (Füllfaktor).

Die mittlere Orbitalbesetzung ($n$) ist a priori nicht bekannt. Die erste Iteration der Grundzustandsschätzung beginnt mit Konfigurationen, die nur korrekte Teilchenzahlen in beiden Spinspezies haben. Nach der ersten Iteration haben wir eine Schätzung des Grundzustands, und anhand dieser können wir den ersten Schätzwert für $n$ konstruieren. Dieser Schätzwert von $n$ wird verwendet, um Konfigurationen wiederherzustellen, die nächste Iteration der Grundzustandsschätzung durchzuführen und den Schätzwert von $n$ selbstkonsistent zu verfeinern. Der Prozess wiederholt sich, bis ein Abbruchkriterium erfüllt ist.

Betrachte folgendes Beispiel für $N = 2$ und $x = |1000\rangle$ ($N_x = 1$). Wir müssen eine der Nullen auf eins kippen, um die Teilchenzahl zu korrigieren, und die Möglichkeiten sind 1100, 1010 und 1001. Basierend auf der Flip-Wahrscheinlichkeit wird eine der Möglichkeiten als wiederhergestellte Konfiguration (oder Bitstring mit korrekter Teilchenzahl) ausgewählt.

Angenommen, wir führen in der ersten Iteration zwei Chargen durch, und die geschätzten Grundzustände aus ihnen sind:

$$\begin{align}\nonumber \text{Batch0: } \vert \psi \rangle &= 0.8 \times \vert 1001 \rangle + 0.6 \times \vert 0110 \rangle \\ \nonumber \text{Batch1: } \vert \psi \rangle &= \frac{1}{\sqrt{3}} \left( \vert 1001 \rangle + \vert 0101 \rangle + \vert 0110 \rangle \right) \nonumber \end{align}$$

Anhand der Rechenbasis-Zustände und ihrer Amplituden können wir die Wahrscheinlichkeit der Elektronenbesetzung (kurz Besetzung) pro Spin-Orbital (Qubit) berechnen (beachte: Wahrscheinlichkeit = |Amplitude|$^2$). Im Folgenden tabellieren wir qubitweise Besetzungen für jeden Bitstring im geschätzten Grundzustand und berechnen die Gesamtorbitalbesetzung für eine Charge. Beachte, dass das am weitesten rechts stehende Bit gemäß der Qiskit-Konvention Qubit-0 (Q0) und das am weitesten links stehende Bit Q3 darstellt.

Besetzung (Batch0):

	Q3	Q2	Q1	Q0
1001	0.64	0.0	0.0	0.64
0110	0.0	0.36	0.36	0.0
n (Batch0)	0.64	0.36	0.36	0.64

Besetzung (Batch1)

	Q3	Q2	Q1	Q0
1001	0.33	0.00	0.00	0.33
0101	0.0	0.33	0.00	0.33
0110	0.0	0.33	0.33	0.00
n (Batch1)	0.33	0.66	0.33	0.66

Besetzung (Durchschnitt über alle Chargen)

	Q3	Q2	Q1	Q0
n (Batch0)	0.64	0.36	0.36	0.64
n (Batch1)	0.33	0.66	0.33	0.66
n (average)	0.49	0.51	0.35	0.65

Anhand der oben berechneten mittleren Orbitalbesetzung können wir die Flip-Wahrscheinlichkeiten für verschiedene Orbitale in der Konfiguration $x = \vert 1000 \rangle$ bestimmen. Da das von Q3 repräsentierte Orbital bereits besetzt ist und nicht umgekippt werden muss, setzen wir seine p(flip) auf $0$. Für die verbleibenden unbesetzten Orbitale ist die Flip-Wahrscheinlichkeit jeweils $\vert x[i] - \text{n}[i] \vert$. Zusammen mit p(flip) berechnen wir auch das Wahrscheinlichkeitsgewicht des Kippens mit der oben beschriebenen modifizierten ReLU-Funktion.

Flip-Wahrscheinlichkeit ($x = \vert 1000 \rangle$, $\delta = 0.01$, $h = N/M = 2/4 = 0.50$)

	Q3	Q2	Q1	Q0
p(flip) ($\vert x[i] - \text{n}[i] \vert$)	0	0.51	0.35	0.65
w(p(flip))	0	0.03	0.007	0.31

Schließlich können wir mit den obigen gewichteten Wahrscheinlichkeiten eines der unbesetzten Orbitale Q2, Q1 und Q0 umkippen. Basierend auf den obigen Werten wird Q0 am wahrscheinlichsten umgekippt, und eine mögliche wiederhergestellte Konfiguration wäre $\vert \text{1001} \rangle$. Der vollständige selbstkonsistente Konfigurationswiederherstellungsprozess lässt sich wie folgt zusammenfassen:

Erste Iteration: Angenommen, die vom Quantencomputer erzeugten Bitstrings (Konfigurationen oder Slater-Determinanten) bilden eine Menge $\widetilde{\chi}$, die sowohl Konfigurationen mit korrekter ($\widetilde{\chi}_{correct}$) als auch inkorrekter ($\widetilde{\chi}_{incorrect}$) Teilchenzahl in jedem Spinbereich enthält.

1. Konfigurationen aus ($\widetilde{\chi}_{correct}$) werden zufällig zu Chargen $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ von Vektoren für die Unterraumprojektion zusammengestellt. Die Anzahl der Chargen und die Stichprobenanzahl in jeder Charge sind benutzerdefinierte Parameter. Je mehr Stichproben in jeder Charge, desto größer die Unterraumdimension und desto rechenaufwendiger wird die Diagonalisierung. Andererseits können zu wenige Stichproben dazu führen, dass Grundzustandsstützvektoren fehlen und eine inkorrekte Schätzung entsteht.
2. Den Eigenzustandslöser (also Projektion auf den Unterraum und Diagonalisierung) auf die Chargen anwenden und näherungsweise Eigenzustände $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$ erhalten.
3. Aus den genäherten Eigenzuständen den ersten Schätzwert für $n$ konstruieren.

Folgende Iterationen:

1. Mit $n$ die Konfigurationen mit falscher Teilchenzahl in $\widetilde{\chi}_{incorrect}$ korrigieren. Diese nennen wir $\widetilde{\chi}_{correct\_new}$. Dann bildet $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ die neue Menge von Konfigurationen mit korrekter Teilchenzahl.
2. $\widetilde{\chi}_{R}$ wird abgetastet, um Chargen $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$ zu erstellen.
3. Der Eigenzustandslöser läuft mit neuen Chargen und erzeugt neue Schätzungen der Grundzustände $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. Aus den genäherten Eigenzuständen den verfeinerten Schätzwert für $n$ konstruieren.
5. Falls das Abbruchkriterium nicht erfüllt ist, zurück zu Schritt 2.1.

4.2 Grundzustandsschätzung

Zunächst transformieren wir die Zähler in eine Bitstring-Matrix und ein Wahrscheinlichkeitsarray für die Nachverarbeitung.

Jede Zeile in der Matrix repräsentiert einen eindeutigen Bitstring. Da Qubits in Qiskit von rechts im Bitstring indiziert werden, repräsentiert Spalte 0 Qubit N-1 und Spalte N-1 repräsentiert Qubit 0, wobei N die Anzahl der Qubits ist.

Die Alpha-Orbitale werden im Spaltenindexbereich (N, N/2] (rechte Hälfte) dargestellt, und die Beta-Orbitale im Bereich (N/2, 0] (linke Hälfte).

from qiskit_addon_sqd.counts import counts_to_arrays

# Convert counts into bitstring and probability arrays
bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)

Es gibt einige benutzerkontrollierbare Optionen, die für diese Methode wichtig sind:

• iterations: Anzahl der selbstkonsistenten Konfigurationswiederherstellungsiterationen
• n_batches: Anzahl der Konfigurationschargen, die von den verschiedenen Aufrufen des Eigenzustandslösers verwendet werden
• samples_per_batch: Anzahl der eindeutigen Konfigurationen in jeder Charge
• max_davidson_cycles: Maximale Anzahl von Davidson-Zyklen, die von jedem Eigenwertlöser ausgeführt werden

import numpy as np
from qiskit_addon_sqd.configuration_recovery import recover_configurations
from qiskit_addon_sqd.fermion import (
    bitstring_matrix_to_ci_strs,
    solve_fermion,
)
from qiskit_addon_sqd.subsampling import postselect_and_subsample

rng = np.random.default_rng(24)
# SQD options
iterations = 5

# Eigenstate solver options
n_batches = 5
samples_per_batch = 500
max_davidson_cycles = 300

# Self-consistent configuration recovery loop
e_hist = np.zeros((iterations, n_batches))  # energy history
s_hist = np.zeros((iterations, n_batches))  # spin history
occupancy_hist = []
avg_occupancy = None
for i in range(iterations):
    print(f"Starting configuration recovery iteration {i}")
    # On the first iteration, we have no orbital occupancy information from the
    # solver, so we begin with the full set of noisy configurations.
    if avg_occupancy is None:
        bs_mat_tmp = bitstring_matrix_full
        probs_arr_tmp = probs_arr_full

    # If we have average orbital occupancy information, we use it to refine the full set of noisy configurations
    else:
        bs_mat_tmp, probs_arr_tmp = recover_configurations(
            bitstring_matrix_full,
            probs_arr_full,
            avg_occupancy,
            num_elec_a,
            num_elec_b,
            rand_seed=rng,
        )

    # Create batches of subsamples. We post-select here to remove configurations
    # with incorrect hamming weight during iteration 0, since no config recovery was performed.
    batches = postselect_and_subsample(
        bs_mat_tmp,
        probs_arr_tmp,
        hamming_right=num_elec_a,
        hamming_left=num_elec_b,
        samples_per_batch=samples_per_batch,
        num_batches=n_batches,
        rand_seed=rng,
    )

    # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.
    e_tmp = np.zeros(n_batches)
    s_tmp = np.zeros(n_batches)
    occs_tmp = []
    coeffs = []
    for j in range(n_batches):
        strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])
        print(f"  Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}")
        energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(
            batches[j],
            hcore,
            eri,
            open_shell=open_shell,
            spin_sq=spin_sq,
            max_davidson=max_davidson_cycles,
        )
        energy_sci += nuclear_repulsion_energy
        e_tmp[j] = energy_sci
        s_tmp[j] = spin
        occs_tmp.append(avg_occs)
        coeffs.append(coeffs_sci)

    # Combine batch results
    avg_occupancy = tuple(np.mean(occs_tmp, axis=0))

    # Track optimization history
    e_hist[i, :] = e_tmp
    s_hist[i, :] = s_tmp
    occupancy_hist.append(avg_occupancy)

Starting configuration recovery iteration 0
  Batch 0 subspace dimension: 21609
  Batch 1 subspace dimension: 21609
  Batch 2 subspace dimension: 21609
  Batch 3 subspace dimension: 21609
  Batch 4 subspace dimension: 21609
Starting configuration recovery iteration 1
  Batch 0 subspace dimension: 609961
  Batch 1 subspace dimension: 616225
  Batch 2 subspace dimension: 627264
  Batch 3 subspace dimension: 633616
  Batch 4 subspace dimension: 624100
Starting configuration recovery iteration 2
  Batch 0 subspace dimension: 564001
  Batch 1 subspace dimension: 605284
  Batch 2 subspace dimension: 582169
  Batch 3 subspace dimension: 559504
  Batch 4 subspace dimension: 591361
Starting configuration recovery iteration 3
  Batch 0 subspace dimension: 550564
  Batch 1 subspace dimension: 549081
  Batch 2 subspace dimension: 531441
  Batch 3 subspace dimension: 527076
  Batch 4 subspace dimension: 531441
Starting configuration recovery iteration 4
  Batch 0 subspace dimension: 544644
  Batch 1 subspace dimension: 580644
  Batch 2 subspace dimension: 527076
  Batch 3 subspace dimension: 531441
  Batch 4 subspace dimension: 537289

4.3 Diskussion der Ergebnisse

Die erste Darstellung zeigt, dass wir nach einigen Iterationen die Grundzustandsenergie auf etwa ~24 mH genau schätzen (die chemische Genauigkeit wird typischerweise auf 1 kcal/mol $\approx$ 1,6 mH festgelegt). Die zweite Darstellung zeigt die mittlere Besetzung jedes räumlichen Orbitals nach der letzten Iteration. Wir sehen, dass sowohl die Spin-auf- als auch die Spin-ab-Elektronen in unseren Lösungen mit hoher Wahrscheinlichkeit die ersten fünf Orbitale besetzen.

Obwohl die geschätzte Grundzustandsenergie angemessen ist, liegt sie nicht im Bereich der chemischen Genauigkeit ($\pm \approx 1,6$ mH). Diese Abweichung lässt sich auf die kleine Unterraumdimension zurückführen, die wir oben für die Projektion und Diagonalisierung verwendet haben. Da wir samples_per_batch=500 verwenden, wird der Unterraum von maximal $500$ Vektoren aufgespannt, dem Grundzustandsstützvektoren fehlen. Eine Erhöhung des Parameters samples_per_batch sollte die Genauigkeit auf Kosten von mehr klassischen Rechenressourcen und Laufzeit verbessern.

# Data for energies plot
x1 = range(iterations)
min_e = [np.min(e) for e in e_hist]
e_diff = [abs(e - exact_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]

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

# Data for avg spatial orbital occupancy
y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]
x2 = range(len(y2))

import matplotlib.pyplot as plt

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

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

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

print(f"Exact energy: {exact_energy:.5f} Ha")
print(f"SQD energy: {min_e[-1]:.5f} Ha")
print(f"Absolute error: {e_diff[-1]:.5f} Ha")
plt.tight_layout()
plt.show()

Exact energy: -109.04667 Ha
SQD energy: -109.02234 Ha
Absolute error: 0.02434 Ha

Übung für den Leser

Erhöhe den Parameter samples_per_batch schrittweise (zum Beispiel von $1000$ auf $10000$ in Schritten von $1000$, soweit es der Arbeitsspeicher deines Computers erlaubt) und vergleiche die geschätzten Grundzustandsenergien.

Referenzen

[1] M. Motta et al., "Bridging physical intuition and hardware efficiency for correlated electronic states: the local unitary cluster Jastrow ansatz for electronic structure" (2023). Chem. Sci., 2023, 14, 11213.

[2] J. Robledo-Moreno et al., "Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer" (2024). arXiv:quant-ph/2405.05068.