{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "9e40af77-7f0f-4dd6-ab0a-420cf396050e",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Diagonalisation quantique basée sur l'échantillonnage d'un hamiltonien chimique\"\n",
        "description: \"Utilisez l'algorithme de diagonalisation quantique basé sur des échantillons pour simuler une molécule d'azote à l'aide d'un matériel quantique bruité.\"\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",
        "# Diagonalisation quantique basée sur l'échantillonnage d'un hamiltonien chimique\n",
        "\n",
        "*Estimation de l'utilisation : moins d'une minute sur un processeur Heron r2 (NOTE : Il s'agit uniquement d'une estimation. Votre durée d'exécution peut varier.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9701917b",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## Résultats d'apprentissage\n",
        "\n",
        "À l'issue de ce tutoriel, les utilisateurs devraient avoir compris :\n",
        "\n",
        "* Comment utiliser le [module complémentaire SQD Qiskit](/docs/addons/qiskit-addon-sqd) pour estimer l'énergie de l'état fondamental d'un système moléculaire à l'aide de chaînes de bits échantillonnées à partir d'une unité de traitement quantique (QPU).\n",
        "* Comment utiliser [ffsim](https://github.com/qiskit-community/ffsim) pour construire un circuit Jastrow à grappes unitaires locales (LUCJ) destiné à la simulation en chimie quantique.\n",
        "\n",
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## Prérequis\n",
        "\n",
        "Nous recommandons aux utilisateurs de se familiariser avec les sujets suivants avant de suivre ce tutoriel :\n",
        "\n",
        "* Chimie quantique et seconde quantification\n",
        "* Utilisation de la primitive « Sampler » pour échantillonner des circuits quantiques\n",
        "\n",
        "<span id=\"background\" />\n",
        "\n",
        "## Arrière-plan\n",
        "\n",
        "Dans ce tutoriel, nous montrons comment traiter des échantillons quantiques bruités afin d'approximer l'état fondamental de la molécule d'azote $\\text{N}_2$ à la longueur de liaison d'équilibre, en utilisant [l](https://github.com/Qiskit/qiskit-addon-sqd) 'extension SQD de Qiskit pour mettre en œuvre [l](https://arxiv.org/abs/2405.05068) 'algorithme de diagonalisation quantique par échantillonnage (SQD). Vous trouverez plus de détails sur ce logiciel dans la [documentation](/docs/addons/qiskit-addon-sqd) correspondante, qui comprend notamment un [exemple simple](/docs/addons/qiskit-addon-sqd/guides/quickstart) pour vous aider à démarrer.\n",
        "\n",
        "Ce tutoriel est recommandé aux utilisateurs familiarisés avec la chimie quantique, et plus particulièrement à ceux qui savent déterminer l'énergie de l'état fondamental d'une molécule. Pour un guide détaillé sur le déroulement du processus, consultez le [cours](/learning/courses/quantum-diagonalization-algorithms) sur l'algorithme de diagonalisation quantique.\n",
        "\n",
        "La SQD est une technique permettant de déterminer les valeurs propres et les vecteurs propres d'opérateurs quantiques, tels que l'hamiltonien d'un système quantique, en combinant le calcul quantique et le calcul classique distribué. Le calcul distribué classique sert à traiter les échantillons obtenus à partir d'un processeur quantique, ainsi qu'à projeter et à diagonaliser un hamiltonien cible dans un sous-espace qu'ils génèrent. Un processus basé sur le SQD comprend les étapes suivantes :\n",
        "\n",
        "1. Choisissez un ansatz de circuit et appliquez-le sur un ordinateur quantique à un état de référence (dans ce cas, l'état [Hartree-Fock](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method) ).\n",
        "2. Echantillons de chaînes de bits de l'état quantique résultant.\n",
        "3. Exécutez la procédure *de récupération de configuration auto-cohérente* sur les chaînes de bits afin d'obtenir l'approximation de l'état fondamental.\n",
        "\n",
        "On sait que la SQD fonctionne bien lorsque l'état propre cible est peu dense : la fonction d'onde est supportée par un ensemble d'états de base $\\mathcal{S} = \\{|x\\rangle \\}$ dont la taille n'augmente pas de manière exponentielle avec la taille du problème.\n",
        "\n",
        "<span id=\"quantum-chemistry\" />\n",
        "\n",
        "### Chimie quantique\n",
        "\n",
        "L'hamiltonien d'un système moléculaire peut être écrit comme suit\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",
        "où $h_{pr}$ et $h_{prqs}$ sont des nombres complexes appelés intégrales moléculaires qui peuvent être calculées à partir des spécifications de la molécule à l'aide d'un programme informatique. Dans ce tutoriel, nous calculons les intégrales à l'aide du logiciel [PySCF](https://pyscf.org/) logiciel.\n",
        "\n",
        "Pour plus de détails sur la manière dont le hamiltonien moléculaire est dérivé, consultez un manuel de chimie quantique (par exemple, *Modern Quantum Chemistry* de Szabo et Ostlund). Pour une explication de haut niveau de la manière dont les problèmes de chimie quantique sont transposés sur les ordinateurs quantiques, consultez la conférence [*Mapping Problems to Qubits*](https://youtube.com/watch?v=TyFU6r8uEsE\\&t=900) de l'université d'été mondiale Qiskit 2024.\n",
        "\n",
        "<span id=\"local-unitary-cluster-jastrow-lucj-ansatz\" />\n",
        "\n",
        "### Approche du cluster unitaire local de Jastrow (LUCJ)\n",
        "\n",
        "La méthode SQD nécessite un modèle de circuit quantique à partir duquel prélever des échantillons. Dans ce tutoriel, nous utiliserons l'approche [LUCJ (Local Unitary Cluster Jastrow)](https://pubs.rsc.org/en/content/articlelanding/2023/sc/d3sc02516k) en raison de sa combinaison entre justification physique et facilité de mise en œuvre matérielle. Nous utiliserons [ffsim](https://qiskit-community.github.io/ffsim/) pour construire le circuit de référence.\n",
        "\n",
        "L'approche LUCJ s'adapte aux QPU dont la connectivité des qubits est limitée. Les orbitales de spin sont mappées sur des qubits de telle sorte que l'ansatz ne nécessite pas de routage à l'aide de portes SWAP. IBM® Le matériel présente une topologie de qubits en réseau hexagonal dense; dans ce cas, nous pouvons adopter un motif en « zigzag », illustré ci-dessous. Dans ce schéma, les orbitales de même spin sont associées à des qubits selon une topologie en ligne (cercles rouges et bleus), et une connexion entre des orbitales de spin différent est présente tous les quatre orbitales spatiales, cette connexion étant assurée par un qubit auxiliaire (cercles violets).\n",
        "\n",
        "![Diagramme de cartographie du Qubit pour l'ansatz LUCJ sur un réseau lourd-hex](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",
        "### Récupération de configuration auto-cohérente\n",
        "\n",
        "La procédure de récupération de la configuration autoconsistante est conçue pour extraire autant de signaux que possible d'échantillons quantiques bruyants. Comme l'hamiltonien moléculaire conserve le nombre de particules et le spin Z, il est logique de choisir un ansatz de circuit qui conserve également ces symétries. Lorsqu'il est appliqué à l'état de Hartree-Fock, l'état résultant a un nombre de particules et un spin Z fixes dans un environnement sans bruit. Par conséquent, les moitiés de spin $\\alpha$ et de spin $\\beta$ de toute chaîne de bits échantillonnée à partir de cet état devraient avoir le même [poids de Hamming](https://en.wikipedia.org/wiki/Hamming_weight) que dans l'état Hartree-Fock. En raison de la présence de bruit dans les processeurs quantiques actuels, certaines chaînes de bits mesurées ne respectent pas cette propriété. Une forme simple de postsélection permettrait d'écarter ces chaînes de bits, mais c'est un gaspillage car ces chaînes de bits peuvent encore contenir un certain signal. La procédure de récupération autoconsistante tente de récupérer une partie de ce signal lors du post-traitement. La procédure est itérative et nécessite en entrée une estimation de l'occupation moyenne de chaque orbitale dans l'état fondamental, qui est d'abord calculée à partir des échantillons bruts. La procédure est exécutée en boucle et chaque itération comporte les étapes suivantes :\n",
        "\n",
        "1. Pour chaque chaîne de bits qui ne respecte pas les symétries spécifiées, les bits sont retournés selon une procédure probabiliste conçue pour rapprocher la chaîne de bits de l'estimation actuelle des occupations orbitales moyennes, afin d'obtenir une nouvelle chaîne de bits.\n",
        "2. Rassembler toutes les chaînes de bits anciennes et nouvelles qui satisfont aux symétries et sous-échantillonner des sous-ensembles d'une taille fixe, choisie à l'avance.\n",
        "3. Pour chaque sous-ensemble de chaînes de bits, projetez l'hamiltonien dans le sous-espace couvert par les vecteurs de base correspondants (voir la [section précédente](#quantum-chemistry) pour une description de ces vecteurs de base) et calculez une estimation de l'état fondamental de l'hamiltonien projeté sur un ordinateur classique.\n",
        "4. Mettre à jour l'estimation de l'occupation moyenne des orbitales avec l'estimation de l'état fondamental ayant l'énergie la plus basse.\n",
        "\n",
        "<span id=\"sqd-workflow-diagram\" />\n",
        "\n",
        "### Diagramme du flux de travail SQD\n",
        "\n",
        "Le flux de travail de la SQD est décrit dans le diagramme suivant :\n",
        "\n",
        "![Diagramme de flux de travail de l'algorithme SQD](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",
        "## Exigences\n",
        "\n",
        "Avant de commencer ce tutoriel, assurez-vous que les éléments suivants sont installés :\n",
        "\n",
        "* Qiskit SDK v1.0 ou plus tard, avec prise en charge de [la visualisation](/docs/api/qiskit/visualization)\n",
        "* Qiskit Runtime v0.22 ou plus tard (`pip install qiskit-ibm-runtime`)\n",
        "* Module complémentaire SQD Qiskit v0.11 ou version ultérieure (`pip install qiskit-addon-sqd`)\n",
        "* ffsim v0.0.75 ou version ultérieure (`pip install ffsim`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c6e44a31",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## Configuration\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",
        "## Exemple de simulateur à petite échelle\n",
        "\n",
        "Dans ce tutoriel, nous allons déterminer une approximation de l'état fondamental d'une molécule d'azote à une distance de liaison proche de son équilibre. Nous utilisons d'abord un petit ensemble de bases d' STO-6G s afin de pouvoir simuler l'expérience et vérifier qu'elle fonctionne.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "afeb054c",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique\n",
        "\n",
        "Tout d'abord, nous définissons la molécule et ses propriétés.\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)], [\"N\", (1.0, 0, 0)]],\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": [
        "Avant de construire le circuit de l'ansatz LUCJ, nous effectuons d'abord un calcul CCSD dans la cellule de code suivante. Les [amplitudes $t_1$ et $t_2$](https://en.wikipedia.org/wiki/Coupled_cluster#Cluster_operator) de ce calcul seront utilisées pour initialiser les paramètres de l'ansatz.\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.1283731437052354\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": [
        "Nous utilisons maintenant [ffsim](https://github.com/qiskit-community/ffsim) pour créer le circuit de référence. Comme notre molécule présente un état de Hartree-Fock à couche fermée, nous utilisons la variante à spin équilibré de l'approche UCJ, [UCJOpSpinBalanced](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced). Nous avons défini `optimize=True` dans la `from_t_amplitudes` méthode afin d'activer la double factorisation « compressée » des amplitudes de l' $t_2$ e (pour plus de détails, voir [la section « The local unitary cluster Jastrow (LUCJ) ansatz](https://qiskit-community.github.io/ffsim/explanations/lucj.html#Parameter-initialization-from-CCSD) » dans la documentation de ffsim).\n",
        "\n",
        "Étant donné que l'ansatz LUCJ s'adapte à la connectivité disponible du QPU, nous devons initialiser le backend du QPU avant de créer l'ansatz. Pour l'instant, nous allons créer un backend générique doté d'une carte à couplage hexagonal fort et d'un ensemble de portes auquel l'approche LUCJ se décompose naturellement. Ensuite, nous utiliserons `ffsim.qiskit.generate_lucj_pass_manager` pour créer un gestionnaire de passes spécialisé dans la transcompilation de l'approche LUCJ vers le backend spécifié, conformément à la structure en « zigzag » décrite dans la [section consacrée au contexte de l'approche LUCJ](#local-unitary-cluster-jastrow-lucj-ansatz). Cette fonction utilise une heuristique de notation pour minimiser les erreurs associées à la configuration sélectionnée, ce qui est important si votre backend est un véritable QPU ou un simulateur doté d'un modèle de bruit. Outre le gestionnaire de passes, cette fonction renvoie également les paires de couplage alpha-bêta pouvant être mises en œuvre sur le matériel. Si toutes les paires ne peuvent pas être mises en œuvre, un avertissement s'affiche.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "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 = 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",
        "coupling_map = CouplingMap.from_heavy_hex(3)\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=\"heavy-hex\",\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",
        "    # Limit the number of optimization iterations to prevent the code cell\n",
        "    # from running too long. Removing this line may improve results.\n",
        "    options=dict(maxiter=1000),\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",
        "### Étape 2 : Optimisation pour l'exécution sur du matériel quantique\n",
        "\n",
        "Nous optimisons ensuite le circuit pour le matériel cible. En général, cette étape consiste à initialiser le backend matériel et un gestionnaire de passes pour ce backend. Cependant, comme l'approche LUCJ est adaptée à la connectivité matérielle, nous avons déjà effectué ces opérations à l'étape précédente. Il ne reste plus qu'à exécuter le gestionnaire de passes sur le circuit pour le transcompiler en un circuit ISA pouvant être exécuté directement sur le QPU.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "7d554aa5",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Gate counts: OrderedDict({'xx_plus_yy': 86, 'p': 16, 'measure': 16, 'cp': 15, 'x': 10, 'swap': 2, '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",
        "### Étape 3 : Exécutez à l'aide d' Qiskit primitives\n",
        "\n",
        "Après avoir optimisé le circuit en vue de son exécution sur le matériel, nous sommes prêts à l'exécuter sur le matériel cible et à collecter des échantillons afin d'estimer l'énergie de l'état fondamental. Comme nous ne disposons que d'un seul circuit, nous allons utiliser [le mode d'exécution « Job »](/docs/guides/execution-modes) du service de calcul d' IBM Quantum pour exécuter notre circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "93c1cef3-298e-4deb-8512-769fe94cd5a5",
      "metadata": {},
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Warning: Trying to add QuantumRegister to a QuantumCircuit having a layout\n"
          ]
        }
      ],
      "source": [
        "rng = np.random.default_rng()\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",
        "### Étape 4 : Post-traitement et restitution du résultat dans le format classique souhaité\n",
        "\n",
        "Un indicateur utile pour évaluer la qualité du résultat du QPU est le nombre de configurations valides renvoyées. Une configuration valide possède le nombre correct de particules et le spin Z correct, ce qui signifie que la moitié droite de la chaîne binaire a un poids de Hamming égal au nombre d'électrons à spin up, et que la moitié gauche a un poids de Hamming égal au nombre d'électrons à spin down. La cellule suivante calcule la fraction des configurations échantillonnées qui sont valides.\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": [
        "Toutes les chaînes de bits sont valides, car nous effectuons un échantillonnage du circuit sur un simulateur sans bruit. Lorsqu'on exécute le programme sur un QPU bruyant, cette fraction sera inférieure à un, mais on espère qu'elle sera supérieure à celle à laquelle on s'attendrait si les chaînes de bits avaient été échantillonnées de manière aléatoire et uniforme, ce qui est calculé dans la cellule suivante.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "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": [
        "Nous estimons maintenant l'énergie de l'état fondamental de l'hamiltonien à l'aide de la fonction `diagonalize_fermionic_hamiltonian` . Cette fonction exécute la procédure de récupération de la configuration autoconsistante pour affiner de manière itérative les échantillons quantiques bruyants afin d'améliorer l'estimation de l'énergie. Nous passons une fonction de rappel afin de pouvoir enregistrer les résultats intermédiaires pour une analyse ultérieure. Voir la [documentation de l'API](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) pour des explications sur les arguments de `diagonalize_fermionic_hamiltonian`.\n",
        "\n",
        "Ici, nous utilisons `initial_occupancies` l'argument pour `diagonalize_fermionic_hamiltonian` spécifier la configuration Hartree-Fock comme estimation initiale pour les occupations orbitales dans l'état fondamental. Cette approche est judicieuse pour les systèmes dont l'état fondamental repose largement sur la configuration Hartree-Fock, mais elle peut ne pas convenir dans d'autres situations, même si des méthodes de calcul plus avancées peuvent permettre d'obtenir de meilleures estimations initiales dans ces cas-là. La spécification `initial_occupancies` permet également d'exécuter la récupération de configuration même si aucune configuration valide n'a été échantillonnée, comme cela peut être le cas lors de l'échantillonnage d'un grand circuit sur un QPU bruité. Sans cet argument, la récupération de la configuration échouerait et générerait une erreur si aucune configuration valide n'était fournie.\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.59275573641656\n",
            "\t\tSubspace dimension: 900\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -108.59275573641656\n",
            "\t\tSubspace dimension: 900\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -108.59275573641656\n",
            "\t\tSubspace dimension: 900\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -108.59275573641656\n",
            "\t\tSubspace dimension: 900\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -108.59275573641656\n",
            "\t\tSubspace dimension: 900\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -108.59275573641656\n",
            "\t\tSubspace dimension: 900\n",
            "Final energy: -108.59275573641656\n",
            "Final energy error: 0.0032316145694579745\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 = 300\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",
        "    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",
        "#### Visualisez les résultats\n",
        "\n",
        "Le premier graphique montre que, dans cette simulation, nous sommes déjà très `1 mH` proches de la réponse exacte après la première itération (on considère généralement que la précision chimique est de l'ordre de `1 kcal/mol`$\\approx$`1.6 mH`). Il s'agit toutefois d'un petit système et, comme les échantillons sont exempts de bruit, il n'est pas nécessaire de procéder à une reconstruction de la configuration. Sur un système plus important fonctionnant avec un QPU bruyant, plusieurs itérations de récupération de la configuration peuvent s'avérer nécessaires, et la précision finale peut s'en trouver réduite. En général, on peut améliorer l'énergie en autorisant davantage d'itérations de récupération de configuration ou en augmentant le nombre d'échantillons par lot.\n",
        "\n",
        "Le deuxième graphique montre l'occupation moyenne de chaque orbite spatiale après la dernière itération. Nous pouvons constater que les électrons de spin-up et de spin-down occupent les cinq premières orbitales avec une forte probabilité dans nos solutions.\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",
        "yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]\n",
        "\n",
        "# Chemical accuracy (+/- 1 milli-Hartree)\n",
        "chem_accuracy = 0.001\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_yticks(yt1)\n",
        "axs[0].set_yticklabels(yt1)\n",
        "axs[0].set_yscale(\"log\")\n",
        "axs[0].set_ylim(1e-4)\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",
        "## Exemple de matériel à grande échelle\n",
        "\n",
        "Nous allons maintenant exécuter un exemple plus complexe sur du matériel quantique réel. Nous allons ici dériver un espace actif pour la molécule d'azote à partir de la base de données « cc-pVDZ ».\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "24ca3090-3f3b-4efb-a482-70b2e1b5d062",
      "metadata": {},
      "source": [
        "<span id=\"steps-1-4\" />\n",
        "\n",
        "### Étapes 1 à 4\n",
        "\n",
        "Nous regroupons ici toutes ces étapes au sein d'un processus unique à plus grande échelle, qui est ensuite exécuté sur du matériel quantique réel.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "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.2879500329450045\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': 7039, 'rz': 6990, 'cz': 1858, 'x': 61, 'measure': 52, 'barrier': 1})\n",
            "Fraction of sampled configurations that are valid: 0.02124\n",
            "Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07\n",
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.13889134249762\n",
            "\t\tSubspace dimension: 120409\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.11785470455858\n",
            "\t\tSubspace dimension: 110889\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.13234360554011\n",
            "\t\tSubspace dimension: 130321\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.16392179579177\n",
            "\t\tSubspace dimension: 223729\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.16281938332986\n",
            "\t\tSubspace dimension: 223729\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.16955816711932\n",
            "\t\tSubspace dimension: 233289\n",
            "Iteration 3\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.17905772999075\n",
            "\t\tSubspace dimension: 324900\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.17532445048462\n",
            "\t\tSubspace dimension: 357604\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.1733168689756\n",
            "\t\tSubspace dimension: 348100\n",
            "Iteration 4\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.18437778820451\n",
            "\t\tSubspace dimension: 474721\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.18450164209159\n",
            "\t\tSubspace dimension: 476100\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.18493571190754\n",
            "\t\tSubspace dimension: 487204\n",
            "Iteration 5\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -109.18616522497996\n",
            "\t\tSubspace dimension: 622521\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -109.18652868888333\n",
            "\t\tSubspace dimension: 644809\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -109.18753326484406\n",
            "\t\tSubspace dimension: 585225\n",
            "Final energy: -109.18753326484406\n",
            "Final energy error: 0.040495951813099396\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)], [\"N\", (1.0, 0, 0)]],\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(\n",
        "    operational=True, simulator=False, min_num_qubits=133\n",
        ")\n",
        "print(f\"Using backend {backend.name}\")\n",
        "\n",
        "# Create pass manager\n",
        "pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(\n",
        "    backend=backend,\n",
        "    norb=norb,\n",
        "    connectivity=\"heavy-hex\",\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",
        "    # Limit the number of optimization iterations to prevent the code cell\n",
        "    # from running too long. Removing this line may improve results.\n",
        "    options=dict(maxiter=1000),\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 = 300\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",
        "    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",
        "yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]\n",
        "\n",
        "# Chemical accuracy (+/- 1 milli-Hartree)\n",
        "chem_accuracy = 0.001\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_yticks(yt1)\n",
        "axs[0].set_yticklabels(yt1)\n",
        "axs[0].set_yscale(\"log\")\n",
        "axs[0].set_ylim(1e-4)\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",
        "## Etapes suivantes\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Recommandations\">\n",
        "  Si ce travail vous a intéressé, les documents suivants pourraient vous intéresser :\n",
        "\n",
        "  * [Diagonalisation quantique de Krylov par échantillonnage d'un modèle de réseau fermionique](/docs/tutorials/sample-based-krylov-quantum-diagonalization) – tutoriel associé utilisant des circuits d'évolution temporelle à la place d'une approche variationnelle\n",
        "  * [Optimiser les flux de travail chimiques SQD grâce au solveur Dice](/docs/addons/qiskit-addon-sqd/guides/integrate-dice-solver) – une page expliquant comment utiliser le logiciel Dice, plus performant, pour la diagonalisation\n",
        "  * [Documentation de l'API de l'extension SQD](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) - référence pour la `diagonalize_fermionic_hamiltonian` fonction\n",
        "  * [*La chimie au-delà de l'échelle de la diagonalisation exacte sur un supercalculateur quantique*](https://www.science.org/doi/10.1126/sciadv.adu9991) – l'article sur lequel s'appuie ce tutoriel\n",
        "</Admonition>\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"
    },
    "hours": 1.5,
    "qpuSeconds": 60
  },
  "nbformat": 4,
  "nbformat_minor": 5
}