{
"cells": [
{
"cell_type": "markdown",
"id": "090b6884",
"metadata": {},
"source": [
"---\n",
"title: \"SQD-Implementierung\"\n",
"description: \"Die samplebasierte Quantendiagonalisierung (SQD) wird im Zusammenhang mit der Lösung des Grundzustands eines Stickstoffmoleküls implementiert.\"\n",
"---\n",
"\n",
"{/* cspell:ignore fontdict fontsize milli hcore ccsd Motta */}\n",
"\n",
"\n",
"\n",
"# SQD zur Energieabschätzung eines chemischen Hamilton-Operators\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ee7de1bc",
"metadata": {},
"source": [
"In dieser Lektion werden wir SQD anwenden, um die Grundzustandsenergie eines Moleküls zu schätzen.\n",
"\n",
"Insbesondere werden wir die folgenden Themen anhand des $4$ -step Qiskit-Musteransatzes diskutieren:\n",
"\n",
"1. Schritt 1: Abbildung des Problems auf Quantenschaltungen und -operatoren\n",
" * Stellen Sie den molekularen Hamiltonian für $N_2$ auf.\n",
" * Erläutern Sie den von der Chemie inspirierten und hardwarefreundlichen lokalen unitären Cluster Jastrow (LUCJ) [\\[1\\]](#references)\n",
"2. Schritt 2: Optimieren für die Zielhardware\n",
" * Optimierung der Gatteranzahl und des Layouts des Ansatzes für die Hardware-Ausführung\n",
"3. Schritt 3: Ausführen auf der Zielhardware\n",
" * Führen Sie die optimierte Schaltung auf einer echten QPU aus, um Stichproben des Unterraums zu erzeugen.\n",
"4. Schritt 4: Nachbearbeitung der Ergebnisse\n",
" * Einführung der Schleife zur Wiederherstellung der selbstkonsistenten Konfiguration [\\[2\\]](#references)\n",
" * 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.\n",
" * Wahrscheinliche Erstellung von Stapeln von Unterstichproben aus wiederhergestellten Bitstrings.\n",
" * Projizieren und diagonalisieren Sie den molekularen Hamiltonian über jeden abgetasteten Unterraum.\n",
" * Speichern Sie die minimale Grundzustandsenergie, die in allen Stapeln gefunden wurde, und aktualisieren Sie die durchschnittliche Orbitalbesetzung.\n",
"\n",
"Während des Unterrichts werden wir verschiedene Softwarepakete verwenden.\n",
"\n",
"* `PySCF` um das Molekül zu definieren und den Hamiltonian aufzustellen.\n",
"* `ffsim` paket, um den LUCJ-Ansatz zu konstruieren.\n",
"* `Qiskit` zum Transpilieren des Ansatzes für die Hardware-Ausführung.\n",
"* `Qiskit IBM Runtime` um die Schaltung auf einer QPU auszuführen und Proben zu sammeln.\n",
"* `Qiskit addon SQD` wiederherstellung der Konfiguration und Schätzung der Energie des Grundzustands durch Unterraumprojektion und Matrixdiagonalisierung.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "719a9c0e-8c00-4fab-ba23-f2c8a7ebb573",
"metadata": {},
"source": [
"\n",
"\n",
"## 1. Problem auf Quantenschaltungen und Operatoren abbilden\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d02f97af",
"metadata": {},
"source": [
"\n",
"\n",
"### Molekularer Hamilton-Operator\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "5dac27ce-56b1-4f7f-ab83-ac153c004e82",
"metadata": {},
"source": [
"Ein molekularer Hamiltonian hat die allgemeine Form:\n",
"\n",
"$$\n",
"\\hat{H} = \\sum_{ \\substack{pr\\\\\\sigma} } h_{pr} \\, \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{r\\sigma}\n",
"+\n",
"\\sum_{ \\substack{prqs\\\\\\sigma\\tau} }\n",
"\\frac{(pr|qs)}{2} \\,\n",
"\\hat{a}^\\dagger_{p\\sigma}\n",
"\\hat{a}^\\dagger_{q\\tau}\n",
"\\hat{a}_{s\\tau}\n",
"\\hat{a}_{r\\sigma}\n",
"$$\n",
"\n",
"$\\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`.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cc987c08-7261-4c4c-a06b-609d7003efe9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"converged SCF energy = -108.835236570774\n",
"CASCI E = -109.046671778080 E(CI) = -32.8155692383188 S^2 = 0.0000000\n"
]
}
],
"source": [
"import warnings\n",
"import pyscf\n",
"import pyscf.cc\n",
"import pyscf.mcscf\n",
"\n",
"warnings.filterwarnings(\"ignore\")\n",
"\n",
"# Specify molecule properties\n",
"open_shell = False\n",
"spin_sq = 0\n",
"\n",
"# Build N2 molecule\n",
"mol = pyscf.gto.Mole()\n",
"mol.build(\n",
" atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]], # Two N atoms 1 angstrom apart\n",
" basis=\"6-31g\",\n",
" symmetry=\"Dooh\",\n",
")\n",
"\n",
"# Define active space\n",
"n_frozen = 2\n",
"active_space = range(n_frozen, mol.nao_nr())\n",
"\n",
"# Get molecular integrals\n",
"scf = pyscf.scf.RHF(mol).run()\n",
"num_orbitals = len(active_space)\n",
"n_electrons = int(sum(scf.mo_occ[active_space]))\n",
"num_elec_a = (n_electrons + mol.spin) // 2\n",
"num_elec_b = (n_electrons - mol.spin) // 2\n",
"cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))\n",
"mo = cas.sort_mo(active_space, base=0)\n",
"hcore, nuclear_repulsion_energy = cas.get_h1cas(mo) # hcore: one-body integrals\n",
"eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals) # eri: two-body integrals\n",
"\n",
"# Compute exact energy for comparison\n",
"exact_energy = cas.run().e_tot"
]
},
{
"cell_type": "markdown",
"id": "3aa4e0f0",
"metadata": {},
"source": [
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "f6a89c89-85e1-4c0d-abeb-8ab8b154a7ba",
"metadata": {},
"source": [
"\n",
"\n",
"### 1.1 Quantenschaltung zur Erzeugung von Proben: Der LUCJ-Ansatz\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "a464c865-1528-45c2-8ede-325583b15976",
"metadata": {},
"source": [
"In dieser Lektion werden wir den LUCJ-Ansatz (local unitary coupled cluster Jastrow) \\ \\[[1\\]](#references) 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.\n",
"\n",
"Der UCJ-Ansatz hat die folgende Form (für ein Produkt aus $L$ Schichten oder Wiederholungen des UCJ-Operators)\n",
"\n",
"$$\n",
"|\\psi\\rangle = \\prod_{\\mu=1}^{L}{(e^{K^{\\mu}} \\times {e^{iJ^{\\mu}}} \\times {e^{-K^{\\mu}}})} |\\Phi_{0}\\rangle\n",
"$$\n",
"\n",
"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:\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "0f0bb614",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "4a32d88b",
"metadata": {},
"source": [
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "963e8386-39b6-40b7-9740-fffcf1573fe6",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "271b1825-bc78-4957-bbf8-21a714e45ed3",
"metadata": {},
"source": [
"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$ ).\n",
"\n",
"Die 2-Qubit-Gatter wirken auf benachbarte Spin-Orbitale (nächstgelegene Qubits) und sind daher auf IBM® QPUs ohne SWAP-Gatter implementierbar.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "2e8ac1d2-8f04-4591-921b-8ba0174e4ad0",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ebc9b48c-e77d-4af5-b8cf-460daeeadcdb",
"metadata": {},
"source": [
"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)}}$ ).\n",
"\n",
"Alle Blöcke in $e^{iJ^{(\\mu)}}$ bestehen aus Zahlen-Zahlen-Gattern $U_{nn}(\\phi)$ [\\[1\\]](#references). 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.\n",
"\n",
"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.\n",
"\n",
"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).\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d3e24a20-1c86-4aea-8300-13bd2ead4e8f",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ef3f1d36-e8ab-46f0-bddb-fa8fb884e531",
"metadata": {},
"source": [
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9afe0036-318c-43b9-9b2c-39a808bda82c",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "6219aa8e-8a63-4bfc-b52b-64507292471d",
"metadata": {},
"source": [
"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.\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "24f54a55-4a09-4508-997a-96306835c7e6",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "165d0a19-0366-4af8-8ea0-20df351f49f5",
"metadata": {},
"source": [
"Als nächstes kann die LUCJ-Version des $J_{\\alpha \\beta}$ -Blocks, der zwischen verschiedenen Elektronenspezies arbeitet, je nach Gerätetopologie unterschiedliche Formen annehmen.\n",
"\n",
"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.\n",
"\n",
"* **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.\n",
"* **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.\n",
"* **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.\n",
"* **Linear** : Nur ein $\\alpha$ und ein $\\beta$ Orbital sind verbunden, was bedeutet, dass der $J_{\\alpha \\beta}$ Block nur ein Tor hat.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8ec1e433-bc00-42bc-bdbb-288c56c32f9d",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "3e070615-2b4a-4866-8c6c-d2d038447f45",
"metadata": {},
"source": [
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "04367dac",
"metadata": {},
"source": [
"\n",
"\n",
"### 1.2 LUCJ-Ansatz-Initialisierung\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "f5eae55f",
"metadata": {},
"source": [
"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.\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "1835e74e-3354-425f-8596-c574f03e7a6e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"E(CCSD) = -109.0398256929733 E_corr = -0.20458912219883\n"
]
}
],
"source": [
"# Get CCSD t2 amplitudes for initializing the ansatz\n",
"ccsd = pyscf.cc.CCSD(\n",
" scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]\n",
")\n",
"ccsd.run()\n",
"\n",
"t1 = ccsd.t1\n",
"t2 = ccsd.t2"
]
},
{
"cell_type": "markdown",
"id": "9f0842fa",
"metadata": {},
"source": [
"\n",
"\n",
"### 1.3 Konstruktion des LUCJ-Ansatzes unter Verwendung von `ffsim`\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "e5f7f7e8-30e3-40e8-9b85-bf057b35b766",
"metadata": {},
"source": [
"Wir werden das [ffsim-Paket](https://github.com/qiskit-community/ffsim/tree/main) 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](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced).\n",
"\n",
"Da die Hardware von IBM eine Heavy-Hex-Topologie hat, werden wir das in [\\[1\\]](#references) 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).\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "6799421a-4425-404a-9994-47215957054d",
"metadata": {},
"outputs": [],
"source": [
"import ffsim\n",
"from qiskit import QuantumCircuit, QuantumRegister\n",
"\n",
"n_reps = 2\n",
"alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]\n",
"alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]\n",
"\n",
"ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
" t2=t2,\n",
" t1=t1,\n",
" n_reps=n_reps,\n",
" interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),\n",
")\n",
"\n",
"nelec = (num_elec_a, num_elec_b)\n",
"\n",
"# create an empty quantum circuit\n",
"qubits = QuantumRegister(2 * num_orbitals, name=\"q\")\n",
"circuit = QuantumCircuit(qubits)\n",
"\n",
"# prepare Hartree-Fock state as the reference state and append it to the quantum circuit\n",
"circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)\n",
"\n",
"# apply the UCJ operator to the reference state\n",
"circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n",
"circuit.measure_all()\n",
"# circuit.decompose().draw(\"mpl\", scale=0.5, fold=-1)"
]
},
{
"cell_type": "markdown",
"id": "3b2de316",
"metadata": {},
"source": [
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "7cb99cd9",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "e9ad460d",
"metadata": {},
"source": [
"\n",
"\n",
"## 2. Für die Zielhardware optimieren\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "cecf2994",
"metadata": {},
"source": [
"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.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4c8dbbef-378e-4ae3-9290-bde58b72024d",
"metadata": {},
"outputs": [],
"source": [
"from qiskit_ibm_runtime import QiskitRuntimeService\n",
"\n",
"service = QiskitRuntimeService()\n",
"# Use the least-busy backend or specify a quantum computer using the syntax commented out below.\n",
"backend = service.least_busy(operational=True, simulator=False)\n",
"# backend = service.backend(\"ibm_brisbane\")"
]
},
{
"cell_type": "markdown",
"id": "1b8ba388-6904-4725-ad97-2a31c0adaa07",
"metadata": {},
"source": [
"Als Nächstes empfehlen wir die folgenden Schritte, um den Ansatz zu optimieren und ihn hardwarekompatibel zu machen.\n",
"\n",
"* 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.\n",
"* Erzeugen Sie einen gestaffelten Passmanager mit der [`generate_preset_pass_manager`](/docs/api/qiskit/qiskit.transpiler.generate_preset_pass_manager) funktion von Qiskit mit Ihrer Wahl von `backend` und `initial_layout`.\n",
"* 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.\n",
"* Führen Sie den Passmanager auf Ihrer Schaltung aus.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "ecd2675c-a302-43fb-80f4-d658d56360d5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Gate counts (w/o pre-init passes): OrderedDict({'rz': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})\n",
"Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})\n"
]
}
],
"source": [
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
"\n",
"spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]\n",
"spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]\n",
"\n",
"initial_layout = spin_a_layout + spin_b_layout\n",
"\n",
"pass_manager = generate_preset_pass_manager(\n",
" optimization_level=3, backend=backend, initial_layout=initial_layout\n",
")\n",
"\n",
"# without PRE_INIT passes\n",
"isa_circuit = pass_manager.run(circuit)\n",
"print(f\"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}\")\n",
"\n",
"# with PRE_INIT passes\n",
"# We will use the circuit generated by this pass manager for hardware execution\n",
"pass_manager.pre_init = ffsim.qiskit.PRE_INIT\n",
"isa_circuit = pass_manager.run(circuit)\n",
"print(f\"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}\")"
]
},
{
"cell_type": "markdown",
"id": "72a780f9",
"metadata": {},
"source": [
"\n",
"\n",
"## 3. Auf der Zielhardware ausführen\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "4cd9164c-697a-4aa6-b71f-86720e3d5b66",
"metadata": {},
"source": [
"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](/docs/guides/execution-modes) verwenden und unseren Schaltkreis ausführen.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "9f542448-5767-4e12-85c8-dd8584545dbb",
"metadata": {},
"outputs": [],
"source": [
"from qiskit_ibm_runtime import SamplerV2 as Sampler\n",
"\n",
"sampler = Sampler(mode=backend)\n",
"sampler.options.dynamical_decoupling.enable = True\n",
"\n",
"job = sampler.run([isa_circuit], shots=10_000) # Takes approximately 5sec of QPU time"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8c5eb3e9-2a5a-423b-ac18-7bba7f69f18f",
"metadata": {},
"outputs": [],
"source": [
"# Run cell after IQX job completion\n",
"primitive_result = job.result()\n",
"pub_result = primitive_result[0]\n",
"counts = pub_result.data.meas.get_counts()"
]
},
{
"cell_type": "markdown",
"id": "e3ebf77a-8efd-43d6-bd19-e26423921c6f",
"metadata": {},
"source": [
"\n",
"\n",
"## 4. Nachbearbeitungsergebnisse\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "1e7f0ecd",
"metadata": {},
"source": [
"Der Post-Processing-Teil des SQD-Workflows lässt sich anhand des folgenden Diagramms zusammenfassen.\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8c51a527",
"metadata": {},
"source": [
"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:\n",
"\n",
"$$\n",
"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\n",
"$$\n",
"\n",
"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.\n",
"\n",
"$$\n",
"H_{\\mathcal{S}^{(k)}} \\vert \\psi^{(k)} \\rangle = E^{(k)} \\vert \\psi^{(k)} \\rangle\n",
"$$\n",
"\n",
"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.\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "39997b5d-8bd3-4bc8-9ade-11fa5cfb34f9",
"metadata": {},
"source": [
"\n",
"\n",
"### 4.1 Konfigurationswiederherstellung: Übersicht\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "b05881ac-80ce-47a6-ac28-c1237f85bc0a",
"metadata": {},
"source": [
"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.\n",
"\n",
"Die Wahrscheinlichkeit des Flippens kann unter $|x[i] - avg\\_occupancy[i]|$ für `i`-th spin orbital. In [\\[2\\]](#references) verwendeten die Autoren eine gewichtete Umdrehungswahrscheinlichkeit unter Verwendung der modifizierten Funktion ReLU.\n",
"\n",
"$$\n",
"\\begin{align}\n",
" w(y) = \\begin{cases}\n",
"\n",
" \\delta \\frac{y}{h} & \\text{if } y \\leq h\\\\ \\nonumber\n",
"\n",
" \\delta + (1 - \\delta) \\frac{y - h}{1 - h} & \\text{if } y > h\n",
"\n",
"\\end{cases}\n",
"\\end{align}\n",
"$$\n",
"\n",
"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).\n",
"\n",
"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.\n",
"\n",
"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).\n",
"\n",
"Nehmen wir an, dass wir in der ersten Iteration zwei Stapel laufen lassen, und die geschätzten Grundzustände daraus sind:\n",
"\n",
"$$\n",
"\\begin{align}\\nonumber\n",
" \\text{Batch0: } \\vert \\psi \\rangle &= 0.8 \\times \\vert 1001 \\rangle + 0.6 \\times \\vert 0110 \\rangle \\\\ \\nonumber\n",
" \\text{Batch1: } \\vert \\psi \\rangle &= \\frac{1}{\\sqrt{3}} \\left( \\vert 1001 \\rangle + \\vert 0101 \\rangle + \\vert 0110 \\rangle \\right) \\nonumber\n",
"\\end{align}\n",
"$$\n",
"\n",
"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.\n",
"\n",
"Belegung ( Batch0 ):\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :--------------: | :------: | :------: | :------: | :------: |\n",
"| 1001 | 0.64 | 0.0 | 0.0 | 0.64 |\n",
"| 0110 | 0.0 | 0.36 | 0.36 | 0.0 |\n",
"| **n** *(Batch0)* | **0.64** | **0.36** | **0.36** | **0.64** |\n",
"\n",
"Belegung ( Batch1 )\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :--------------: | :------: | :------: | :------: | :------: |\n",
"| 1001 | 0.33 | 0.00 | 0.00 | 0.33 |\n",
"| 0101 | 0.0 | 0.33 | 0.00 | 0.33 |\n",
"| 0110 | 0.0 | 0.33 | 0.33 | 0.00 |\n",
"| **n** *(Batch1)* | **0.33** | **0.66** | **0.33** | **0.66** |\n",
"\n",
"Belegung (Durchschnitt aller Chargen)\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :--------------------: | :------: | :------: | :------: | :------: |\n",
"| **n** *(Batch0)* | 0.64 | 0.36 | 0.36 | 0.64 |\n",
"| **n** *(Batch1)* | 0.33 | 0.66 | 0.33 | 0.66 |\n",
"| **n** *(Durchschnitt)* | **0.49** | **0.51** | **0.35** | **0.65** |\n",
"\n",
"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.\n",
"\n",
"Wahrscheinlichkeit des Umdrehens ( $x = \\vert 1000 \\rangle$, $\\delta = 0.01$, $h = N/M = 2/4 = 0.50$ )\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :------------------------------------------: | :-: | :--: | :---: | :--: |\n",
"| p(flip) ( $\\vert x[i] - \\text{n}[i] \\vert$ ) | 0 | 0.51 | 0.35 | 0.65 |\n",
"| w(p(flip)) | 0 | 0.03 | 0.007 | 0.31 |\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "c1940235",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "5c614290",
"metadata": {},
"source": [
"Der gesamte Prozess der Wiederherstellung einer selbstkonsistenten Konfiguration lässt sich wie folgt zusammenfassen:\n",
"\n",
"**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.\n",
"\n",
"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.\n",
"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$.\n",
"3. Konstruieren Sie aus den angenäherten Eigenzuständen die erste Schätzung für $n$.\n",
"\n",
"**Nachfolgende Iterationen:**\n",
"\n",
"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.\n",
"2. $\\widetilde{\\chi}_{R}$ wird eine Stichprobe gezogen, um Lose zu erstellen $\\mathcal{S}^{(1)}, \\cdots, \\mathcal{S}^{(K)}$.\n",
"3. Der Eigenzustandslöser läuft mit neuen Stapeln und erzeugt neue Schätzungen der Grundzustände $|\\psi^{(1)}\\rangle, \\cdots, |\\psi^{(K)}\\rangle$.\n",
"4. Konstruieren Sie aus den angenäherten Eigenzuständen eine verfeinerte Schätzung für $n$.\n",
"5. Wenn das Abbruchkriterium nicht erfüllt ist, gehen Sie zurück zu Schritt `2.1`.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "5eea548a",
"metadata": {},
"source": [
"\n",
"\n",
"### 4.2 Grundzustandsschätzung\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "31dfe00e",
"metadata": {},
"source": [
"Zunächst werden die Zählungen in eine Bitstring-Matrix und ein Wahrscheinlichkeitsfeld für die Nachbearbeitung umgewandelt.\n",
"\n",
"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.\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "71550274",
"metadata": {},
"outputs": [],
"source": [
"from qiskit_addon_sqd.counts import counts_to_arrays\n",
"\n",
"# Convert counts into bitstring and probability arrays\n",
"bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)"
]
},
{
"cell_type": "markdown",
"id": "acbdbdc5",
"metadata": {},
"source": [
"Es gibt ein paar benutzergesteuerte Optionen, die für diese Technik wichtig sind:\n",
"\n",
"* `iterations`: Anzahl der Iterationen zur Wiederherstellung der selbstkonsistenten Konfiguration\n",
"* `n_batches`: Anzahl der Stapel von Konfigurationen, die von den verschiedenen Aufrufen des Eigenzustandslösers verwendet werden\n",
"* `samples_per_batch`: Anzahl der eindeutigen Konfigurationen, die in jeder Charge enthalten sein sollen\n",
"* `max_davidson_cycles`: Maximale Anzahl von Davidson-Zyklen, die von jedem Eigensolver ausgeführt werden\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a2d22e0b-8f51-42a0-858f-ad0297cb0bae",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting configuration recovery iteration 0\n",
" Batch 0 subspace dimension: 21609\n",
" Batch 1 subspace dimension: 21609\n",
" Batch 2 subspace dimension: 21609\n",
" Batch 3 subspace dimension: 21609\n",
" Batch 4 subspace dimension: 21609\n",
"Starting configuration recovery iteration 1\n",
" Batch 0 subspace dimension: 609961\n",
" Batch 1 subspace dimension: 616225\n",
" Batch 2 subspace dimension: 627264\n",
" Batch 3 subspace dimension: 633616\n",
" Batch 4 subspace dimension: 624100\n",
"Starting configuration recovery iteration 2\n",
" Batch 0 subspace dimension: 564001\n",
" Batch 1 subspace dimension: 605284\n",
" Batch 2 subspace dimension: 582169\n",
" Batch 3 subspace dimension: 559504\n",
" Batch 4 subspace dimension: 591361\n",
"Starting configuration recovery iteration 3\n",
" Batch 0 subspace dimension: 550564\n",
" Batch 1 subspace dimension: 549081\n",
" Batch 2 subspace dimension: 531441\n",
" Batch 3 subspace dimension: 527076\n",
" Batch 4 subspace dimension: 531441\n",
"Starting configuration recovery iteration 4\n",
" Batch 0 subspace dimension: 544644\n",
" Batch 1 subspace dimension: 580644\n",
" Batch 2 subspace dimension: 527076\n",
" Batch 3 subspace dimension: 531441\n",
" Batch 4 subspace dimension: 537289\n"
]
}
],
"source": [
"import numpy as np\n",
"from qiskit_addon_sqd.configuration_recovery import recover_configurations\n",
"from qiskit_addon_sqd.fermion import (\n",
" bitstring_matrix_to_ci_strs,\n",
" solve_fermion,\n",
")\n",
"from qiskit_addon_sqd.subsampling import postselect_and_subsample\n",
"\n",
"rng = np.random.default_rng(24)\n",
"# SQD options\n",
"iterations = 5\n",
"\n",
"# Eigenstate solver options\n",
"n_batches = 5\n",
"samples_per_batch = 500\n",
"max_davidson_cycles = 300\n",
"\n",
"# Self-consistent configuration recovery loop\n",
"e_hist = np.zeros((iterations, n_batches)) # energy history\n",
"s_hist = np.zeros((iterations, n_batches)) # spin history\n",
"occupancy_hist = []\n",
"avg_occupancy = None\n",
"for i in range(iterations):\n",
" print(f\"Starting configuration recovery iteration {i}\")\n",
" # On the first iteration, we have no orbital occupancy information from the\n",
" # solver, so we begin with the full set of noisy configurations.\n",
" if avg_occupancy is None:\n",
" bs_mat_tmp = bitstring_matrix_full\n",
" probs_arr_tmp = probs_arr_full\n",
"\n",
" # If we have average orbital occupancy information, we use it to refine\n",
" # the full set of noisy configurations.\n",
" else:\n",
" bs_mat_tmp, probs_arr_tmp = recover_configurations(\n",
" bitstring_matrix_full,\n",
" probs_arr_full,\n",
" avg_occupancy,\n",
" num_elec_a,\n",
" num_elec_b,\n",
" rand_seed=rng,\n",
" )\n",
"\n",
" # Create batches of subsamples. We postselect here to remove configurations\n",
" # with incorrect hamming weight during iteration 0, since no config recovery was performed.\n",
" batches = postselect_and_subsample(\n",
" bs_mat_tmp,\n",
" probs_arr_tmp,\n",
" hamming_right=num_elec_a,\n",
" hamming_left=num_elec_b,\n",
" samples_per_batch=samples_per_batch,\n",
" num_batches=n_batches,\n",
" rand_seed=rng,\n",
" )\n",
"\n",
" # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.\n",
" e_tmp = np.zeros(n_batches)\n",
" s_tmp = np.zeros(n_batches)\n",
" occs_tmp = []\n",
" coeffs = []\n",
" for j in range(n_batches):\n",
" strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])\n",
" print(f\" Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}\")\n",
" energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(\n",
" batches[j],\n",
" hcore,\n",
" eri,\n",
" open_shell=open_shell,\n",
" spin_sq=spin_sq,\n",
" max_cycle=max_davidson_cycles,\n",
" )\n",
" energy_sci += nuclear_repulsion_energy\n",
" e_tmp[j] = energy_sci\n",
" s_tmp[j] = spin\n",
" occs_tmp.append(avg_occs)\n",
" coeffs.append(coeffs_sci)\n",
"\n",
" # Combine batch results\n",
" avg_occupancy = tuple(np.mean(occs_tmp, axis=0))\n",
"\n",
" # Track optimization history\n",
" e_hist[i, :] = e_tmp\n",
" s_hist[i, :] = s_tmp\n",
" occupancy_hist.append(avg_occupancy)"
]
},
{
"cell_type": "markdown",
"id": "e4f6ec2c-032d-4e18-8ea4-cf867cba5054",
"metadata": {},
"source": [
"\n",
"\n",
"### 4.3 Erörterung der Ergebnisse\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "592eabb7-18f3-4622-8710-bfc43ad6cdec",
"metadata": {},
"source": [
"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.\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "61959e69-a182-4636-abcb-a32349fc9076",
"metadata": {},
"outputs": [],
"source": [
"# Data for energies plot\n",
"x1 = range(iterations)\n",
"min_e = [np.min(e) for e in e_hist]\n",
"e_diff = [abs(e - exact_energy) for e in min_e]\n",
"yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]\n",
"\n",
"# Chemical accuracy (+/- 1 milli-Hartree)\n",
"chem_accuracy = 0.001\n",
"\n",
"# Data for avg spatial orbital occupancy\n",
"y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]\n",
"x2 = range(len(y2))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "8cd90034-6ef3-41bd-a847-c115cade82f7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact energy: -109.04667 Ha\n",
"SQD energy: -109.02234 Ha\n",
"Absolute error: 0.02434 Ha\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n",
"\n",
"# Plot energies\n",
"axs[0].plot(x1, e_diff, label=\"energy error\", marker=\"o\")\n",
"axs[0].set_xticks(x1)\n",
"axs[0].set_xticklabels(x1)\n",
"axs[0].set_yticks(yt1)\n",
"axs[0].set_yticklabels(yt1)\n",
"axs[0].set_yscale(\"log\")\n",
"axs[0].set_ylim(1e-6)\n",
"axs[0].axhline(\n",
" y=chem_accuracy, color=\"#BF5700\", linestyle=\"--\", label=\"chemical accuracy\"\n",
")\n",
"axs[0].set_title(\"Approximated Ground State Energy Error vs SQD Iterations\")\n",
"axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n",
"axs[0].set_ylabel(\"Energy Error (Ha)\", fontdict={\"fontsize\": 12})\n",
"axs[0].legend()\n",
"\n",
"# Plot orbital occupancy\n",
"axs[1].bar(x2, y2, width=0.8)\n",
"axs[1].set_xticks(x2)\n",
"axs[1].set_xticklabels(x2)\n",
"axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n",
"axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n",
"axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n",
"\n",
"print(f\"Exact energy: {exact_energy:.5f} Ha\")\n",
"print(f\"SQD energy: {min_e[-1]:.5f} Ha\")\n",
"print(f\"Absolute error: {e_diff[-1]:.5f} Ha\")\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "e1530472",
"metadata": {},
"source": [
"\n",
"\n",
"#### Übung für den Leser\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d2b2241f",
"metadata": {},
"source": [
"\n",
"\n",
"## Referenzen\n",
"\n",
"\\[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](https://pubs.rsc.org/en/content/articlehtml/2023/sc/d3sc02516k).\n",
"\n",
"\\[2] J. Robledo-Moreno et al, \"Chemie jenseits exakter Lösungen auf einem quantenzentrierten Supercomputer\" (2024). [arXiv:quant-ph/2405.05068](https://arxiv.org/abs/2405.05068).\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"id": "a1b8767d",
"source": "© IBM Corp., 2017-2026"
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}