{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "9e40af77-7f0f-4dd6-ab0a-420cf396050e",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Quantendiagonalisierung eines Chemie-Hamiltonians auf der Grundlage von Stichproben\"\n",
        "description: \"Verwenden Sie den samplebasierten Quanten-Diagonalisierungsalgorithmus, um ein Stickstoffmolekül mit Hilfe von rauschbehafteter Quantenhardware zu simulieren.\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore fontdict fontsize milli LUCJ CCSD ccsd hcore pvdz */}\n",
        "\n",
        "<span id=\"sample-based-quantum-diagonalization-of-a-chemistry-hamiltonian\" />\n",
        "\n",
        "# Quantendiagonalisierung eines Chemie-Hamiltonians auf der Grundlage von Stichproben\n",
        "\n",
        "*Geschätzte Nutzungsdauer: weniger als eine Minute auf einem Heron r2 Prozessor (HINWEIS: Dies ist nur eine Schätzung. Ihre Laufzeit kann variieren.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9701917b",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## Lernziele\n",
        "\n",
        "* Wie man das [SQD-Qiskit-Add-on](/docs/addons/qiskit-addon-sqd) verwendet, um die Grundzustandsenergie eines molekularen Systems anhand von Bitstrings zu approximieren, die von einer Quantenverarbeitungseinheit (QPU) abgetastet wurden.\n",
        "* Wie man mit [ffsim](https://github.com/qiskit-community/ffsim) eine LUCJ-Schaltung (Local Unitary Cluster Jastrow) für die quantenchemische Simulation erstellt.\n",
        "\n",
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## Voraussetzungen\n",
        "\n",
        "* Quantenchemie und zweite Quantisierung\n",
        "* Verwendung der Sampler-Primitive zum Abtasten von Quantenschaltungen\n",
        "\n",
        "<span id=\"background\" />\n",
        "\n",
        "## Hintergrund\n",
        "\n",
        "In diesem Tutorial zeigen wir, wie man verrauschte Quantenstichproben nachbearbeitet, um den Grundzustand des Stickstoffmoleküls $\\text{N}_2$ bei Gleichgewichtsbindungslänge zu approximieren, indem man das [SQD-Qiskit-Add-on](https://github.com/Qiskit/qiskit-addon-sqd) zur Implementierung des [Algorithmus der](https://arxiv.org/abs/2405.05068) stichprobenbasierten Quantendiagonalisierung (SQD) verwendet. Weitere Informationen zur Software finden Sie in der entsprechenden [Dokumentation](/docs/addons/qiskit-addon-sqd), einschließlich eines [einfachen Beispiels](/docs/addons/qiskit-addon-sqd/guides/quickstart) für den Einstieg.\n",
        "\n",
        "Die Stickstoffchemie verbindet dieses Beispiel mit dem Thema Nachhaltigkeit durch die Herstellung von Ammoniak für Düngemittel. Das Verständnis, wie Katalysatoren Stickstoff aktivieren, ist Teil der Forschung zu weniger energieintensiven Produktionsverfahren. Die Berechnung des Grundzustands eines isolierten Stickstoffmoleküls führt zu Methoden der Elektronenstruktur, die für diese umfassendere Herausforderung von Bedeutung sind.\n",
        "\n",
        "Dieses Tutorial wird für Nutzer empfohlen, die mit Quantenchemie vertraut sind, insbesondere mit der Bestimmung der Grundzustandsenergien eines Moleküls. Eine ausführliche Anleitung zum Arbeitsablauf finden Sie im [Kurs](/learning/courses/quantum-diagonalization-algorithms) zum Quanten-Diagonalisierungsalgorithmus.\n",
        "\n",
        "SQD ist eine Methode zur Bestimmung von Eigenwerten und Eigenvektoren von Quantenoperatoren, wie beispielsweise dem Hamilton-Operator eines Quantensystems, durch die kombinierte Nutzung von Quanten- und verteilter klassischer Datenverarbeitung. Im klassischen verteilten Rechnen werden die von einem Quantenprozessor gewonnenen Samples verarbeitet, um einen Ziel-Hamilton-Operator in einen von ihnen aufgespannten Unterraum zu projizieren und zu diagonalisieren. Ein auf SQD basierender Arbeitsablauf umfasst die folgenden Schritte:\n",
        "\n",
        "1. Wählen Sie einen Schaltungsansatz und wenden Sie ihn auf einem Quantencomputer auf einen Referenzzustand (in diesem Fall den [Hartree-Fock-Zustand](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method) ) an.\n",
        "2. Beispiel-Bitstrings aus dem resultierenden Quantenzustand.\n",
        "3. Führen Sie die *selbstkonsistente Konfigurationswiederherstellung* für die Bitstrings durch, um die Grundzustandsnäherung zu erhalten.\n",
        "\n",
        "Es ist bekannt, dass SQD gut funktioniert, wenn der Ziel-Eigenzustand spärlich ist: Die Wellenfunktion ist in einer Menge von Basiszuständen $\\mathcal{S} = \\{|x\\rangle \\}$ enthalten, deren Größe nicht exponentiell mit der Größe des Problems zunimmt.\n",
        "\n",
        "<span id=\"quantum-chemistry\" />\n",
        "\n",
        "### Quantenchemie\n",
        "\n",
        "Der Hamiltonian eines molekularen Systems kann geschrieben werden als\n",
        "\n",
        "$$\n",
        "\\hat{H} = \\sum_{ \\substack{pr\\\\\\sigma} } h_{pr} \\, \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{r\\sigma}\n",
        "+ \\frac12\n",
        "\\sum_{ \\substack{prqs\\\\\\sigma\\tau} }\n",
        "h_{prqs} \\,\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",
        "wobei $h_{pr}$ und $h_{prqs}$ komplexe Zahlen sind, die als molekulare Integrale bezeichnet werden und mit Hilfe eines Computerprogramms aus der Spezifikation des Moleküls berechnet werden können. In diesem Tutorium berechnen wir die Integrale mit Hilfe des [PySCF](https://pyscf.org/) software-Paket.\n",
        "\n",
        "Wie der molekulare Hamiltonian abgeleitet wird, ist in einem Lehrbuch der Quantenchemie nachzulesen (z. B. *Modern Quantum Chemistry* von Szabo und Ostlund). Eine ausführliche Erklärung, wie Probleme der Quantenchemie auf Quantencomputer abgebildet werden, finden Sie in der Vorlesung [*Mapping Problems to Qubits*](https://youtube.com/watch?v=TyFU6r8uEsE\\&t=900) von der Qiskit Global Summer School 2024.\n",
        "\n",
        "<span id=\"local-unitary-cluster-jastrow-lucj-ansatz\" />\n",
        "\n",
        "### Lokaler unitäre Cluster Jastrow (LUCJ) Ansatz\n",
        "\n",
        "SQD benötigt einen Quantenkreis-Ansatz, aus dem Stichproben entnommen werden können. In diesem Tutorial verwenden wir den Ansatz [des](https://pubs.rsc.org/en/content/articlelanding/2023/sc/d3sc02516k) lokalen unitären Clusters Jastrow (LUCJ), da er sowohl physikalisch fundiert als auch hardwarefreundlich ist. Wir werden [ffsim](https://qiskit-community.github.io/ffsim/) verwenden, um die Ansatzschaltung zu erstellen.\n",
        "\n",
        "Der LUCJ-Ansatz lässt sich an QPUs mit eingeschränkter Qubit-Konnektivität anpassen. Die Spinorbitale werden so auf Qubits abgebildet, dass der Ansatz keine Weiterleitung mit SWAP-Gattern erfordert. IBM® Die Hardware weist entweder eine Qubit-Topologie mit quadratischem Gitter oder mit dichtem hexagonalen Gitter auf. Für das Heavy-Hex-Gitter können wir das unten dargestellte „Zick-Zack“-Muster verwenden. In diesem Modell werden Orbitalen mit gleichem Spin auf Qubits mit einer Linientopologie abgebildet (rote und blaue Kreise), und alle vier Orbitalen besteht eine Verbindung zwischen Orbitalen mit unterschiedlichem Spin, wobei diese Verbindung durch ein Ancilla-Qubit (violette Kreise) ermöglicht wird.\n",
        "\n",
        "![Qubit-Abbildungsdiagramm für den LUCJ-Ansatz auf einem schweren Hex-Gitter](https://quantum.cloud.ibm.com/docs/images/tutorials/improving-energy-estimation-of-a-fermionic-hamiltonian-with-sqd/7e0ee7e1-2d24-417f-ac59-25c58db79aa9.avif)\n",
        "\n",
        "<span id=\"self-consistent-configuration-recovery\" />\n",
        "\n",
        "### Selbstkonsistente Konfigurationswiederherstellung\n",
        "\n",
        "Das Verfahren zur Wiederherstellung der selbstkonsistenten Konfiguration ist so konzipiert, dass möglichst viele Signale aus verrauschten Quantenproben extrahiert werden. Da der molekulare Hamiltonian die Teilchenzahl und den Spin Z beibehält, ist es sinnvoll, einen Schaltkreis-Ansatz zu wählen, der diese Symmetrien ebenfalls beibehält. Bei Anwendung auf den Hartree-Fock-Zustand hat der resultierende Zustand eine feste Teilchenzahl und einen festen Spin Z in der geräuschlosen Einstellung. Daher sollten die Spin- $\\alpha$ und Spin- $\\beta$ Hälften jeder Bitkette, die aus diesem Zustand abgetastet wird, das gleiche [Hamming-Gewicht](https://en.wikipedia.org/wiki/Hamming_weight) wie im Hartree-Fock-Zustand haben. Aufgrund des Rauschens in aktuellen Quantenprozessoren werden einige gemessene Bitstrings diese Eigenschaft verletzen. Eine einfache Form der Nachselektion würde diese Bitstrings verwerfen, aber das ist verschwenderisch, weil diese Bitstrings noch ein Signal enthalten könnten. Das Verfahren der selbstkonsistenten Wiederherstellung versucht, einen Teil dieses Signals in der Nachbearbeitung wiederherzustellen. Das Verfahren ist iterativ und erfordert als Eingabe eine Schätzung der durchschnittlichen Besetzungen der einzelnen Orbitale im Grundzustand, die zunächst aus den Rohproben berechnet wird. Das Verfahren wird in einer Schleife ausgeführt, wobei jede Iteration die folgenden Schritte umfasst:\n",
        "\n",
        "1. Für jeden Bitstring, der gegen die angegebenen Symmetrien verstößt, werden die Bits mit einem probabilistischen Verfahren umgedreht, um den Bitstring näher an die aktuelle Schätzung der durchschnittlichen Orbitalbesetzungen heranzuführen und so einen neuen Bitstring zu erhalten.\n",
        "2. Sammeln Sie alle alten und neuen Bitstrings, die die Symmetrien erfüllen, und nehmen Sie eine Teilstichprobe mit einer festen Größe, die im Voraus festgelegt wird.\n",
        "3. Projizieren Sie für jede Teilmenge von Bitstrings den Hamiltonian in den Unterraum, der von den entsprechenden Basisvektoren aufgespannt wird (eine Beschreibung dieser Basisvektoren finden Sie im [vorherigen Abschnitt](#quantum-chemistry) ), und berechnen Sie eine Grundzustandsschätzung des projizierten Hamiltonian auf einem klassischen Computer.\n",
        "4. Aktualisieren Sie die Schätzung der durchschnittlichen Orbitalbesetzungen mit der Grundzustandsschätzung mit der niedrigsten Energie.\n",
        "\n",
        "<span id=\"sqd-workflow-diagram\" />\n",
        "\n",
        "### SQD-Workflow-Diagramm\n",
        "\n",
        "Der SQD-Arbeitsablauf ist im folgenden Diagramm dargestellt:\n",
        "\n",
        "![Arbeitsablaufdiagramm des SQD-Algorithmus](https://quantum.cloud.ibm.com/docs/images/tutorials/improving-energy-estimation-of-a-fermionic-hamiltonian-with-sqd/fd7e816f-4e2e-4dd7-a7da-f71afb9ca68d.avif)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "88422c4b",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## Anforderungen\n",
        "\n",
        "Vergewissern Sie sich vor Beginn dieses Tutorials, dass Sie die folgenden Komponenten installiert haben:\n",
        "\n",
        "* Qiskit SDK v1.0 oder später, mit [Visual](/docs/api/qiskit/visualization) isierungsunterstützung\n",
        "* Qiskit Runtime v0.22 oder später (`pip install qiskit-ibm-runtime`)\n",
        "* SQD-Qiskit-Add-on v0.12 oder höher (`pip install qiskit-addon-sqd`)\n",
        "* ffsim v0.0.75 oder höher (`pip install ffsim`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c6e44a31",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## Einrichtung\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "6e51c3d8",
      "metadata": {},
      "outputs": [],
      "source": [
        "import math\n",
        "\n",
        "import ffsim\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import pyscf\n",
        "import pyscf.cc\n",
        "import pyscf.mcscf\n",
        "from qiskit import QuantumCircuit, QuantumRegister\n",
        "from qiskit.primitives import StatevectorSampler\n",
        "from qiskit.providers.fake_provider import GenericBackendV2\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "from qiskit_ibm_runtime import SamplerV2 as Sampler"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4bc6ee26-4371-4cd2-80a7-60752bf8775d",
      "metadata": {},
      "source": [
        "<span id=\"small-scale-simulator-example\" />\n",
        "\n",
        "## Beispiel für einen Simulator im kleinen Maßstab\n",
        "\n",
        "In diesem Tutorial werden wir eine Näherung für den Grundzustand eines Stickstoffmoleküls in der Nähe seines Gleichgewichtsbindungsabstands finden. Zunächst verwenden wir einen kleinen Basis-Satz von „ STO-6G “, um das Experiment zu simulieren und sicherzustellen, dass es funktioniert.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "afeb054c",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### Schritt 1: Klassische Eingaben einem Quantenproblem zuordnen\n",
        "\n",
        "Zunächst legen wir das Molekül und seine Eigenschaften fest.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "b821e660",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "converged SCF energy = -108.464957764796\n",
            "CASCI E = -108.595987350986  E(CI) = -32.4115475088426  S^2 = 0.0000000\n",
            "norb = 8\n",
            "nelec = (5, 5)\n"
          ]
        }
      ],
      "source": [
        "# Specify molecule properties\n",
        "spin_sq = 0\n",
        "\n",
        "# Build N2 molecule\n",
        "mol = pyscf.gto.Mole()\n",
        "mol.build(\n",
        "    atom=[[\"N\", (0, 0, -0.5)], [\"N\", (0, 0, 0.5)]],\n",
        "    basis=\"sto-6g\",\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",
        "norb = len(active_space)\n",
        "n_electrons = int(sum(scf.mo_occ[active_space]))\n",
        "n_alpha = (n_electrons + mol.spin) // 2\n",
        "n_beta = (n_electrons - mol.spin) // 2\n",
        "nelec = (n_alpha, n_beta)\n",
        "cas = pyscf.mcscf.CASCI(scf, norb, nelec)\n",
        "mo = cas.sort_mo(active_space, base=0)\n",
        "hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n",
        "eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)\n",
        "\n",
        "# Compute exact energy using FCI\n",
        "reference_energy = cas.run().e_tot\n",
        "\n",
        "print(f\"norb = {norb}\")\n",
        "print(f\"nelec = {nelec}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "96bfe018",
      "metadata": {},
      "source": [
        "Bevor wir die LUCJ-Ansatzschaltung konstruieren, führen wir zunächst eine CCSD-Berechnung in der folgenden Codezelle durch. Die [$t_1$ und $t_2$ Amplituden](https://en.wikipedia.org/wiki/Coupled_cluster#Cluster_operator) aus dieser Berechnung werden zur Initialisierung der Parameter des Ansatzes verwendet.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "efe83d98",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "E(CCSD) = -108.5933309085008  E_corr = -0.1283731437052352\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",
        ").run()\n",
        "t1 = ccsd.t1\n",
        "t2 = ccsd.t2"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "f4d882fa",
      "metadata": {},
      "source": [
        "Nun verwenden wir [ffsim](https://github.com/qiskit-community/ffsim), um das Ansatz-Schaltbild zu erstellen. Da unser Molekül einen Hartree-Fock-Zustand mit geschlossener Schale aufweist, verwenden wir die spinausgeglichene Variante des UCJ-Ansatzes, [UCJOpSpinBalanced](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced). Wir setzen `optimize=True` in der `from_t_amplitudes` Methode, um die „komprimierte“ Doppelfaktorisierung der Amplituden von $t_2$ s zu aktivieren (Einzelheiten finden Sie [im Abschnitt](https://qiskit-community.github.io/ffsim/explanations/lucj.html#Parameter-initialization-from-CCSD) „The local unitary cluster Jastrow (LUCJ) ansatz“ in der Dokumentation von ffsim).\n",
        "\n",
        "Da sich der LUCJ-Ansatz an die verfügbare Konnektivität der QPU anpasst, müssen wir das QPU-Backend initialisieren, bevor wir den Ansatz erstellen. Zunächst erstellen wir ein generisches Backend mit einer quadratischen Gitterkopplungsabbildung und einem Gatter-Set, in das sich der LUCJ-Ansatz auf natürliche Weise zerlegen lässt. Anschließend verwenden wir, `ffsim.qiskit.generate_lucj_pass_manager` um einen Pass-Manager zu erstellen, der darauf spezialisiert ist, den LUCJ-Ansatz gemäß einem „Ladder“-Layout (siehe die [Erläuterung](https://qiskit-community.github.io/ffsim/explanations/lucj.html#The-local-UCJ-\\(LUCJ\\)-ansatz) in der Dokumentation von ffsim) in das angegebene Backend zu transpilieren. Diese Funktion nutzt eine Bewertungsheuristik, um die mit dem ausgewählten Layout verbundenen Fehler zu minimieren. Dies ist wichtig, wenn es sich bei Ihrem Backend um eine echte QPU oder um einen Simulator mit einem Rauschmodell handelt. Diese Funktion gibt nicht nur den Pass-Manager zurück, sondern auch die Alpha-Beta-Kopplungspaare, die auf der Hardware implementiert werden können. Wenn nicht alle Paare implementiert werden können, wird eine Warnung ausgegeben.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "dd69a86c",
      "metadata": {},
      "outputs": [],
      "source": [
        "import warnings\n",
        "\n",
        "from qiskit.transpiler import CouplingMap\n",
        "\n",
        "warnings.formatwarning = lambda msg, *args, **kwargs: f\"Warning: {msg}\\n\"\n",
        "\n",
        "# Set ansatz properties\n",
        "n_reps = 4\n",
        "pairs_aa = [(p, p + 1) for p in range(norb - 1)]\n",
        "\n",
        "# Let generate_lucj_pass_manager determine the alpha-beta interactions\n",
        "pairs_ab = None\n",
        "\n",
        "# Initialize backend\n",
        "coupling_map = CouplingMap.from_grid(norb, 2)\n",
        "backend = GenericBackendV2(\n",
        "    coupling_map.size(),\n",
        "    coupling_map=coupling_map,\n",
        "    basis_gates=[\"cp\", \"xx_plus_yy\", \"p\", \"x\", \"swap\"],\n",
        ")\n",
        "\n",
        "# Create pass manager\n",
        "pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(\n",
        "    backend=backend,\n",
        "    norb=norb,\n",
        "    connectivity=\"square\",\n",
        "    interaction_pairs=(pairs_aa, pairs_ab),\n",
        "    optimization_level=3,\n",
        ")\n",
        "\n",
        "# Create the LUCJ ansatz operator\n",
        "ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
        "    t2=t2,\n",
        "    t1=t1,\n",
        "    n_reps=n_reps,\n",
        "    interaction_pairs=(pairs_aa, pairs_ab),\n",
        "    # Setting optimize=True enables the \"compressed\" factorization\n",
        "    optimize=True,\n",
        ")\n",
        "\n",
        "# create an empty quantum circuit\n",
        "qubits = QuantumRegister(2 * norb, name=\"q\")\n",
        "circuit = QuantumCircuit(qubits)\n",
        "\n",
        "# prepare Hartree-Fock state as the reference state and append it\n",
        "# to the quantum circuit\n",
        "circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, 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()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "db11bf6d",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-for-quantum-hardware-execution\" />\n",
        "\n",
        "### Schritt 2: Optimierung für die Ausführung auf Quantenhardware\n",
        "\n",
        "Als Nächstes optimieren wir die Schaltung für eine bestimmte Hardware. In der Regel umfasst dieser Schritt die Initialisierung des Hardware-Backends und eines Pass-Managers für dieses Backend. Da der LUCJ-Ansatz jedoch an die Hardware-Konnektivität angepasst ist, haben wir diese Schritte bereits im vorherigen Schritt durchgeführt. Nun muss nur noch der Pass-Manager auf der Schaltung ausgeführt werden, um sie in eine ISA-Schaltung zu transpilieren, die direkt auf der QPU ausgeführt werden kann.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "7d554aa5",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Gate counts: OrderedDict({'xx_plus_yy': 254, 'cp': 88, 'p': 64, 'measure': 16, 'x': 10, 'barrier': 1})\n"
          ]
        }
      ],
      "source": [
        "isa_circuit = pass_manager.run(circuit)\n",
        "print(f\"Gate counts: {isa_circuit.count_ops()}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0cc1edef",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-qiskit-primitives\" />\n",
        "\n",
        "### Schritt 3: Mit „ Qiskit primitives “ ausführen\n",
        "\n",
        "Nachdem wir die Schaltung für die Hardware-Ausführung optimiert haben, können wir sie nun auf der Zielhardware ausführen und Messwerte für die Schätzung der Grundzustandsenergie erfassen. Da wir nur einen Circuit haben, werden wir [den Ausführungsmodus](/docs/guides/execution-modes) „ IBM Quantum “ des Compute Service verwenden und unseren Circuit ausführen.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "93c1cef3-298e-4deb-8512-769fe94cd5a5",
      "metadata": {},
      "outputs": [],
      "source": [
        "rng = np.random.default_rng(12345)\n",
        "sampler = StatevectorSampler(seed=rng)\n",
        "job = sampler.run([isa_circuit], shots=100_000)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "332ecab3-77e6-473f-b0e7-af30f983393a",
      "metadata": {},
      "outputs": [],
      "source": [
        "primitive_result = job.result()\n",
        "pub_result = primitive_result[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6df05b6e",
      "metadata": {},
      "source": [
        "<span id=\"step-4-post-process-and-return-result-in-desired-classical-format\" />\n",
        "\n",
        "### Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format\n",
        "\n",
        "Eine nützliche Kennzahl zur Beurteilung der Qualität der QPU-Ausgabe ist die Anzahl der zurückgegebenen gültigen Konfigurationen. Eine gültige Konfiguration hat die richtige Partikelanzahl und den richtigen Spin Z, was bedeutet, dass die rechte Hälfte der Bitfolge ein Hamming-Gewicht hat, das der Anzahl der Spin-up-Elektronen entspricht, und die linke Hälfte ein Hamming-Gewicht hat, das der Anzahl der Spin-down-Elektronen entspricht. Die folgende Zelle berechnet den Anteil der gültigen Konfigurationen aus der Stichprobe.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "718f8517",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Fraction of sampled configurations that are valid: 1.0\n"
          ]
        }
      ],
      "source": [
        "def is_valid_bitstring(\n",
        "    bitstring: str, norb: int, nelec: tuple[int, int]\n",
        ") -> bool:\n",
        "    n_alpha, n_beta = nelec\n",
        "    return (\n",
        "        len(bitstring) == 2 * norb\n",
        "        and bitstring[norb:].count(\"1\") == n_alpha\n",
        "        and bitstring[:norb].count(\"1\") == n_beta\n",
        "    )\n",
        "\n",
        "\n",
        "bit_array = pub_result.data.meas\n",
        "num_valid = sum(\n",
        "    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n",
        ")\n",
        "valid_fraction = num_valid / bit_array.num_shots\n",
        "print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7b126f3a",
      "metadata": {},
      "source": [
        "Alle Bitstrings sind gültig, da wir die Schaltung auf einem rauschfreien Simulator abtasten. Bei Ausführung auf einer verrauschten QPU wird dieser Anteil kleiner als eins sein, aber hoffentlich größer als der Anteil, den man erwarten würde, wenn die Bitstrings gleichmäßig zufällig ausgewählt würden – dieser Wert wird in der folgenden Zelle berechnet.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "6b3e4bca",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625\n"
          ]
        }
      ],
      "source": [
        "expected_fraction_random = (\n",
        "    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n",
        ")\n",
        "print(\n",
        "    f\"Expected fraction of valid configurations from uniformly random bitstrings: \"\n",
        "    f\"{expected_fraction_random}\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "eb704101-0fe8-4d12-b572-b1d844e35a90",
      "metadata": {},
      "source": [
        "Nun schätzen wir die Grundzustandsenergie des Hamiltonianers mit Hilfe der Funktion `diagonalize_fermionic_hamiltonian` . Diese Funktion führt das Verfahren der selbstkonsistenten Konfigurationswiederherstellung durch, um die verrauschten Quantenproben iterativ zu verfeinern und die Energieschätzung zu verbessern. Wir übergeben eine Callback-Funktion, damit wir die Zwischenergebnisse für eine spätere Analyse speichern können. In der [API-Dokumentation](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) finden Sie Erläuterungen zu den Argumenten für `diagonalize_fermionic_hamiltonian`.\n",
        "\n",
        "Hier verwenden wir das `initial_occupancies` Argument, `diagonalize_fermionic_hamiltonian` um die Hartree-Fock-Konfiguration als anfängliche Schätzung für die Orbitalbelegungen im Grundzustand anzugeben. Dieser Ansatz ist sinnvoll für Systeme, bei denen der Grundzustand eine signifikante Unterstützung für die Hartree-Fock-Konfiguration aufweist, aber in anderen Situationen möglicherweise nicht geeignet ist, obwohl fortgeschrittenere Berechnungsmethoden in diesen Fällen zu besseren Anfangsschätzungen führen könnten. Durch die Angabe von kann die `initial_occupancies` Konfigurationswiederherstellung auch dann ausgeführt werden, wenn keine gültigen Konfigurationen erfasst wurden, wie dies beispielsweise bei der Erfassung eines großen Schaltkreises auf einer verrauschten QPU der Fall sein kann. Ohne dieses Argument würde die Konfigurationswiederherstellung fehlschlagen und einen Fehler auslösen, wenn keine gültigen Konfigurationen angegeben wurden.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "2f32a352",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -108.59589540253083\n",
            "\t\tSubspace dimension: 2500\n",
            "Final energy: -108.59589540253083\n",
            "Final energy error: 9.19484552071026e-05\n"
          ]
        }
      ],
      "source": [
        "from functools import partial\n",
        "\n",
        "from qiskit_addon_sqd.fermion import (\n",
        "    SCIResult,\n",
        "    diagonalize_fermionic_hamiltonian,\n",
        "    solve_sci_batch,\n",
        ")\n",
        "\n",
        "# SQD options\n",
        "energy_tol = 1e-3\n",
        "occupancies_tol = 1e-3\n",
        "max_iterations = 5\n",
        "\n",
        "# Eigenstate solver options\n",
        "num_batches = 3\n",
        "samples_per_batch = 1000\n",
        "max_dim = 50\n",
        "symmetrize_spin = True\n",
        "carryover_threshold = 1e-4\n",
        "max_cycle = 200\n",
        "\n",
        "# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies\n",
        "initial_occupancies = (\n",
        "    np.array([1] * n_alpha + [0] * (norb - n_alpha)),\n",
        "    np.array([1] * n_beta + [0] * (norb - n_beta)),\n",
        ")\n",
        "\n",
        "# Pass options to the built-in eigensolver. If you just want to use the defaults,\n",
        "# you can omit this step, in which case you would not specify the sci_solver argument\n",
        "# in the call to diagonalize_fermionic_hamiltonian below.\n",
        "sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)\n",
        "\n",
        "# List to capture intermediate results\n",
        "result_history = []\n",
        "\n",
        "\n",
        "def callback(results: list[SCIResult]):\n",
        "    result_history.append(results)\n",
        "    iteration = len(result_history)\n",
        "    print(f\"Iteration {iteration}\")\n",
        "    for i, result in enumerate(results):\n",
        "        print(f\"\\tSubsample {i}\")\n",
        "        print(f\"\\t\\tEnergy: {result.energy + nuclear_repulsion_energy}\")\n",
        "        print(\n",
        "            f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
        "        )\n",
        "\n",
        "\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    hcore,\n",
        "    eri,\n",
        "    bit_array,\n",
        "    samples_per_batch=samples_per_batch,\n",
        "    norb=norb,\n",
        "    nelec=nelec,\n",
        "    num_batches=num_batches,\n",
        "    energy_tol=energy_tol,\n",
        "    occupancies_tol=occupancies_tol,\n",
        "    max_iterations=max_iterations,\n",
        "    sci_solver=sci_solver,\n",
        "    symmetrize_spin=symmetrize_spin,\n",
        "    max_dim=max_dim,\n",
        "    initial_occupancies=initial_occupancies,\n",
        "    carryover_threshold=carryover_threshold,\n",
        "    callback=callback,\n",
        "    seed=rng,\n",
        ")\n",
        "\n",
        "final_energy = result.energy + nuclear_repulsion_energy\n",
        "energy_error = final_energy - reference_energy\n",
        "print(f\"Final energy: {final_energy}\")\n",
        "print(f\"Final energy error: {energy_error}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9d78906b-4759-4506-9c69-85d4e67766b3",
      "metadata": {},
      "source": [
        "<span id=\"visualize-the-results\" />\n",
        "\n",
        "#### Visualisieren Sie die Ergebnisse\n",
        "\n",
        "Die erste Grafik zeigt, dass wir in dieser Simulation bereits nach der ersten Iteration nahe `1 mH` am exakten Ergebnis liegen (die chemische Genauigkeit wird üblicherweise mit `1 kcal/mol`$\\approx$ angegeben `1.6 mH`). Es handelt sich hierbei jedoch um ein kleines System, und da die Stichproben rauschfrei sind, ist eine Konfigurationswiederherstellung nicht erforderlich. Bei einem größeren System, das auf einer rauschbehafteten QPU läuft, sind möglicherweise mehrere Iterationen zur Konfigurationswiederherstellung erforderlich, und die endgültige Genauigkeit könnte geringer ausfallen. Im Allgemeinen lässt sich die Energie verbessern, indem man mehr Iterationen zur Konfigurationswiederherstellung zulässt oder die Anzahl der Stichproben pro Stapel erhöht.\n",
        "\n",
        "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"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "caffd888-e89c-4aa9-8bae-4d1bb723b35e",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/caffd888-e89c-4aa9-8bae-4d1bb723b35e-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Data for energies plot\n",
        "x1 = range(len(result_history))\n",
        "min_e = [\n",
        "    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy\n",
        "    for result in result_history\n",
        "]\n",
        "e_diff = [abs(e - reference_energy) for e in min_e]\n",
        "\n",
        "# Chemical accuracy (+/- 1.6 milli-Hartree)\n",
        "chem_accuracy = 0.0016\n",
        "\n",
        "# Data for avg spatial orbital occupancy\n",
        "y2 = np.sum(result.orbital_occupancies, axis=0)\n",
        "x2 = range(len(y2))\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_yscale(\"log\")\n",
        "axs[0].set_ylim(1e-5, 1e-2)\n",
        "axs[0].axhline(\n",
        "    y=chem_accuracy,\n",
        "    color=\"#BF5700\",\n",
        "    linestyle=\"--\",\n",
        "    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",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ce0eecb3-8a23-4118-aa1e-a28afcec6334",
      "metadata": {},
      "source": [
        "<span id=\"large-scale-hardware-example\" />\n",
        "\n",
        "## Beispiel für groß angelegte Hardware\n",
        "\n",
        "Nun führen wir ein umfangreicheres Beispiel auf echter Quantenhardware durch. Hier leiten wir einen aktiven Raum für das Stickstoffmolekül aus dem Basissatz „ cc-pVDZ “ ab. Außerdem werden wir folgende Änderungen vornehmen:\n",
        "\n",
        "* Verringern Sie den Wert `n_reps` , `1` um Störgeräusche durch Gate-Fehler zu reduzieren\n",
        "* Stellen Sie die Ansatz-Konnektivität entsprechend der Backend-Prozessorfamilie ein („heavy-hex“ für Heron, „square“ für Nighthawk)\n",
        "* Erhöhen Sie auf `max_dim` , `1000` um größere Unterraumdimensionen zu ermöglichen.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "24ca3090-3f3b-4efb-a482-70b2e1b5d062",
      "metadata": {},
      "source": [
        "<span id=\"steps-1-4\" />\n",
        "\n",
        "### Schritte 1–4\n",
        "\n",
        "Hier fassen wir alle Schritte zu einem einzigen, groß angelegten Arbeitsablauf zusammen, der dann auf echter Quantenhardware ausgeführt wird.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "3858949c-a55d-4ff8-a0fc-54fb53e131b5",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "converged SCF energy = -108.929838385609\n",
            "norb = 26\n",
            "nelec = (5, 5)\n",
            "E(CCSD) = -109.2177884185544  E_corr = -0.2879500329450041\n",
            "Using backend ibm_boston\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].\n",
            "Removing interaction (24, 24) from the end.\n",
            "Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].\n",
            "Removing interaction (20, 20) from the end.\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Gate counts: OrderedDict({'sx': 7036, 'rz': 6965, 'cz': 1858, 'measure': 52, 'x': 50, 'barrier': 1})\n",
            "Fraction of sampled configurations that are valid: 0.01677\n",
            "Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07\n",
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.17922180754668\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.17713541195594\n",
            "\t\tSubspace dimension: 962361\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.17407767345108\n",
            "\t\tSubspace dimension: 968256\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.1873108379026\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.1876970140745\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.18255107152576\n",
            "\t\tSubspace dimension: 1000000\n",
            "Iteration 3\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.19052148187927\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.19316375313544\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.19087348337145\n",
            "\t\tSubspace dimension: 1000000\n",
            "Iteration 4\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.19374847699348\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.19517039328333\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.19422416032815\n",
            "\t\tSubspace dimension: 1000000\n",
            "Iteration 5\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.19580768221108\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.19608553581007\n",
            "\t\tSubspace dimension: 1000000\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.19548131833729\n",
            "\t\tSubspace dimension: 1000000\n",
            "Final energy: -109.19608553581007\n",
            "Final energy error: 0.03194368084709254\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/sample-based-quantum-diagonalization/extracted-outputs/3858949c-a55d-4ff8-a0fc-54fb53e131b5-3.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# ------------------------------ Step 1 ------------------------------\n",
        "# Build N2 molecule\n",
        "mol = pyscf.gto.Mole()\n",
        "mol.build(\n",
        "    atom=[[\"N\", (0, 0, -0.5)], [\"N\", (0, 0, 0.5)]],\n",
        "    basis=\"cc-pvdz\",\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",
        "norb = len(active_space)\n",
        "n_electrons = int(sum(scf.mo_occ[active_space]))\n",
        "n_alpha = (n_electrons + mol.spin) // 2\n",
        "n_beta = (n_electrons - mol.spin) // 2\n",
        "nelec = (n_alpha, n_beta)\n",
        "cas = pyscf.mcscf.CASCI(scf, norb, nelec)\n",
        "mo = cas.sort_mo(active_space, base=0)\n",
        "hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n",
        "eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)\n",
        "\n",
        "# Store reference energy from SCI calculation performed separately\n",
        "reference_energy = -109.22802921665716\n",
        "\n",
        "print(f\"norb = {norb}\")\n",
        "print(f\"nelec = {nelec}\")\n",
        "\n",
        "# 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",
        ").run()\n",
        "t1 = ccsd.t1\n",
        "t2 = ccsd.t2\n",
        "\n",
        "# Set ansatz properties\n",
        "n_reps = 1\n",
        "pairs_aa = [(p, p + 1) for p in range(norb - 1)]\n",
        "\n",
        "# Let generate_lucj_pass_manager determine the alpha-beta interactions\n",
        "pairs_ab = None\n",
        "\n",
        "# Initialize backend\n",
        "service = QiskitRuntimeService()\n",
        "backend = service.least_busy(operational=True, simulator=False)\n",
        "print(f\"Using backend {backend.name}\")\n",
        "\n",
        "# Create pass manager\n",
        "connectivity_map = {\"Heron\": \"heavy-hex\", \"Nighthawk\": \"square\"}\n",
        "pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(\n",
        "    backend=backend,\n",
        "    norb=norb,\n",
        "    connectivity=connectivity_map[\n",
        "        backend.configuration().processor_type[\"family\"]\n",
        "    ],\n",
        "    interaction_pairs=(pairs_aa, pairs_ab),\n",
        "    optimization_level=3,\n",
        ")\n",
        "\n",
        "# Create the LUCJ ansatz operator\n",
        "ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
        "    t2=t2,\n",
        "    t1=t1,\n",
        "    n_reps=n_reps,\n",
        "    interaction_pairs=(pairs_aa, pairs_ab),\n",
        "    # Setting optimize=True enables the \"compressed\" factorization\n",
        "    optimize=True,\n",
        ")\n",
        "\n",
        "# create an empty quantum circuit\n",
        "qubits = QuantumRegister(2 * norb, name=\"q\")\n",
        "circuit = QuantumCircuit(qubits)\n",
        "\n",
        "# prepare Hartree-Fock state as the reference state and append it\n",
        "# to the quantum circuit\n",
        "circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, 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",
        "\n",
        "\n",
        "# ------------------------------ Step 2 ------------------------------\n",
        "\n",
        "isa_circuit = pass_manager.run(circuit)\n",
        "print(f\"Gate counts: {isa_circuit.count_ops()}\")\n",
        "\n",
        "\n",
        "# ------------------------------ Step 3 ------------------------------\n",
        "sampler = Sampler(mode=backend)\n",
        "sampler.options.environment.job_tags = [\"TUT_SQD\"]\n",
        "job = sampler.run([isa_circuit], shots=100_000)\n",
        "primitive_result = job.result()\n",
        "pub_result = primitive_result[0]\n",
        "\n",
        "\n",
        "# ------------------------------ Step 4 ------------------------------\n",
        "\n",
        "bit_array = pub_result.data.meas\n",
        "num_valid = sum(\n",
        "    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n",
        ")\n",
        "valid_fraction = num_valid / bit_array.num_shots\n",
        "print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")\n",
        "expected_fraction_random = (\n",
        "    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n",
        ")\n",
        "print(\n",
        "    f\"Expected fraction of valid configurations from uniformly random bitstrings: \"\n",
        "    f\"{expected_fraction_random}\"\n",
        ")\n",
        "# SQD options\n",
        "energy_tol = 1e-3\n",
        "occupancies_tol = 1e-3\n",
        "max_iterations = 5\n",
        "\n",
        "# Eigenstate solver options\n",
        "num_batches = 3\n",
        "samples_per_batch = 1000\n",
        "max_dim = 1000\n",
        "symmetrize_spin = True\n",
        "carryover_threshold = 1e-4\n",
        "max_cycle = 200\n",
        "\n",
        "# Use the Hartree-Fock configuration as an initial guess for the\n",
        "# orbital occupancies\n",
        "initial_occupancies = (\n",
        "    np.array([1] * n_alpha + [0] * (norb - n_alpha)),\n",
        "    np.array([1] * n_beta + [0] * (norb - n_beta)),\n",
        ")\n",
        "\n",
        "# Pass options to the built-in eigensolver. If you just want to use the defaults,\n",
        "# you can omit this step, in which case you would not specify the\n",
        "# sci_solver argument in the call to diagonalize_fermionic_hamiltonian below.\n",
        "sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)\n",
        "\n",
        "# List to capture intermediate results\n",
        "result_history = []\n",
        "\n",
        "\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    hcore,\n",
        "    eri,\n",
        "    bit_array,\n",
        "    samples_per_batch=samples_per_batch,\n",
        "    norb=norb,\n",
        "    nelec=nelec,\n",
        "    num_batches=num_batches,\n",
        "    energy_tol=energy_tol,\n",
        "    occupancies_tol=occupancies_tol,\n",
        "    max_iterations=max_iterations,\n",
        "    sci_solver=sci_solver,\n",
        "    symmetrize_spin=symmetrize_spin,\n",
        "    max_dim=max_dim,\n",
        "    initial_occupancies=initial_occupancies,\n",
        "    carryover_threshold=carryover_threshold,\n",
        "    callback=callback,\n",
        "    seed=rng,\n",
        ")\n",
        "\n",
        "final_energy = result.energy + nuclear_repulsion_energy\n",
        "energy_error = final_energy - reference_energy\n",
        "print(f\"Final energy: {final_energy}\")\n",
        "print(f\"Final energy error: {energy_error}\")\n",
        "\n",
        "# Data for energies plot\n",
        "x1 = range(len(result_history))\n",
        "min_e = [\n",
        "    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy\n",
        "    for result in result_history\n",
        "]\n",
        "e_diff = [abs(e - reference_energy) for e in min_e]\n",
        "\n",
        "# Chemical accuracy (+/- 1.6 milli-Hartree)\n",
        "chem_accuracy = 0.0016\n",
        "\n",
        "# Data for avg spatial orbital occupancy\n",
        "y2 = np.sum(result.orbital_occupancies, axis=0)\n",
        "x2 = range(len(y2))\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_yscale(\"log\")\n",
        "axs[0].set_ylim(1e-3, 1e-1)\n",
        "axs[0].axhline(\n",
        "    y=chem_accuracy,\n",
        "    color=\"#BF5700\",\n",
        "    linestyle=\"--\",\n",
        "    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",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "405e89ea-57da-4021-bb18-91e8d583d310",
      "metadata": {},
      "source": [
        "<span id=\"next-steps\" />\n",
        "\n",
        "## Nächste Schritte\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Empfehlungen\">\n",
        "  Wenn Sie diese Arbeit interessant fanden, könnten Sie auch an folgenden Materialien interessiert sein:\n",
        "\n",
        "  * [Beispielbasierte Krylov-Quanten-Diagonalisierung eines fermionischen Gittermodells](/docs/tutorials/sample-based-krylov-quantum-diagonalization) – ein dazugehöriges Tutorial, bei dem Zeitentwicklungsschaltungen anstelle eines variationalen Ansatzes verwendet werden\n",
        "  * [Skalierung von SQD-Chemie-Workflows mit dem Dice-Solver](/docs/addons/qiskit-addon-sqd/guides/integrate-dice-solver) – eine Seite, die zeigt, wie man die effizientere Dice-Software zur Diagonalisierung einsetzt\n",
        "  * `diagonalize_fermionic_hamiltonian`[ SQD-Add-on-API-Dokumentation](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) – Referenz zur Funktion\n",
        "  * [*Chemie jenseits der Grenzen der exakten Diagonalisierung auf einem quantenzentrierten Supercomputer*](https://www.science.org/doi/10.1126/sciadv.adu9991) – der Artikel, auf dem dieses Tutorial basiert\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "hours": 1.5,
    "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"
    },
    "qpuSeconds": 60
  },
  "nbformat": 4,
  "nbformat_minor": 5
}