SQD zur Energieabschätzung eines chemischen Hamilton-Operators

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

Insbesondere werden wir die folgenden Themen anhand des $4$ -step Qiskit-Musteransatzes diskutieren:

1. Schritt 1: Abbildung des Problems auf Quantenschaltungen und -operatoren
   • Stellen Sie den molekularen Hamiltonian für $N_2$ auf.
   • Erläutern Sie den von der Chemie inspirierten und hardwarefreundlichen lokalen unitären Cluster Jastrow (LUCJ) [1]
2. Schritt 2: Optimieren für die Zielhardware
   • Optimierung der Gatteranzahl und des Layouts des Ansatzes für die Hardware-Ausführung
3. Schritt 3: Ausführen auf der Zielhardware
   • Führen Sie die optimierte Schaltung auf einer echten QPU aus, um Stichproben des Unterraums zu erzeugen.
4. Schritt 4: Nachbearbeitung der Ergebnisse
   • Einführung der Schleife zur Wiederherstellung der selbstkonsistenten Konfiguration [2]
     • Nachbearbeitung des vollständigen Satzes von Bitstring-Stichproben unter Verwendung des Vorwissens über die Partikelanzahl und die durchschnittliche Orbitalbesetzung, die bei der letzten Iteration berechnet wurde.
     • Wahrscheinliche Erstellung von Stapeln von Unterstichproben aus wiederhergestellten Bitstrings.
     • Projizieren und diagonalisieren Sie den molekularen Hamiltonian über jeden abgetasteten Unterraum.
     • Speichern Sie die minimale Grundzustandsenergie, die in allen Stapeln gefunden wurde, und aktualisieren Sie die durchschnittliche Orbitalbesetzung.

Während des Unterrichts werden wir verschiedene Softwarepakete verwenden.

• PySCF um das Molekül zu definieren und den Hamiltonian aufzustellen.
• ffsim paket, um den LUCJ-Ansatz zu konstruieren.
• Qiskit zum Transpilieren des Ansatzes für die Hardware-Ausführung.
• Qiskit IBM Runtime um die Schaltung auf einer QPU auszuführen und Proben zu sammeln.
• Qiskit addon SQD wiederherstellung der Konfiguration und Schätzung der Energie des Grundzustands durch Unterraumprojektion und Matrixdiagonalisierung.

1. Problem auf Quantenschaltungen und Operatoren abbilden

Molekularer Hamilton-Operator

Ein molekularer Hamiltonian hat die allgemeine Form:

$\hat{a}^\dagger_{p\sigma}$ / $\hat{a}_{p\sigma}$ sind die fermionischen Erzeugungs-/Annihilationsoperatoren, die mit dem $p$ -ten Basiselement und dem Spin $\sigma$ verbunden sind. $h_{pr}$ und $(pr|qs)$ sind die elektronischen Ein- und Zweikörperintegrale. Mit pySCF, definieren wir das Molekül und berechnen die Ein- und Zweikörperintegrale des Hamiltonianers für den Basissatz 6-31g.

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

Output:

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

In dieser Lektion werden wir die Jordan-Wigner-Transformation (JW) verwenden, um eine fermionische Wellenfunktion auf eine Qubit-Wellenfunktion abzubilden, so dass sie mit einem Quantenschaltkreis vorbereitet werden kann. Die JW-Transformation bildet den Fock-Raum der Fermionen in M räumlichen Orbitalen auf den Hilbert-Raum von 2M Qubits ab, d. h. ein räumliches Orbital wird in zwei Spin-Orbitale aufgeteilt, von denen eines mit einem Elektron mit Spin nach oben ( $\alpha$ ) und das andere mit Spin nach unten ( $\beta$ ) verbunden ist. Ein Spin-Orbital kann besetzt oder unbesetzt sein. Wenn wir von der Anzahl der Orbitale sprechen, meinen wir normalerweise die Anzahl der räumlichen Orbitale. Die Anzahl der Spin-Orbitale wird verdoppelt. In Quantenschaltungen wird jedes Spin-Orbital mit einem Qubit dargestellt. Eine Gruppe von Qubits steht also für Spin-Up- oder $\alpha$ -Orbitale und eine andere Gruppe für Spin-Down- oder $\beta$ -Orbitale. Zum Beispiel hat das $N_2$ Molekül für den 6-31g Basissatz $16$ räumliche Orbitale (d.h. $16$ $\alpha$ + $16$ $\beta$ = $32$ Spin-Orbitale). Wir brauchen also einen $32$ -Quantenschaltkreis (eventuell brauchen wir zusätzliche Ancilla-Qubits, wie später besprochen). Die Qubits werden auf rechnerischer Basis gemessen, um Bitstrings zu erzeugen, die elektronische Konfigurationen oder (Slater-)Determinanten darstellen. In dieser Lektion werden wir die Begriffe Bitstrings, Konfigurationen und Determinanten synonym verwenden. Die Bitstrings geben Aufschluss über die Elektronenbelegung in den Spin-Orbitalen: Ein $1$ an einer Bitposition bedeutet, dass das entsprechende Spin-Orbital besetzt ist, während ein $0$ bedeutet, dass das Spin-Orbital leer ist. Da Probleme der elektronischen Struktur teilchenerhaltend sind, muss nur eine bestimmte Anzahl von Spin-Orbitalen besetzt sein. Das Molekül $N_2$ hat $5$ Elektronen mit aufsteigendem Spin ( $\alpha$ ) und $5$ Elektronen mit absteigendem Spin ( $\beta$ ). Daher muss jede Bitfolge, die die Orbitale $\alpha$ und $\beta$ repräsentiert, jeweils fünf $1\text{s}$ für das Molekül $N_2$ enthalten.

1.1 Quantenschaltung zur Erzeugung von Proben: Der LUCJ-Ansatz

In dieser Lektion werden wir den LUCJ-Ansatz (local unitary coupled cluster Jastrow) \ [1] für die Vorbereitung von Quantenzuständen und die anschließende Probenahme verwenden. Zunächst werden wir die verschiedenen Bausteine des vollständigen UCJ-Ansatzes und die Näherungen in der lokalen Version davon erklären. Als nächstes werden wir mit Hilfe des ffsim-Pakets den LUCJ-Ansatz konstruieren und ihn mit dem Qiskit-Transpiler für die Hardware-Ausführung optimieren.

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

wobei $\vert \Phi_{0} \rangle$ ein Referenzzustand ist, der üblicherweise als Hartree-Fock-Zustand (HF) angenommen wird. Da der Hartree-Fock-Zustand so definiert ist, dass die Orbitale mit der niedrigsten Nummer besetzt sind, müssen bei der Vorbereitung des HF-Zustands X-Gates angewendet werden, um die den besetzten Orbitalen entsprechenden Qubits auf Eins zu setzen. Der HF-Zustandsvorbereitungsblock für 4 räumliche Orbitale und 2 Up- und 2 Down-Spins kann zum Beispiel 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-Evolution ( $e^{iJ^{(\mu)}}$ ), die von orbitalen Rotationen ( $e^{K^{(\mu)}}$ und $e^{-K^{(\mu)}}$ ) unterbrochen wird.

Orbitalrotationsblöcke arbeiten mit einer einzigen Spinart ( $\alpha$ (up-spin)/ $\beta$ (down-spin)). Für jede Elektronenart besteht die Orbitalrotation aus einer Schicht von Ein-Qubit-Gattern ( $R_{z}$ ), gefolgt von einer Folge von 2-Qubit-Gattern ( $XX + YY$ ).

Die 2-Qubit-Gatter wirken auf benachbarte Spin-Orbitale (nächstgelegene Qubits) und sind daher auf IBM® QPUs ohne SWAP-Gatter implementierbar.

Der $e^{iJ^{(\mu)}}$, auch bekannt als diagonaler Coulomb-Operator, besteht aus drei Blöcken. Zwei von ihnen arbeiten auf demselben Spin-Sektor ( $e^{iJ_{\alpha \alpha}^{(\mu)}}$ und $e^{iJ_{\beta \beta}^{(\mu)}}$ ), und einer arbeitet zwischen zwei Spin-Sektoren ( $e^{iJ_{\alpha \beta}^{(\mu)}}$ ).

Alle Blöcke in $e^{iJ^{(\mu)}}$ bestehen aus Zahlen-Zahlen-Gattern $U_{nn}(\phi)$ [1]. Ein $U_{nn}(\phi)$ Gatter kann weiter unterteilt werden in ein $R_{ZZ}(\frac{\phi}{2})$ Gatter, gefolgt von zwei Einzel-Qubit $Rz(-\frac{\phi}{2})$ Gattern, die auf zwei separate Qubits wirken.

Gleichspin-Komponenten ( $J_{\alpha \alpha}$ und $J_{\beta \beta}$ ) haben $U_{nn}$ Gatter zwischen allen möglichen Qubit-Paaren. Da supraleitende QPUs jedoch nur über eine begrenzte Konnektivität verfügen, müssen Qubits ausgetauscht werden, um Gatter zwischen nicht benachbarten Qubits zu realisieren.

Betrachten Sie zum Beispiel den folgenden $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (oder $e^{iJ_{\beta \beta}^{(\mu)}}$ ) Block für $N = 4$ räumliche Orbitale. Bei einer linearen Qubit-Konnektivität sind die letzten drei Gatter nicht direkt implementierbar, da sie zwischen nicht benachbarten Qubits funktionieren (zum Beispiel sind Q0 und Q2 nicht direkt verbunden). Daher benötigen wir SWAP-Gatter, um sie nebeneinander zu schalten (die folgende Abbildung zeigt ein Beispiel mit $3$ SWAP-Gattern).

Als Nächstes implementiert $J_{\alpha \beta}$ Gates zwischen gleichen indizierten Orbitalen aus verschiedenen Spin-Sektoren (z. B. zwischen $0\alpha$ und $0\beta$ ). Wenn die Qubits auf einer QPU nicht physisch benachbart sind, erfordern diese Gatter ebenfalls SWAPs.

Aus der obigen Diskussion geht hervor, dass der UCJ-Ansatz einige Hürden für die HW-Ausführung aufweist, da er aufgrund der Wechselwirkungen zwischen nicht benachbarten Qubits SWAP-Gatter benötigt. Die lokale Variante des UCJ-Ansatzes, LUCJ, löst diese Herausforderung, indem sie einige $U_{nn}$ aus dem diagonalen Coulomb-Operator entfernt.

In denselben Elektronenspezies-Blöcken ( $J_{\alpha \alpha}$ und $J_{\beta \beta}$ ) behalten wir nur die $U_{nn}$ Gatter bei, die mit der Konnektivität der nächsten Nachbarn kompatibel sind, und entfernen die Gatter zwischen nicht benachbarten Qubits in der LUCJ-Version. Die folgende Abbildung zeigt den LUCJ-Block nach Entfernung der nicht benachbarten Gatter.

Als nächstes kann die LUCJ-Version des $J_{\alpha \beta}$ -Blocks, der zwischen verschiedenen Elektronenspezies arbeitet, je nach Gerätetopologie unterschiedliche Formen annehmen.

Auch hier beseitigt die LUCJ-Version die nicht kompatiblen 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 Wechselwirkungen $\alpha$ - $\beta$ werden zwischen jedem $4$ -ten indizierten (z.B. dem 0., 4. und 8.) Spin-Orbital aufrechterhalten und sind Ancilla-vermittelt, d.h. wir brauchen Ancilla-Qubits zwischen den linearen Ketten, die $\alpha$ und $\beta$ Orbitale darstellen. Diese Regelung erfordert eine begrenzte Anzahl von SWAPs.
• Hexagonal : Jedes andere Orbital, wie das 0., 2. und 4. indizierte Orbital, wird zum 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 Tor hat.

Die Entfernung von Gattern aus dem UCJ-Ansatz zur Konstruktion der LUCJ-Version macht ihn zwar HW-kompatibler, aber der Ansatz verliert an Ausdruckskraft. Daher können bei Verwendung des LUCJ-Ansatzes mehr Wiederholungen ( $L$ ) des modifizierten UCJ-Operators erforderlich sein.

1.2 LUCJ-Ansatz-Initialisierung

Der LUCJ ist ein parametrisierter Ansatz, und wir müssen die Parameter vor der Hardware-Ausführung initialisieren. Eine Möglichkeit, den Ansatz zu initialisieren, ist die Verwendung von t1 und t2 Amplituden aus der klassischen gekoppelten Cluster-Singles and Doubles (CCSD)-Methode, wobei t1 Amplituden die Koeffizienten von Einzelanregungsoperatoren und t2 Amplituden für Doppelanregungsoperatoren sind.

Beachten Sie, dass die Initialisierung des LUCJ-Ansatzes mit den Amplituden von t1 und t2 zwar zu brauchbaren Ergebnissen führt, die Ansatzparameter aber möglicherweise weiter optimiert werden müssen.

# 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

Output:

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 Konstruktion des LUCJ-Ansatzes unter Verwendung von ffsim

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

Da die Hardware von IBM eine Heavy-Hex-Topologie hat, werden wir das in [1] verwendete und oben erläuterte Zick-Zack-Muster für Qubit-Interaktionen übernehmen. In diesem Muster sind Orbitale (Qubits) mit gleichem Spin durch eine Linientopologie verbunden (rote und blaue Kreise). Aufgrund der Heavy-Hex-Topologie haben Orbitale für verschiedene Spins Verbindungen zwischen jedem 4. Orbital, d. h. 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-Satz mit wiederholten Schichten kann durch Zusammenlegung einiger benachbarter Blöcke optimiert werden. Betrachten Sie einen Fall für n_reps=2. Die beiden Orbitalrotationsblöcke in der Mitte können zu einem einzigen Orbitalrotationsblock zusammengefasst werden. Das Paket ffsim verfügt über einen Pass-Manager namens ffsim.qiskit.PRE_INIT , der die Schaltung durch Zusammenführen solcher benachbarter Blöcke optimiert.

2. Für die Zielhardware optimieren

Zunächst holen wir uns ein Backend unserer Wahl. Wir werden unsere Schaltung für das Backend optimieren und dann die optimierte Schaltung auf demselben Backend ausführen, um Proben 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 ihn hardwarekompatibel zu machen.

• Wählen Sie physikalische Qubits (initial_layout) aus der Zielhardware aus, die dem oben beschriebenen Zick-Zack-Muster (zwei lineare Ketten mit einem Ancilla-Qubit dazwischen) entsprechen. Die Anordnung der Qubits nach diesem Muster führt zu einer effizienten, hardwarekompatiblen Schaltung mit weniger Gattern.
• Erzeugen Sie einen gestaffelten Passmanager mit der generate_preset_pass_manager funktion von Qiskit mit Ihrer Wahl von backend und initial_layout.
• Setzen Sie die Stufe pre_init Ihres Staged Pass Managers auf ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT enthält Qiskit-Transpiler-Passagen, die Gatter in orbitale Rotationen zerlegen und dann die orbitalen Rotationen zusammenführen, was zu weniger Gattern in der endgültigen Schaltung führt.
• Führen Sie den Passmanager auf Ihrer Schaltung aus.

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

Output:

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

Nachdem wir die Schaltung für die Ausführung auf der Hardware optimiert haben, können wir sie nun auf der Zielhardware ausführen und Messwerte zur Schätzung der Grundzustandsenergie erfassen. Da wir nur einen Schaltkreis haben, werden wir den qiskit-ibm-runtime Job-Ausführungsmodus verwenden und unseren Schaltkreis ausführen.

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. Nachbearbeitungsergebnisse

Der Post-Processing-Teil des SQD-Workflows lässt sich anhand des folgenden Diagramms zusammenfassen.

Die Abtastung des LUCJ-Ansatzes in der Berechnungsbasis erzeugt einen Pool von verrauschten Konfigurationen $\tilde{\mathcal{\chi}}$, die in der Nachbearbeitungsroutine verwendet werden. Es handelt sich um eine Methode, die Konfigurationswiederherstellung genannt wird (Einzelheiten werden später erörtert), um Konfigurationen mit falschen Elektronenzahlen probabilistisch zu korrigieren. Nur Konfigurationen mit korrekten Elektronennummern $\tilde{\mathcal{\chi}}_{R}$ werden dann in eine Teilstichprobe aufgenommen und auf der Grundlage der Häufigkeit des Auftretens jeder einzelnen Konfiguration auf mehrere Lose verteilt. Jede Gruppe von Proben definiert einen Unterraum ( $\mathcal{S^{(k)}}$ ). Als nächstes wird der molekulare Hamiltonian $H$ auf Unterräume projiziert:

Jeder projizierte Hamiltonian $H_{\mathcal{S}^{(k)}}$ wird dann in einen Eigensolver eingespeist, wo er diagonalisiert wird, um Eigenwerte und Eigenvektoren zu berechnen und einen Eigenzustand zu rekonstruieren. In dieser Lektion projizieren und diagonalisieren wir den Hamiltonian mit Hilfe des Pakets qiskit-addon-sqd , das die Davidson-Methode aus PySCF zur Diagonalisierung verwendet.

Wir sammeln dann den niedrigsten Eigenwert (Energie) aus den Stapeln und berechnen auch die durchschnittliche Orbitalbelegung $\text{n}$. Die Informationen über die durchschnittliche Belegung werden im Schritt der Konfigurationswiederherstellung verwendet, um die Rauschkonfigurationen probabilistisch zu korrigieren.

Als nächstes wird die Schleife zur Wiederherstellung der selbstkonsistenten Konfiguration im Detail erläutert und es werden konkrete Code-Beispiele zur Implementierung der oben genannten Schritte zur Abschätzung der Grundzustandsenergie des $N_2$ -Hamiltonianers gezeigt.

4.1 Konfigurationswiederherstellung: Übersicht

Jedes Bit in einer Bitkette (Slater-Determinante) steht für ein Spin-Orbital. Die rechte Hälfte eines Bitstrings steht für Spin-up-Orbitale, die linke Hälfte für Spin-down-Orbitale. Ein 1 bedeutet, dass das Orbital mit einem Elektron besetzt ist, und ein 0 bedeutet, dass das Orbital leer ist. Wir kennen die richtige Anzahl von Teilchen (sowohl Up-Spin-Elektron als auch Down-Spin-Elektron) a priori. Nehmen wir an, wir haben eine Determinante $x$ mit $N_x$ Elektronen (d. h. es gibt $N_x$ Anzahl von $1$ s in der Bitkette) darin. Die korrekte Anzahl von Partikeln ist $N$. Wenn $N_x \neq N$, dann wissen wir, dass der Bitstring durch Rauschen verfälscht ist. Die selbstkonsistente Konfigurationsroutine versucht, den Bitstring zu korrigieren, indem sie die $|N_x - N|$ Bits mit Hilfe von Informationen über die durchschnittliche Orbitalbesetzung probabilistisch umkehrt. Die durchschnittliche Orbitalbesetzung ( $n$ ) gibt an, wie wahrscheinlich es ist, dass ein Orbital von einem Elektron besetzt ist. Wenn $N_x < N$, haben wir weniger Elektronen und müssen einige $0$ s in $1$ s umwandeln und umgekehrt.

Die Wahrscheinlichkeit des Flippens kann unter $|x[i] - avg\_occupancy[i]|$ für i-th spin orbital. In [2] verwendeten die Autoren eine gewichtete Umdrehungswahrscheinlichkeit unter Verwendung der modifizierten Funktion ReLU.

Hier definiert $h$ die Position der "Ecke" der Funktion ReLU, und der Parameter $\delta$ definiert den Wert der Funktion ReLU an der Ecke. Für $\delta = 0$ wird $w$ zur echten Funktion ReLU und für $\delta >0$ zur modifizierten Funktion ReLU. In dem Papier 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 durchschnittliche Orbitalbesetzung ( $n$ ) ist nicht von vornherein bekannt. Die erste Iteration der Grundzustandsschätzung beginnt mit Konfigurationen, die nur korrekte Teilchenzahlen in beiden Spinarten aufweisen. Nach der ersten Iteration haben wir eine Schätzung des Grundzustands, und mit Hilfe der Schätzung können wir die erste Schätzung von $n$ konstruieren. Diese Schätzung von $n$ wird verwendet, um Konfigurationen wiederherzustellen, die nächste Iteration der Grundzustandsschätzung durchzuführen und die Schätzung von $n$ selbstkonsistent zu verfeinern. Der Prozess wird wiederholt, bis ein Stoppkriterium erfüllt ist.

Betrachten Sie das folgende Beispiel für $N = 2$ und $x = |1000\rangle$ ( $N_x = 1$ ). Wir müssen eines der 0s auf 1 umdrehen, um es für die Partikelzahlen zu korrigieren, und die Auswahlmöglichkeiten sind 1100, 1010 und 1001. Basierend auf der Wahrscheinlichkeit des Umkippens wird eine der Möglichkeiten als wiederhergestellte Konfiguration ausgewählt (oder der Bitstring mit der richtigen Anzahl von Partikeln).

Nehmen wir an, dass wir in der ersten Iteration zwei Stapel laufen lassen, und die geschätzten Grundzustände daraus sind:

Mit Hilfe der Basiszustände und ihrer Amplituden können wir die Wahrscheinlichkeit der Elektronenbelegung (kurz: Belegung ) pro Spin-Orbital (Qubit) berechnen (man beachte, dass Wahrscheinlichkeit = Amplitude = $^2$ ). Im Folgenden werden die Qubit-weisen Besetzungen für jeden Bitstring im geschätzten Grundzustand tabellarisch dargestellt und die Gesamtorbitalbesetzung für eine Charge berechnet. Beachten Sie, dass gemäß der Qiskit-Ordnungskonvention das rechteste Bit für qubit-0 ( Q0 ) und das linkste Bit für Q3 steht.

Belegung ( Batch0 ):

Belegung ( Batch1 )

Belegung (Durchschnitt aller Chargen)

Anhand der oben berechneten durchschnittlichen Orbitalbesetzung können wir die Umkehrwahrscheinlichkeiten für verschiedene Orbitale in der Konfiguration $x = \vert 1000 \rangle$ ermitteln. Da das durch Q3 dargestellte Orbital bereits besetzt ist und nicht umgedreht werden muss, setzen wir sein p(flip) auf $0$. Für die übrigen Orbitale, die unbesetzt sind, beträgt die Umkippwahrscheinlichkeit jeweils $\vert x[i] - \text{n}[i] \vert$. Zusammen mit p(flip) berechnen wir auch die Wahrscheinlichkeitsgewichtung, die mit dem Umdrehen verbunden ist, indem wir die oben beschriebene modifizierte Funktion ReLU verwenden.

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

Schließlich können wir mit Hilfe der oben genannten gewichteten Wahrscheinlichkeiten eines der unbesetzten Orbitale Q2, Q1 und Q0 umdrehen. Basierend auf den obigen Werten wird Q0 höchstwahrscheinlich umgedreht, und eine mögliche wiederhergestellte Konfiguration kann $\vert \text{1001} \rangle$ sein.

Der gesamte Prozess der Wiederherstellung einer selbstkonsistenten Konfiguration 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 richtiger ( $\widetilde{\chi}_{correct}$ ) als auch mit falscher ( $\widetilde{\chi}_{incorrect}$ ) Anzahl von Teilchen in jedem Spin-Sektor enthält.

1. Die Konfigurationen aus ( $\widetilde{\chi}_{correct}$ ) werden nach dem Zufallsprinzip ausgewählt, um Stapel $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ von Vektoren für die Unterraumprojektion zu erstellen. Die Anzahl der Chargen und Proben in jeder Charge sind benutzerdefinierte Parameter. Je größer die Anzahl der Proben in jeder Charge ist, desto größer ist die Dimension des Unterraums und desto rechenintensiver wird die Diagonalisierung. Andererseits können bei einer zu geringen Anzahl von Stichproben die Stützvektoren des Grundzustands fehlen, was zu einer falschen Schätzung führt.
2. Führen Sie den Eigenzustandslöser (d. h. Projektion auf den Unterraum und Diagonalisierung) auf den Stapeln aus und erhalten Sie ungefähre Eigenzustände. $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. Konstruieren Sie aus den angenäherten Eigenzuständen die erste Schätzung für $n$.

Nachfolgende Iterationen:

1. Korrigieren Sie mit $n$ die Konfigurationen mit falscher Teilchenzahl in $\widetilde{\chi}_{incorrect}$. Nehmen wir an, wir nennen sie $\widetilde{\chi}_{correct\_new}$. Dann bildet $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ den neuen Satz von Konfigurationen mit korrekten Teilchenzahlen.
2. $\widetilde{\chi}_{R}$ wird eine Stichprobe gezogen, um Lose zu erstellen $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. Der Eigenzustandslöser läuft mit neuen Stapeln und erzeugt neue Schätzungen der Grundzustände $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. Konstruieren Sie aus den angenäherten Eigenzuständen eine verfeinerte Schätzung für $n$.
5. Wenn das Abbruchkriterium nicht erfüllt ist, gehen Sie zurück zu Schritt 2.1.

4.2 Grundzustandsschätzung

Zunächst werden die Zählungen in eine Bitstring-Matrix und ein Wahrscheinlichkeitsfeld für die Nachbearbeitung umgewandelt.

Jede Zeile in der Matrix steht für einen eindeutigen Bitstring. Da die Qubits in Qiskit von der rechten Seite eines Bitstrings indiziert werden, steht die Spalte 0 für das Qubit N-1 und die Spalte N-1 für das Qubit 0, wobei N die Anzahl der Qubits ist.

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

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 ein paar benutzergesteuerte Optionen, die für diese Technik wichtig sind:

• iterations: Anzahl der Iterationen zur Wiederherstellung der selbstkonsistenten Konfiguration
• n_batches: Anzahl der Stapel von Konfigurationen, die von den verschiedenen Aufrufen des Eigenzustandslösers verwendet werden
• samples_per_batch: Anzahl der eindeutigen Konfigurationen, die in jeder Charge enthalten sein sollen
• max_davidson_cycles: Maximale Anzahl von Davidson-Zyklen, die von jedem Eigensolver 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 postselect 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_cycle=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)

Output:

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 Erörterung der Ergebnisse

Die erste Grafik zeigt, dass wir nach einigen Iterationen die Grundzustandsenergie mit einer Genauigkeit von ~24 mH schätzen (die chemische Genauigkeit wird normalerweise mit 1 kcal/mol angenommen $\approx$ 1.6 mH ). Das zweite Diagramm zeigt die durchschnittliche Belegung der einzelnen räumlichen Orbitale nach der letzten Iteration. Wir können sehen, dass sowohl die Spin-up- als auch die Spin-down-Elektronen die ersten fünf Orbitale mit hoher Wahrscheinlichkeit in unseren Lösungen besetzen.

Obwohl die geschätzte Grundzustandsenergie gut ist, liegt sie nicht innerhalb der chemischen Genauigkeitsgrenze ( $\pm \approx 1.6$ mH ). Diese Lücke kann auf die kleine Unterraumdimension zurückgeführt werden, die wir oben für Projektion und Diagonalisierung verwendet haben. Da wir samples_per_batch=500 verwendet haben, wird der Unterraum von maximal $500$ Vektoren aufgespannt, d. h. von fehlenden Vektoren der Grundzustandsunterstützung. 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()

Output:

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

Übung für den Leser

Erhöhen Sie schrittweise den Parameter samples_per_batch (z. B. von $1000$ auf $10000$ in einem Schritt von $1000$; dies ist mit dem Speicher Ihres Computers möglich) und vergleichen Sie die geschätzten Grundzustandsenergien.

Referenzen

[1] M. Motta et al, "Brückenschlag zwischen physikalischer Intuition und Hardware-Effizienz für korrelierte elektronische Zustände: der lokale unitäre Cluster-Jastrow-Ansatz für die elektronische Struktur" (2023). Chem. Naturwissenschaft, 2023, 14, 11213.

[2] J. Robledo-Moreno et al, "Chemie jenseits exakter Lösungen auf einem quantenzentrierten Supercomputer" (2024). arXiv:quant-ph/2405.05068.