{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "de193554-b271-4295-95e4-8904f0f6ee8a",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Geometria molecolare\"\n",
        "description: \"In questa lezione variamo la geometria di una molecola semplice, minimizzando l'energia ad ogni passo.\"\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",
        "# Determinazione della geometria molecolare\n",
        "\n",
        "Nella sezione precedente, abbiamo implementato il VQE per determinare l'energia di stato fondamentale di una molecola. Questo è un uso valido dell'informatica quantistica, ma ancora più utile sarebbe determinare la struttura di una molecola.\n",
        "\n",
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "## Fase 1: mappare gli input classici su un problema quantistico\n",
        "\n",
        "Restando al nostro esempio di base dell'idrogeno biatomico, l'unico parametro geometrico da variare è la lunghezza del legame. A tal fine, si procede come in precedenza, ma utilizzando una variabile nella costruzione iniziale della molecola (una lunghezza di legame, *x*, nell'argomento). Si tratta di una modifica abbastanza semplice, ma richiede che la variabile sia inclusa nelle funzioni durante tutto il processo, poiché inizia nella costruzione dell'hamiltoniana fermionica e si propaga attraverso la mappatura e infine alla funzione di costo.\n",
        "\n",
        "Per prima cosa, carichiamo alcuni dei pacchetti utilizzati in precedenza e definiamo la funzione 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": [
        "Ora, per definire la nostra hamiltoniana, useremo PySCF esattamente come nell'esempio precedente, ma ora includeremo una variabile, `x`, per svolgere il ruolo di distanza interatomica. Questo restituirà l'energia del nucleo, l'energia di un singolo elettrone e l'energia di due elettroni come in precedenza.\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": [
        "Ricordiamo che la costruzione di cui sopra consiste nel creare un'hamiltoniana fermionica basata sulla specie atomica, sulla geometria e sugli orbitali elettronici. Di seguito, mappiamo questa hamiltoniana fermionica sugli operatori di Pauli. Questa funzione `build_hamiltonian` includerà anche una variabile geometrica come argomento.\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": [
        "Caricheremo i pacchetti rimanenti per l'esecuzione di VQE stesso, come l'ansatz di efficient\\_su2 e i minimizzatori di 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": [
        "Definiremo ancora una volta la funzione di costo, ma questa prende sempre come argomento un'hamiltoniana completamente costruita e mappata, quindi non cambia nulla in questa funzione.\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",
        "## Fase 2: Ottimizzazione del problema per l'esecuzione quantistica\n",
        "\n",
        "Poiché l'hamiltoniana cambierà a ogni nuova geometria, il transpiling dell'operatore cambierà a ogni passo. Possiamo tuttavia definire un gestore di passaggi generale da applicare in ogni fase, specifico per l'hardware che vogliamo utilizzare.\n",
        "\n",
        "In questo caso utilizzeremo il backend meno trafficato disponibile. Utilizzeremo quel backend come modello per il nostro AerSimulator,, consentendo al nostro simulatore di imitare, ad esempio, il comportamento del rumore del backend reale. Questi modelli di rumore non sono perfetti, ma possono aiutare a capire cosa aspettarsi dall'hardware reale.\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": [
        "Importiamo il pass manager e i relativi pacchetti per aiutarci a ottimizzare il nostro circuito. Questo passo, e quello che lo precede, sono indipendenti dall'hamiltoniana e quindi sono invariati rispetto alla lezione precedente.\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",
        "## Fase 3: Eseguire utilizzando le primitive di \" IBM Quantum \".\n",
        "\n",
        "Nel blocco di codice riportato di seguito, definiamo un array per memorizzare i risultati ottenuti in ciascuna fase del nostro processo di calcolo della distanza interatomica $x$. Abbiamo scelto l'intervallo $x$ sulla base delle nostre conoscenze relative al valore sperimentale della lunghezza di legame all'equilibrio: 0.74 Angstrom. Per prima cosa lo eseguiremo su un simulatore, quindi importeremo il nostro Estimator ( BackendEstimator ) da `qiskit.primitives`. Per ogni fase geometrica, costruiamo l'hamiltoniano ed eseguiamo un certo numero di passaggi di ottimizzazione (in questo caso 500) utilizzando l'ottimizzatore \"cobyla\". Ad ogni fase della geometria, memorizziamo sia l'energia totale che l'energia elettronica. A causa dell'elevato numero di fasi dell'ottimizzatore, l'operazione potrebbe richiedere un'ora o più. Potresti voler modificare i dati inseriti qui sotto per ridurre il tempo necessario.\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": [
        "I risultati di questa elaborazione vengono discussi più avanti nella sezione dedicata alla post-elaborazione; per ora, basti notare che la simulazione ha avuto esito positivo. Ora sei pronto per eseguire il programma su hardware reale. Imposteremo la resilienza su `1`, indicando che verrà utilizzata la mitigazione degli errori TREX. Ora che stiamo lavorando con hardware reale, utilizzeremo il servizio di calcolo IBM Quantum e le primitive IBM Quantum. Si noti che sia il ciclo \"for\" relativo alla geometria sia le diverse prove variazionali si trovano all'interno della sessione.\n",
        "\n",
        "Poiché i costi e i tempi di esecuzione sono limitati, abbiamo ridotto il numero di fasi della geometria e dell'ottimizzatore. Assicuratevi di adattare questi passaggi in base ai vostri obiettivi di precisione e ai vostri limiti di tempo.\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",
        "## Fase 4: Post-elaborazione\n",
        "\n",
        "Sia per il simulatore che per l'hardware reale, possiamo tracciare le energie di stato fondamentale calcolate per ogni distanza interatomica e vedere dove si raggiunge l'energia più bassa. Questa dovrebbe essere la distanza interatomica riscontrata in natura, e in effetti è vicina. Una curva più dolce potrebbe essere ottenuta provando altre ansaetze, ottimizzatori ed eseguendo il calcolo più volte a ogni passo di geometria e calcolando la media su diverse condizioni iniziali casuali.\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": [
        "Si noti che il semplice aumento del numero di passi di ottimizzazione non è in grado di migliorare i risultati del simulatore, poiché tutte le ottimizzazioni convergono alla tolleranza richiesta in un numero di iterazioni inferiore al massimo.\n",
        "\n",
        "I risultati ottenuti dall'hardware reale sono paragonabili, a parte una gamma leggermente diversa di valori campionati.\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": [
        "Oltre a prevedere una lunghezza di legame H2 di 0.74 Angstrom, l'energia totale dovrebbe essere di -1.17 Hartrees. Vediamo che i risultati dell'hardware reale si sono avvicinati di più a questi valori rispetto al simulatore. Questo è probabilmente dovuto al fatto che il rumore era presente (o simulato) in entrambi i casi, ma solo nel caso dell'hardware reale è stata impiegata la mitigazione degli errori.\n",
        "\n",
        "<span id=\"closing\" />\n",
        "\n",
        "### Chiusura\n",
        "\n",
        "Si conclude così il nostro corso sulle VQE per la chimica quantistica. Se siete interessati a comprendere alcune delle teorie dell'informazione utilizzate nell'informatica quantistica, date un'occhiata al corso di John Watrous sui [fondamenti dell'informazione quantistica](/learning/courses/basics-of-quantum-information). Per un altro breve esempio di flusso di lavoro VQE, si veda il nostro [tutorial sulla stima dell'energia allo stato solido della catena di Heisenberg con VQE](/docs/tutorials/spin-chain-vqe). Oppure sfogliate le [esercitazioni](/docs/tutorials) e i [corsi](/learning) per trovare altro materiale didattico sulle più recenti tecnologie di calcolo quantistico.\n",
        "\n",
        "Non dimenticate di sostenere l'esame di questo corso. Un punteggio pari o superiore all'80% vi farà guadagnare un badge Credly, che vi verrà automaticamente inviato via e-mail. Grazie per essere parte della rete 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
}