{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "de193554-b271-4295-95e4-8904f0f6ee8a",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Géométrie moléculaire\"\n",
        "description: \"Dans cette leçon, nous modifions la géométrie d'une molécule simple, en minimisant l'énergie à chaque étape.\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore pxxr prqs nelecas mcscf chmax Dmax vmax ecore ncas Excp disp */}\n",
        "\n",
        "<span id=\"determining-a-molecular-geometry\" />\n",
        "\n",
        "# Détermination d'une géométrie moléculaire\n",
        "\n",
        "Dans la section précédente, nous avons mis en œuvre VQE pour déterminer l'énergie de l'état fondamental d'une molécule. Il s'agit là d'une utilisation valable de l'informatique quantique, mais il serait encore plus utile de déterminer la structure d'une molécule.\n",
        "\n",
        "<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",
        "En restant dans notre exemple de base de l'hydrogène diatomique, le seul paramètre géométrique à faire varier est la longueur de la liaison. Pour ce faire, nous procédons comme précédemment, mais en utilisant une variable dans notre construction initiale de la molécule (une longueur de liaison, *x*, dans l'argument). Il s'agit d'un changement assez simple, mais qui exige que la variable soit incluse dans les fonctions tout au long du processus, puisqu'elle commence dans la construction du hamiltonien fermionique et se propage à travers la cartographie et enfin dans la fonction de coût.\n",
        "\n",
        "Tout d'abord, nous chargeons certains des paquets que nous avons utilisés précédemment et nous définissons la fonction Cholesky.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "id": "0a5d39bd-8c26-404f-8057-c29e3af70df4",
      "metadata": {},
      "outputs": [],
      "source": [
        "from qiskit.quantum_info import SparsePauliOp\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "\n",
        "#!pip install pyscf==2.4.0\n",
        "from pyscf import ao2mo, gto, mcscf, scf\n",
        "\n",
        "\n",
        "def cholesky(V, eps):\n",
        "    # see https://arxiv.org/pdf/1711.02242.pdf section B2\n",
        "    # see https://arxiv.org/abs/1808.02625\n",
        "    # see https://arxiv.org/abs/2104.08957\n",
        "    no = V.shape[0]\n",
        "    chmax, ng = 20 * no, 0\n",
        "    W = V.reshape(no**2, no**2)\n",
        "    L = np.zeros((no**2, chmax))\n",
        "    Dmax = np.diagonal(W).copy()\n",
        "    nu_max = np.argmax(Dmax)\n",
        "    vmax = Dmax[nu_max]\n",
        "    while vmax > eps:\n",
        "        L[:, ng] = W[:, nu_max]\n",
        "        if ng > 0:\n",
        "            L[:, ng] -= np.dot(L[:, 0:ng], (L.T)[0:ng, nu_max])\n",
        "        L[:, ng] /= np.sqrt(vmax)\n",
        "        Dmax[: no**2] -= L[: no**2, ng] ** 2\n",
        "        ng += 1\n",
        "        nu_max = np.argmax(Dmax)\n",
        "        vmax = Dmax[nu_max]\n",
        "    L = L[:, :ng].reshape((no, no, ng))\n",
        "    print(\n",
        "        \"accuracy of Cholesky decomposition \",\n",
        "        np.abs(np.einsum(\"prg,qsg->prqs\", L, L) - V).max(),\n",
        "    )\n",
        "    return L, ng\n",
        "\n",
        "\n",
        "def identity(n):\n",
        "    return SparsePauliOp.from_list([(\"I\" * n, 1)])\n",
        "\n",
        "\n",
        "def creators_destructors(n, mapping=\"jordan_wigner\"):\n",
        "    c_list = []\n",
        "    if mapping == \"jordan_wigner\":\n",
        "        for p in range(n):\n",
        "            if p == 0:\n",
        "                ell, r = \"I\" * (n - 1), \"\"\n",
        "            elif p == n - 1:\n",
        "                ell, r = \"\", \"Z\" * (n - 1)\n",
        "            else:\n",
        "                ell, r = \"I\" * (n - p - 1), \"Z\" * p\n",
        "            cp = SparsePauliOp.from_list([(ell + \"X\" + r, 0.5), (ell + \"Y\" + r, -0.5j)])\n",
        "            c_list.append(cp)\n",
        "    else:\n",
        "        raise ValueError(\"Unsupported mapping.\")\n",
        "    d_list = [cp.adjoint() for cp in c_list]\n",
        "    return c_list, d_list"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "40e25b54-7927-4dbb-a26f-1c6b33f7f349",
      "metadata": {},
      "source": [
        "Pour définir notre hamiltonien, nous utiliserons PySCF exactement comme dans l'exemple précédent, mais nous inclurons maintenant une variable, `x`, qui jouera le rôle de notre distance interatomique. Cette opération renvoie l'énergie du noyau, l'énergie d'un seul électron et l'énergie de deux électrons comme précédemment.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "dbd10d0c-feb1-4a86-9bd9-b61101a08b95",
      "metadata": {},
      "outputs": [],
      "source": [
        "def ham_terms(x: float):\n",
        "    distance = x\n",
        "    a = distance / 2\n",
        "    mol = gto.Mole()\n",
        "    mol.build(\n",
        "        verbose=0,\n",
        "        atom=[\n",
        "            [\"H\", (0, 0, -a)],\n",
        "            [\"H\", (0, 0, a)],\n",
        "        ],\n",
        "        basis=\"sto-6g\",\n",
        "        spin=0,\n",
        "        charge=0,\n",
        "        symmetry=\"Dooh\",\n",
        "    )\n",
        "\n",
        "    # mf = scf.RHF(mol)\n",
        "    # mx = mcscf.CASCI(mf, ncas=2, nelecas=(1, 1))\n",
        "    # mx.kernel()\n",
        "\n",
        "    mf = scf.RHF(mol)\n",
        "    mf.kernel()\n",
        "    if not mf.converged:\n",
        "        raise RuntimeError(f\"SCF did not converge for distance {x}\")\n",
        "\n",
        "    mx = mcscf.CASCI(mf, ncas=2, nelecas=(1, 1))\n",
        "    casci_energy = mx.kernel()\n",
        "    if casci_energy is None:\n",
        "        raise RuntimeError(f\"CASCI failed for distance {x}\")\n",
        "\n",
        "    # Other variables that might come in handy:\n",
        "    # active_space = range(mol.nelectron // 2 - 1, mol.nelectron // 2 + 1)\n",
        "    #    E1 = mf.kernel()\n",
        "    # mo = mx.sort_mo(active_space, base=0)\n",
        "    #    E2 = mx.kernel(mo)[:2]\n",
        "\n",
        "    h1e, ecore = mx.get_h1eff()\n",
        "    h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)\n",
        "    return ecore, h1e, h2e"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "db49a702-e60e-4c9b-8ff2-94aa2ada022c",
      "metadata": {},
      "source": [
        "Rappelons que la construction ci-dessus consiste à créer un hamiltonien fermionique basé sur l'espèce atomique, la géométrie et les orbitales électroniques. Ci-dessous, nous mettons en correspondance cet hamiltonien fermionique avec les opérateurs de Pauli. Cette fonction `build_hamiltonian` comprendra également une variable géométrique comme argument.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "id": "84e6a56b-eaea-4c0a-8502-567cfc5140a2",
      "metadata": {},
      "outputs": [],
      "source": [
        "def build_hamiltonian(distx: float) -> SparsePauliOp:\n",
        "    ecore = ham_terms(distx)[0]\n",
        "    h1e = ham_terms(distx)[1]\n",
        "    h2e = ham_terms(distx)[2]\n",
        "\n",
        "    ncas, _ = h1e.shape\n",
        "\n",
        "    C, D = creators_destructors(2 * ncas, mapping=\"jordan_wigner\")\n",
        "    Exc = []\n",
        "    for p in range(ncas):\n",
        "        Excp = [C[p] @ D[p] + C[ncas + p] @ D[ncas + p]]\n",
        "        for r in range(p + 1, ncas):\n",
        "            Excp.append(\n",
        "                C[p] @ D[r]\n",
        "                + C[ncas + p] @ D[ncas + r]\n",
        "                + C[r] @ D[p]\n",
        "                + C[ncas + r] @ D[ncas + p]\n",
        "            )\n",
        "        Exc.append(Excp)\n",
        "\n",
        "    # low-rank decomposition of the Hamiltonian\n",
        "    Lop, ng = cholesky(h2e, 1e-6)\n",
        "    t1e = h1e - 0.5 * np.einsum(\"pxxr->pr\", h2e)\n",
        "\n",
        "    H = ecore * identity(2 * ncas)\n",
        "    # one-body term\n",
        "    for p in range(ncas):\n",
        "        for r in range(p, ncas):\n",
        "            H += t1e[p, r] * Exc[p][r - p]\n",
        "    # two-body term\n",
        "    for g in range(ng):\n",
        "        Lg = 0 * identity(2 * ncas)\n",
        "        for p in range(ncas):\n",
        "            for r in range(p, ncas):\n",
        "                Lg += Lop[p, r, g] * Exc[p][r - p]\n",
        "        H += 0.5 * Lg @ Lg\n",
        "\n",
        "    return H.chop().simplify()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ffcb2569-f832-4fa9-a8bb-d7626c2a233d",
      "metadata": {},
      "source": [
        "Nous chargerons les paquets restants pour exécuter VQE lui-même, tels que l'ansatz efficient\\_su2 et les minimiseurs SciPy :\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "5f8dd8e0-6dfb-4be7-af15-1591f569201c",
      "metadata": {},
      "outputs": [],
      "source": [
        "# General imports\n",
        "\n",
        "# Pre-defined ansatz circuit and operator class for Hamiltonian\n",
        "from qiskit.circuit.library import efficient_su2\n",
        "from qiskit.quantum_info import SparsePauliOp\n",
        "\n",
        "# SciPy minimizer routine\n",
        "from scipy.optimize import minimize\n",
        "\n",
        "# Plotting functions\n",
        "\n",
        "# qiskit-ibm-runtime tools\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "\n",
        "service = QiskitRuntimeService()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f94c1680-9c2a-4585-a63e-a49da4eb02f3",
      "metadata": {},
      "source": [
        "Nous allons à nouveau définir la fonction de coût, mais celle-ci prend toujours comme argument un hamiltonien entièrement construit et cartographié, de sorte que rien ne change à propos de cette fonction.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "id": "0a80ad8c-a9cb-4cab-835d-f65131b99c87",
      "metadata": {},
      "outputs": [],
      "source": [
        "def cost_func(params, ansatz, H, estimator):\n",
        "    pub = (ansatz, [H], [params])\n",
        "    result = estimator.run(pubs=[pub]).result()\n",
        "    energy = result[0].data.evs[0]\n",
        "    return energy\n",
        "\n",
        "\n",
        "# def cost_func_sim(params, ansatz, H, estimator):\n",
        "#    energy = estimator.run(ansatz, H, parameter_values=params).result().values[0]\n",
        "#    return energy"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a068e645-8dfb-4e9c-b1b1-7dd936808188",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-problem-for-quantum-execution\" />\n",
        "\n",
        "## Étape 2 : Optimiser le problème pour l'exécution quantique\n",
        "\n",
        "Comme l'hamiltonien change à chaque nouvelle géométrie, le transpiling de l'opérateur change à chaque étape. Nous pouvons néanmoins définir un gestionnaire de passe général à appliquer à chaque étape, spécifique au matériel que nous voulons utiliser.\n",
        "\n",
        "Nous utiliserons ici le backend le moins chargé disponible. Nous utiliserons ce backend comme modèle pour notre AerSimulator,, ce qui permettra à notre simulateur d'imiter, par exemple, le comportement du backend réel en matière de bruit. Ces modèles de bruit ne sont pas parfaits, mais ils peuvent vous aider à savoir à quoi vous attendre avec du matériel réel.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "fa07518f-a2c0-4f7a-b344-8d9576427478",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Here, we select the least busy backend available:\n",
        "backend = service.least_busy(operational=True, simulator=False)\n",
        "print(backend)\n",
        "# Or to select a specific real backend use the line below, and substitute 'ibm_strasbourg'\n",
        "# for your chosen device. backend = service.get_backend('ibm_strasbourg')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "id": "e67fc84a-431b-4efd-937f-49c8b7ac3abb",
      "metadata": {},
      "outputs": [],
      "source": [
        "# To run on a simulator:\n",
        "# -----------\n",
        "from qiskit_aer import AerSimulator\n",
        "\n",
        "backend_sim = AerSimulator.from_backend(backend)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0636c55f-f46f-46ca-acaf-72bcc2f5f663",
      "metadata": {},
      "source": [
        "Nous importons le gestionnaire de passes et les paquets connexes pour nous aider à optimiser notre circuit. Cette étape, et celle qui la précède, sont indépendantes de l'hamiltonien et sont donc inchangées par rapport à la leçon précédente.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "id": "8202332e-69ca-4049-af27-e98a77f15a5d",
      "metadata": {},
      "outputs": [],
      "source": [
        "from qiskit.transpiler import PassManager\n",
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "from qiskit.transpiler.passes import (\n",
        "    ALAPScheduleAnalysis,\n",
        "    PadDynamicalDecoupling,\n",
        "    ConstrainedReschedule,\n",
        ")\n",
        "from qiskit.circuit.library import XGate\n",
        "\n",
        "target = backend.target\n",
        "pm = generate_preset_pass_manager(target=target, optimization_level=3)\n",
        "pm.scheduling = PassManager(\n",
        "    [\n",
        "        ALAPScheduleAnalysis(target=target),\n",
        "        ConstrainedReschedule(\n",
        "            acquire_alignment=target.acquire_alignment,\n",
        "            pulse_alignment=target.pulse_alignment,\n",
        "            target=target,\n",
        "        ),\n",
        "        PadDynamicalDecoupling(\n",
        "            target=target,\n",
        "            dd_sequence=[XGate(), XGate()],\n",
        "            pulse_alignment=target.pulse_alignment,\n",
        "        ),\n",
        "    ]\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4932ef5e-3c81-46a1-92cb-3082399734a0",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-ibm-quantum-primitives\" />\n",
        "\n",
        "## Étape 3 : Exécuter à l'aide des primitives de l' IBM Quantum.\n",
        "\n",
        "Dans le bloc de code ci-dessous, nous définissons un tableau destiné à stocker les résultats obtenus à chaque étape de notre calcul de distance interatomique $x$. Nous avons choisi l'intervalle de valeurs $x$ en nous basant sur la valeur expérimentale de la longueur de liaison à l'équilibre : 0.74 angströms. Nous allons d'abord l'exécuter sur un simulateur, et nous allons donc importer notre Estimator ( BackendEstimator ) depuis `qiskit.primitives`. Pour chaque étape géométrique, nous construisons l'hamiltonien et effectuons un certain nombre d'itérations d'optimisation (ici 500) à l'aide de l'optimiseur « cobyla ». À chaque étape de la modélisation géométrique, nous enregistrons à la fois l'énergie totale et l'énergie électronique. En raison du nombre élevé d'étapes de l'optimiseur, cette opération peut prendre une heure, voire plus. Vous pouvez modifier les données ci-dessous afin de réduire le temps nécessaire.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "id": "4c41a221-8a02-4932-882b-5afcc98d1d8d",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "accuracy of Cholesky decomposition  1.1102230246251565e-15\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/home/porter284/.pyenv/versions/3.11.12/lib/python3.11/site-packages/scipy/_lib/pyprima/common/preproc.py:68: UserWarning: COBYLA: Invalid MAXFUN; it should be at least num_vars + 2; it is set to 34\n",
            "  warn(f'{solver}: Invalid MAXFUN; it should be at least {min_maxfun_str}; it is set to {maxfun}')\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = 1.316011435623847\n",
            "The corresponding X is:\n",
            "[2.32948769 5.39918229 3.03787975 4.11789904 4.97130735 2.68662232\n",
            " 1.76573151 2.48982571 5.40431972 3.65780829 1.33792786 5.48472494\n",
            " 6.18738702 1.78741883 0.78195251 2.96658955 1.35827677 5.599321\n",
            " 4.54850148 1.0939048  4.26158726 0.52100721 0.82318    4.76796961\n",
            " 3.75795507 3.8526447  5.51100375 5.91023075 2.61494836 1.79908918\n",
            " 2.65937756 5.53964148]\n",
            "\n",
            "-0.44791260077615314\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: 1.316011435623847\n",
            "       x: [ 2.329e+00  5.399e+00 ...  2.659e+00  5.540e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  5.551115123125783e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = 0.7235003672327549\n",
            "The corresponding X is:\n",
            "[2.56282915 5.63369524 5.58059887 4.049643   4.2021266  3.06866011\n",
            " 6.01619635 1.52520776 4.35403161 0.33673958 0.32623161 1.2179545\n",
            " 2.84001371 3.98956684 4.89632562 1.38303588 1.96194695 2.13182089\n",
            " 0.29739166 1.77895165 3.29151585 3.54355374 4.49626674 0.95756626\n",
            " 0.87103927 4.53068385 1.31051302 0.37103108 1.02961355 3.13342311\n",
            " 5.65815319 2.24770604]\n",
            "\n",
            "-0.5994426600672451\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: 0.7235003672327549\n",
            "       x: [ 2.563e+00  5.634e+00 ...  5.658e+00  2.248e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  5.551115123125783e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = 0.34960914928810116\n",
            "The corresponding X is:\n",
            "[5.44143165 6.75955835 1.56836472 3.09522093 4.67873235 1.67071481\n",
            " 0.3056494  0.65998337 1.02197668 5.21162959 0.43690354 3.56522934\n",
            " 4.56033119 1.90736037 0.40863891 2.87007312 3.2516952  5.90360196\n",
            " 1.99057799 5.20726456 0.74710237 6.03179202 3.80685028 0.03844391\n",
            " 5.88580196 3.62233258 3.98723567 2.50591888 5.44020267 2.2792993\n",
            " 5.57102303 4.46548617]\n",
            "\n",
            "-0.7087452725518989\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: 0.34960914928810116\n",
            "       x: [ 5.441e+00  6.760e+00 ...  5.571e+00  4.465e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  2.220446049250313e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = 0.10594558882184543\n",
            "The corresponding X is:\n",
            "[5.35675483 2.26629567 1.45430546 5.56758296 5.76309509 0.73239338\n",
            " 5.1216998  3.03258872 4.33624828 1.93197674 0.5292902  3.32274987\n",
            " 3.43247633 0.81490741 0.48060245 1.9944799  5.67519646 5.12534057\n",
            " 0.06510627 2.52989834 6.1699519  0.94828957 5.91634548 1.5994961\n",
            " 4.27902164 2.3129213  1.82353095 2.10634209 1.43740426 4.06988733\n",
            " 0.59624074 4.93925418]\n",
            "\n",
            "-0.7760164293781545\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: 0.10594558882184543\n",
            "       x: [ 5.357e+00  2.266e+00 ...  5.962e-01  4.939e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  1.1102230246251565e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.06473600797229297\n",
            "The corresponding X is:\n",
            "[6.07735568 0.18019501 0.20743128 4.15445985 3.59388894 5.10047555\n",
            " 6.09938474 6.54707528 3.36251167 2.05475223 3.67078456 5.96010605\n",
            " 2.58589996 5.2723619  3.26352977 2.47432334 3.50289983 2.06620525\n",
            " 6.0946056  1.22751903 0.97320057 2.19564095 5.73174941 2.05127682\n",
            " 5.73805165 3.84046105 1.84816963 2.1247504  3.11106736 2.44136052\n",
            " 3.39002685 0.81596991]\n",
            "\n",
            "-0.8207034521437214\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.06473600797229297\n",
            "       x: [ 6.077e+00  1.802e-01 ...  3.390e+00  8.160e-01]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  5.551115123125783e-17\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.19562982094782935\n",
            "The corresponding X is:\n",
            "[-0.02184462  3.67041038  7.25918653  5.89799546  0.63583624  1.84214506\n",
            "  2.84059837  5.31485182  1.6053784   0.04556618  0.32018993 -0.03884066\n",
            "  0.69131496  0.24203727  1.97397262  3.59723495  0.43355775  2.30131056\n",
            "  4.63482292  3.9857415   4.32320753  4.55388437  2.18753433  5.99034987\n",
            "  2.50489913  0.90650534  4.82518088  2.32954849  2.29901832  5.33658863\n",
            "  5.91246716  3.2405013 ]\n",
            "\n",
            "-0.8571013345978292\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.19562982094782935\n",
            "       x: [-2.184e-02  3.670e+00 ...  5.912e+00  3.241e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  1.1102230246251565e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.2833766309947055\n",
            "The corresponding X is:\n",
            "[ 3.1700088   5.05055456  1.2545611   4.28751811  0.6255103   1.67526577\n",
            "  5.48201473  4.83820497  7.34880059  5.99705431  4.2502643   0.32066274\n",
            "  0.41001404  0.27271241  4.15682546  4.22393693  4.35148115  0.64538137\n",
            "  5.26288622  5.03810489  4.62426621  4.74997689  1.09603919  0.34752466\n",
            "  1.8116275   0.7474807   5.31754143  4.11181763  1.58797998  5.6299796\n",
            "  3.0109383  -0.19062772]\n",
            "\n",
            "-0.8713513097947054\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.2833766309947055\n",
            "       x: [ 3.170e+00  5.051e+00 ...  3.011e+00 -1.906e-01]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  1.1102230246251565e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.3527503628484244\n",
            "The corresponding X is:\n",
            "[3.90513622 4.61398739 5.92552705 1.99953405 4.82157369 1.35702441\n",
            " 2.77701782 5.73612247 4.22710527 1.83463189 0.45796297 4.62509318\n",
            " 0.98998668 0.11666217 3.0234641  4.54298546 0.14034033 4.15635797\n",
            " 1.41257357 4.48719602 2.39365535 0.19672041 5.0763044  1.86357581\n",
            " 3.657757   4.60298344 2.49769577 1.88086199 3.00108725 1.84475841\n",
            " 5.24047385 4.91142914]\n",
            "\n",
            "-0.8819275737684243\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.3527503628484244\n",
            "       x: [ 3.905e+00  4.614e+00 ...  5.240e+00  4.911e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  2.7755575615628914e-17\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.4022181851996095\n",
            "The corresponding X is:\n",
            "[6.09453981 3.5109422  3.37216019 4.94732621 1.25662002 5.89645164\n",
            " 5.06403334 2.68073141 4.40385083 1.13638366 1.73347762 6.82932871\n",
            " 1.15265014 2.07145964 4.36520459 1.14960341 1.62288871 4.32315915\n",
            " 5.45622821 0.93554005 3.17418483 0.47230243 1.31535502 5.77698726\n",
            " 2.04927925 2.50663538 5.9706002  5.4984681  2.9421232  1.56636313\n",
            " 1.09394523 4.62582   ]\n",
            "\n",
            "-0.8832883769450639\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.4022181851996095\n",
            "       x: [ 6.095e+00  3.511e+00 ...  1.094e+00  4.626e+00]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "accuracy of Cholesky decomposition  1.1102230246251565e-16\n",
            "Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            "Number of function values = 34   Least value of F = -0.44423031870708934\n",
            "The corresponding X is:\n",
            "[4.05765050e+00 3.99144950e+00 3.13287593e+00 3.28855137e+00\n",
            " 4.32613515e+00 4.91104512e+00 1.86521867e+00 2.18822879e+00\n",
            " 6.01336171e+00 1.82501276e+00 2.64830637e+00 5.53045823e+00\n",
            " 2.36110093e+00 3.98821703e+00 4.69013438e-01 4.38996815e+00\n",
            " 7.78103801e-04 1.72994378e+00 2.24970934e+00 1.11978200e+00\n",
            " 2.24846445e+00 4.90745512e+00 5.38474921e+00 5.03587994e+00\n",
            " 3.54297277e+00 4.78147533e+00 1.25990218e+00 1.99168068e+00\n",
            " 5.89203503e+00 1.77673987e+00 5.37848357e+00 5.60245198e-01]\n",
            "\n",
            "-0.8852113278070892\n",
            " message: Return from COBYLA because the objective function has been evaluated MAXFUN times.\n",
            " success: False\n",
            "  status: 3\n",
            "     fun: -0.44423031870708934\n",
            "       x: [ 4.058e+00  3.991e+00 ...  5.378e+00  5.602e-01]\n",
            "    nfev: 34\n",
            "   maxcv: 0.0\n",
            "All energies have been calculated\n"
          ]
        }
      ],
      "source": [
        "from qiskit.primitives import BackendEstimatorV2\n",
        "\n",
        "estimator = BackendEstimatorV2(backend=backend_sim)\n",
        "\n",
        "distances_sim = np.arange(0.3, 1.3, 0.1)\n",
        "vqe_energies_sim = []\n",
        "vqe_elec_energies_sim = []\n",
        "\n",
        "for dist in distances_sim:\n",
        "    xx = dist\n",
        "\n",
        "    # Random initial state and efficient_su2 ansatz\n",
        "    H = build_hamiltonian(xx)\n",
        "    ansatz = efficient_su2(H.num_qubits)\n",
        "    ansatz_isa = pm.run(ansatz)\n",
        "    x0 = 2 * np.pi * np.random.random(ansatz_isa.num_parameters)\n",
        "    H_isa = H.apply_layout(ansatz_isa.layout)\n",
        "    nuclear_repulsion = ham_terms(xx)[0]\n",
        "\n",
        "    res = minimize(\n",
        "        cost_func,\n",
        "        x0,\n",
        "        args=(ansatz_isa, H_isa, estimator),\n",
        "        method=\"cobyla\",\n",
        "        options={\"maxiter\": 20, \"disp\": True},\n",
        "    )\n",
        "\n",
        "    # Note this returns the total energy, and we are often interested in the electronic energy\n",
        "    tot_energy = getattr(res, \"fun\")\n",
        "    electron_energy = getattr(res, \"fun\") - nuclear_repulsion\n",
        "    print(electron_energy)\n",
        "    vqe_energies_sim.append(tot_energy)\n",
        "    vqe_elec_energies_sim.append(electron_energy)\n",
        "\n",
        "    # Print all results\n",
        "    print(res)\n",
        "\n",
        "print(\"All energies have been calculated\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 39,
      "id": "bc6164ca-6909-4780-8009-6dc274c66268",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "np.float64(1.2000000000000004)"
            ]
          },
          "execution_count": 39,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "xx"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "633a2a89-950f-4d13-a56e-6adc079245ea",
      "metadata": {},
      "source": [
        "Les résultats de cette simulation sont abordés ci-dessous dans la section consacrée au post-traitement; pour l'instant, notons simplement que la simulation s'est déroulée avec succès. Vous êtes désormais prêt à exécuter le programme sur du matériel réel. Nous allons définir la résilience sur `1`, ce qui indique que le mécanisme d'atténuation des erreurs TREX sera utilisé. Maintenant que nous travaillons avec du matériel réel, nous allons utiliser le service de calcul IBM Quantum et les primitives IBM Quantum. Notez que tant la boucle « for » relative à la géométrie que les multiples essais variationnels se trouvent à l'intérieur de la session.\n",
        "\n",
        "En raison des coûts et des délais associés à l'utilisation de matériel réel, nous avons réduit le nombre d'étapes géométriques et d'étapes d'optimisation ci-dessous. Veillez à adapter ces étapes en fonction de vos objectifs de précision et des délais impartis.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b09744f5-87a9-4744-b495-cf993e5ffcb3",
      "metadata": {},
      "outputs": [],
      "source": [
        "# To continue running on real hardware use\n",
        "from qiskit_ibm_runtime import Session\n",
        "from qiskit_ibm_runtime import EstimatorV2 as Estimator\n",
        "from qiskit_ibm_runtime import EstimatorOptions\n",
        "\n",
        "estimator_options = EstimatorOptions(resilience_level=1, default_shots=2000)\n",
        "\n",
        "distances = np.arange(0.5, 0.9, 0.1)\n",
        "vqe_energies = []\n",
        "vqe_elec_energies = []\n",
        "\n",
        "with Session(backend=backend) as session:\n",
        "    estimator = Estimator(mode=session, options=estimator_options)\n",
        "\n",
        "    for dist in distances:\n",
        "        xx = dist\n",
        "\n",
        "        # Random initial state and efficient_su2 ansatz\n",
        "\n",
        "        H = build_hamiltonian(xx)\n",
        "        ansatz = efficient_su2(H.num_qubits)\n",
        "        ansatz_isa = pm.run(ansatz)\n",
        "        H_isa = H.apply_layout(ansatz_isa.layout)\n",
        "        nuclear_repulsion = ham_terms(xx)[0]\n",
        "        x0 = 2 * np.pi * np.random.random(ansatz_isa.num_parameters)\n",
        "\n",
        "        res = minimize(\n",
        "            cost_func,\n",
        "            x0,\n",
        "            args=(ansatz_isa, H_isa, estimator),\n",
        "            method=\"cobyla\",\n",
        "            options={\"maxiter\": 50, \"disp\": True},\n",
        "        )\n",
        "\n",
        "        # Note this returns the total energy, and we are often interested in the electronic energy\n",
        "        tot_energy = getattr(res, \"fun\")\n",
        "        electron_energy = getattr(res, \"fun\") - nuclear_repulsion\n",
        "        print(electron_energy)\n",
        "        vqe_energies.append(tot_energy)\n",
        "        vqe_elec_energies.append(electron_energy)\n",
        "\n",
        "        # Print all results\n",
        "        print(res)\n",
        "\n",
        "print(\"All energies have been calculated\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5b6f1491-88bd-4fb1-a1fe-2e29ffa33f17",
      "metadata": {},
      "source": [
        "<span id=\"step-4-post-processing\" />\n",
        "\n",
        "## Étape 4 : Post-traitement\n",
        "\n",
        "Tant pour le simulateur que pour le matériel réel, nous pouvons tracer les énergies de l'état fondamental calculées pour chaque distance inter-atomique et voir où l'énergie la plus basse est atteinte. C'est la distance interatomique que l'on trouve dans la nature, et elle est en effet très proche. Une courbe plus régulière pourrait être obtenue en essayant d'autres réponses, d'autres optimiseurs et en effectuant le calcul plusieurs fois à chaque étape de la géométrie et en calculant la moyenne sur plusieurs conditions initiales aléatoires.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e81f3ead-27ac-415e-a9f2-64a51d4b7aa3",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/quantum-chem-with-vqe/geometry/extracted-outputs/e81f3ead-27ac-415e-a9f2-64a51d4b7aa3-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Here we can plot the results from this simulation.\n",
        "plt.plot(distances_sim, vqe_energies_sim, label=\"VQE Energy\")\n",
        "plt.xlabel(\"Atomic distance (Angstrom)\")\n",
        "plt.ylabel(\"Energy\")\n",
        "plt.legend()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1efe0d76-3c6c-4369-9ace-fc890084e676",
      "metadata": {},
      "source": [
        "Il est à noter que le simple fait d'augmenter le nombre d'étapes d'optimisation n'est pas susceptible d'améliorer les résultats du simulateur, puisque toutes les optimisations ont convergé vers la tolérance requise en moins de temps que le nombre maximum d'itérations.\n",
        "\n",
        "Les résultats obtenus avec le matériel réel sont comparables, à l'exception d'une gamme de valeurs échantillonnées légèrement différente.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "de8f53cf-9547-4578-b6cb-20d2b5602ee0",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/quantum-chem-with-vqe/geometry/extracted-outputs/de8f53cf-9547-4578-b6cb-20d2b5602ee0-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.plot(distances, vqe_energies, label=\"VQE Energy\")\n",
        "plt.xlabel(\"Atomic distance (Angstrom)\")\n",
        "plt.ylabel(\"Energy\")\n",
        "plt.legend()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7e732366-425a-40e4-b6f2-d26f11f7d86b",
      "metadata": {},
      "source": [
        "En plus d'attendre une longueur de liaison H2 de 0.74 Angstrom, l'énergie totale devrait être de -1.17 Hartrees. Nous constatons que les résultats du matériel réel sont plus proches de ces valeurs que ceux du simulateur. Cela est probablement dû au fait que le bruit était présent (ou simulé) dans les deux cas, mais ce n'est que dans le cas du matériel réel que des mesures d'atténuation des erreurs ont été employées.\n",
        "\n",
        "<span id=\"closing\" />\n",
        "\n",
        "### Fermeture en cours\n",
        "\n",
        "Ceci conclut notre cours sur le VQE pour la chimie quantique. Si vous souhaitez comprendre la théorie de l'information sous-jacente utilisée dans l'informatique quantique, consultez le cours de John Watrous sur les [bases de l'information quantique](/learning/courses/basics-of-quantum-information). Pour un exemple plus court d'un flux de travail VQE, voir notre [tutoriel sur l'estimation de l'énergie de l'état fondamental de la chaîne de Heisenberg à l'aide de VQE](/docs/tutorials/spin-chain-vqe). Vous pouvez également consulter les [didacticiels](/docs/tutorials) et les [cours](/learning) pour trouver d'autres supports pédagogiques sur les dernières technologies en matière d'informatique quantique.\n",
        "\n",
        "N'oubliez pas de passer l'examen de ce cours. Un score de 80 % ou plus vous permettra d'obtenir un badge Credly, qui vous sera automatiquement envoyé par courrier électronique. Merci de faire partie du réseau IBM Quantum®!\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0ac6d249-d210-479d-a792-c8b4e94b8b88",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "1.3.2\n",
            "0.35.0\n"
          ]
        }
      ],
      "source": [
        "import qiskit\n",
        "import qiskit_ibm_runtime\n",
        "\n",
        "print(qiskit.version.get_version_info())\n",
        "print(qiskit_ibm_runtime.version.get_version_info())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}