{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "9e40af77-7f0f-4dd6-ab0a-420cf396050e",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Diagonalizzazione quantistica di un Hamiltoniano chimico basata su campioni\"\n",
        "description: \"Utilizza l'algoritmo di diagonalizzazione quantistica basato su campioni per simulare una molecola di azoto utilizzando hardware quantistico rumoroso.\"\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",
        "# Diagonalizzazione quantistica di un Hamiltoniano chimico basata su campioni\n",
        "\n",
        "*Stima di utilizzo: meno di un minuto su un processore Heron r2 (NOTA: questa è solo una stima. Il tempo di esecuzione potrebbe variare)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9701917b",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## Risultati di apprendimento\n",
        "\n",
        "Dopo aver seguito questo tutorial, gli utenti dovrebbero aver compreso:\n",
        "\n",
        "* Come utilizzare [l](/docs/addons/qiskit-addon-sqd) 'add-on SQD Qiskit per approssimare l'energia dello stato fondamentale di un sistema molecolare utilizzando stringhe di bit campionate da un'unità di elaborazione quantistica (QPU).\n",
        "* Come utilizzare [ffsim](https://github.com/qiskit-community/ffsim) per costruire un circuito Jastrow a cluster unitario locale (LUCJ) per la simulazione di chimica quantistica.\n",
        "\n",
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## Prerequisiti\n",
        "\n",
        "Si consiglia agli utenti di familiarizzarsi con i seguenti argomenti prima di seguire questo tutorial:\n",
        "\n",
        "* Chimica quantistica e seconda quantizzazione\n",
        "* Utilizzo della primitiva Sampler per campionare da circuiti quantistici\n",
        "\n",
        "<span id=\"background\" />\n",
        "\n",
        "## Sfondo\n",
        "\n",
        "In questo tutorial mostriamo come eseguire la post-elaborazione di campioni quantistici soggetti a rumore per approssimare lo stato fondamentale della molecola di azoto $\\text{N}_2$ alla lunghezza di legame di equilibrio, utilizzando [l](https://github.com/Qiskit/qiskit-addon-sqd) 'add-on SQD di Qiskit per implementare [l](https://arxiv.org/abs/2405.05068) 'algoritmo di diagonalizzazione quantistica basata su campioni (SQD). Maggiori dettagli sul software sono disponibili nella [documentazione](/docs/addons/qiskit-addon-sqd) corrispondente, che include anche un [semplice esempio](/docs/addons/qiskit-addon-sqd/guides/quickstart) per iniziare.\n",
        "\n",
        "Questo tutorial è consigliato agli utenti che hanno familiarità con la chimica quantistica: in particolare, è richiesta una certa dimestichezza con il calcolo delle energie dello stato fondamentale di una molecola. Per una guida dettagliata al flusso di lavoro, consulta il [corso sull'algoritmo di diagonalizzazione quantistica](/learning/courses/quantum-diagonalization-algorithms).\n",
        "\n",
        "L'SQD è una tecnica che consente di determinare gli autovalori e gli autovettori di operatori quantistici, come l'hamiltoniano di un sistema quantistico, ricorrendo congiuntamente al calcolo quantistico e al calcolo classico distribuito. Il calcolo distribuito classico viene utilizzato per elaborare i campioni ottenuti da un processore quantistico e per proiettare e diagonalizzare un hamiltoniano di riferimento in un sottospazio da essi generato. Un flusso di lavoro basato su SQD prevede le seguenti fasi:\n",
        "\n",
        "1. Scegliere un'ipotesi di circuito e applicarla su un computer quantistico a uno stato di riferimento (in questo caso, lo stato [Hartree-Fock](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method) ).\n",
        "2. Esempi di bitstings dallo stato quantico risultante.\n",
        "3. Eseguire la procedura *di recupero della configurazione auto-consistente* sulle stringhe di bit per ottenere l'approssimazione dello stato fondamentale.\n",
        "\n",
        "È noto che la SQD funziona bene quando l'autostato bersaglio è rado: la funzione d'onda è supportata in un insieme di stati base $\\mathcal{S} = \\{|x\\rangle \\}$ la cui dimensione non aumenta esponenzialmente con la dimensione del problema.\n",
        "\n",
        "<span id=\"quantum-chemistry\" />\n",
        "\n",
        "### Chimica quantistica\n",
        "\n",
        "L'hamiltoniana di un sistema molecolare può essere scritta come\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",
        "dove $h_{pr}$ e $h_{prqs}$ sono numeri complessi chiamati integrali molecolari che possono essere calcolati dalle specifiche della molecola utilizzando un programma informatico. In questa esercitazione, calcoliamo gli integrali utilizzando il pacchetto software [PySCF](https://pyscf.org/) pacchetto software.\n",
        "\n",
        "Per i dettagli su come viene ricavata l'hamiltoniana molecolare, consultare un testo di chimica quantistica (per esempio, *Modern Quantum Chemistry* di Szabo e Ostlund). Per una spiegazione di alto livello di come i problemi di chimica quantistica vengono mappati sui computer quantistici, si veda la lezione [*Mapping Problems to Qubits*](https://youtube.com/watch?v=TyFU6r8uEsE\\&t=900) della Qiskit Global Summer School 2024.\n",
        "\n",
        "<span id=\"local-unitary-cluster-jastrow-lucj-ansatz\" />\n",
        "\n",
        "### Approccio del cluster unitario locale di Jastrow (LUCJ)\n",
        "\n",
        "SQD richiede un ansatz del circuito quantistico da cui estrarre i campioni. In questo tutorial utilizzeremo l'approccio [del cluster unitario locale di Jastrow (LUCJ)](https://pubs.rsc.org/en/content/articlelanding/2023/sc/d3sc02516k) per la sua combinazione di fondamenti fisici e facilità di implementazione hardware. Useremo [ffsim](https://qiskit-community.github.io/ffsim/) per costruire il circuito di ansatz.\n",
        "\n",
        "L'approccio LUCJ si adatta alle QPU con connettività limitata dei qubit. Gli orbitali di spin vengono mappati sui qubit in modo tale che l'ansatz non richieda l'utilizzo di porte SWAP. IBM® L'hardware presenta una topologia dei qubit a reticolo esagonale denso; in tal caso, possiamo adottare uno schema \"a zig-zag\", illustrato di seguito. In questo schema, gli orbitali con lo stesso spin sono associati a qubit con una topologia lineare (cerchi rossi e blu), mentre tra gli orbitali con spin diverso è presente una connessione ogni quattro orbitali spaziali, facilitata da un qubit ancilla (cerchi viola).\n",
        "\n",
        "![Diagramma di mappatura dei Qubit per l'ansatz LUCJ su un reticolo heavy-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",
        "### Ripristino della configurazione auto-coerente\n",
        "\n",
        "La procedura di recupero della configurazione autoconsistente è progettata per estrarre il maggior segnale possibile da campioni quantistici rumorosi. Poiché l'hamiltoniana molecolare conserva il numero di particelle e lo spin Z, ha senso scegliere un'ipotesi di circuito che conservi anche queste simmetrie. Quando viene applicato allo stato Hartree-Fock, lo stato risultante ha un numero di particelle e uno spin Z fissi nell'impostazione senza rumore. Pertanto, le metà di spin $\\alpha$ e di spin $\\beta$ di qualsiasi stringa di bit campionata da questo stato dovrebbero avere lo stesso [peso di Hamming](https://en.wikipedia.org/wiki/Hamming_weight) dello stato Hartree-Fock. A causa della presenza di rumore negli attuali processori quantistici, alcune stringhe di bit misurate violeranno questa proprietà. Una semplice forma di post-selezione consentirebbe di scartare queste stringhe di bit, ma è uno spreco perché queste stringhe di bit potrebbero contenere ancora qualche segnale. La procedura di recupero autoconsistente tenta di recuperare parte del segnale in post-elaborazione. La procedura è iterativa e richiede come input una stima delle occupazioni medie di ciascun orbitale nello stato fondamentale, che viene prima calcolata dai campioni grezzi. La procedura viene eseguita in un ciclo e ogni iterazione prevede le seguenti fasi:\n",
        "\n",
        "1. Per ogni stringa di bit che viola le simmetrie specificate, si capovolgono i suoi bit con una procedura probabilistica volta ad avvicinare la stringa di bit alla stima corrente delle occupazioni orbitali medie, per ottenere una nuova stringa di bit.\n",
        "2. Raccogliere tutte le bitstringhe vecchie e nuove che soddisfano le simmetrie e sotto-campionare sottoinsiemi di dimensioni fisse, scelte in anticipo.\n",
        "3. Per ogni sottoinsieme di stringhe di bit, proiettare l'hamiltoniana nel sottospazio coperto dai vettori base corrispondenti (si veda la [sezione precedente](#quantum-chemistry) per una descrizione di questi vettori base) e calcolare una stima dello stato fondamentale dell'hamiltoniana proiettata su un computer classico.\n",
        "4. Aggiornare la stima delle occupazioni orbitali medie con la stima dello stato fondamentale con l'energia più bassa.\n",
        "\n",
        "<span id=\"sqd-workflow-diagram\" />\n",
        "\n",
        "### Diagramma del flusso di lavoro SQD\n",
        "\n",
        "Il flusso di lavoro SQD è rappresentato nel seguente diagramma:\n",
        "\n",
        "![Diagramma del flusso di lavoro dell'algoritmo 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",
        "## Requisiti\n",
        "\n",
        "Prima di iniziare questa esercitazione, assicuratevi di aver installato quanto segue:\n",
        "\n",
        "* Qiskit SDK v1.0 o versioni successive, con supporto [alla visualizzazione](/docs/api/qiskit/visualization)\n",
        "* Qiskit Runtime v0.22 o successivamente (`pip install qiskit-ibm-runtime`)\n",
        "* Componente aggiuntivo SQD Qiskit v0.11 o versioni successive (`pip install qiskit-addon-sqd`)\n",
        "* ffsim v0.0.75 o versioni successive (`pip install ffsim`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c6e44a31",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## Configura\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",
        "## Esempio di simulatore su piccola scala\n",
        "\n",
        "In questo tutorial cercheremo di trovare un'approssimazione dello stato fondamentale di una molecola di azoto in prossimità della sua distanza di legame di equilibrio. Per prima cosa utilizziamo un piccolo set di basi \" STO-6G \" in modo da poter simulare l'esperimento e assicurarci che funzioni.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "afeb054c",
      "metadata": {},
      "source": [
        "<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",
        "Per prima cosa, identifichiamo la molecola e le sue proprietà.\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": [
        "Prima di costruire il circuito di ansatz LUCJ, eseguiamo un calcolo CCSD nella seguente cella di codice. Le [ampiezze $t_1$ e $t_2$](https://en.wikipedia.org/wiki/Coupled_cluster#Cluster_operator) di questo calcolo saranno utilizzate per inizializzare i parametri dell'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": [
        "Ora utilizziamo [ffsim](https://github.com/qiskit-community/ffsim) per creare il circuito di ansatz. Poiché la nostra molecola presenta uno stato di Hartree-Fock a guscio chiuso, utilizziamo la variante a spin bilanciato dell’ansatz UCJ, [UCJOpSpinBalanced](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced). Abbiamo impostato `optimize=True` nel metodo `from_t_amplitudes` per abilitare la doppia fattorizzazione \"compressa\" delle ampiezze di $t_2$ (per ulteriori dettagli, consultare [la](https://qiskit-community.github.io/ffsim/explanations/lucj.html#Parameter-initialization-from-CCSD) documentazione di ffsim relativa all'approccio LUCJ (Local Unitary Cluster Jastrow)).\n",
        "\n",
        "Poiché l'ansatz LUCJ si adatta alla connettività disponibile della QPU, è necessario inizializzare il backend della QPU prima di creare l'ansatz. Per ora, creeremo un backend generico con una mappa ad accoppiamento esagonale forte e un insieme di gate in cui l'ansatz LUCJ si scompone naturalmente. Quindi, useremo `ffsim.qiskit.generate_lucj_pass_manager` per creare un gestore di passaggi specializzato nella trasposizione dell'approccio LUCJ al backend specificato, secondo lo schema \"a zig-zag\" descritto nella [sezione](#local-unitary-cluster-jastrow-lucj-ansatz) introduttiva dedicata all'approccio LUCJ. Questa funzione utilizza un algoritmo euristico di valutazione per ridurre al minimo gli errori associati al layout selezionato, il che è importante se il backend è una vera QPU o un simulatore con un modello di rumore. Oltre a restituire il gestore dei pass, questa funzione restituisce anche le coppie di accoppiamento alfa-beta che possono essere implementate sull'hardware. Se non è possibile implementare tutte le coppie, viene emesso un avviso.\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",
        "### Fase 2: Ottimizzazione per l'esecuzione su hardware quantistico\n",
        "\n",
        "Successivamente, ottimizziamo il circuito per un hardware specifico. In genere, questa fase prevede l'inizializzazione del backend hardware e di un gestore di passaggi per quel backend. Tuttavia, poiché l'approccio LUCJ è adattato alla connettività hardware, abbiamo già eseguito queste operazioni nella fase precedente. Non resta che eseguire il gestore di passaggi sul circuito per traspilarlo in un circuito ISA che possa essere eseguito direttamente sulla 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",
        "### Passaggio 3: eseguire utilizzando Qiskit primitives\n",
        "\n",
        "Dopo aver ottimizzato il circuito per l'esecuzione su hardware, siamo pronti a eseguirlo sull'hardware di destinazione e a raccogliere campioni per la stima dell'energia dello stato fondamentale. Poiché disponiamo di un solo circuito, utilizzeremo [la modalità di esecuzione](/docs/guides/execution-modes) “ IBM Quantum ” di Compute Service per eseguire il nostro circuito.\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",
        "### Fase 4: Post-elaborazione e restituzione del risultato nel formato classico desiderato\n",
        "\n",
        "Un parametro utile per valutare la qualità dell'output della QPU è il numero di configurazioni valide restituite. Una configurazione valida ha il numero corretto di particelle e lo spin Z, il che significa che la metà destra della stringa di bit ha un peso di Hamming pari al numero di elettroni con spin up, mentre la metà sinistra ha un peso di Hamming pari al numero di elettroni con spin down. La cella seguente calcola la frazione di configurazioni campionate che sono valide.\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": [
        "Tutte le stringhe di bit sono valide perché stiamo effettuando il campionamento del circuito su un simulatore privo di rumore. Se l'operazione viene eseguita su una QPU soggetta a rumore, la frazione sarà inferiore a uno, ma si spera che sia maggiore di quella che ci si aspetterebbe se le stringhe di bit fossero campionate in modo casuale e uniforme, come calcolato nella cella seguente.\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": [
        "Ora stimiamo l'energia di stato fondamentale dell'hamiltoniana utilizzando la funzione `diagonalize_fermionic_hamiltonian` . Questa funzione esegue la procedura di recupero della configurazione autoconsistente per affinare iterativamente i campioni quantistici rumorosi e migliorare la stima dell'energia. Si passa una funzione di callback in modo da poter salvare i risultati intermedi per un'analisi successiva. Vedere la [documentazione dell'API](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) per le spiegazioni degli argomenti di `diagonalize_fermionic_hamiltonian`.\n",
        "\n",
        "Qui, utilizziamo `initial_occupancies` l'argomento per `diagonalize_fermionic_hamiltonian` specificare la configurazione di Hartree-Fock come ipotesi iniziale per le occupazioni orbitali nello stato fondamentale. Questo approccio è sensato per i sistemi in cui lo stato fondamentale ha un supporto significativo sulla configurazione di Hartree-Fock, ma potrebbe non essere appropriato in altre situazioni, anche se metodi computazionali più avanzati potrebbero fornire ipotesi iniziali migliori in questi casi. Specificando è `initial_occupancies` anche possibile eseguire il ripristino della configurazione anche se non sono state campionate configurazioni valide, come può accadere quando si campiona un circuito di grandi dimensioni su una QPU rumorosa. Senza questo argomento, il ripristino della configurazione non andrebbe a buon fine e genererebbe un errore se non fossero state fornite configurazioni valide.\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",
        "#### Visualizza i risultati\n",
        "\n",
        "Il primo grafico mostra che in questa simulazione siamo già vicini `1 mH` alla risposta esatta dopo la prima iterazione (l'accuratezza chimica è generalmente considerata pari a `1 kcal/mol`$\\approx$`1.6 mH`). Si tratta però di un sistema di piccole dimensioni e, poiché i campioni sono privi di rumore, non è necessario recuperare la configurazione. Su un sistema più grande in esecuzione su una QPU soggetta a rumore, potrebbero essere necessarie più iterazioni di recupero della configurazione e la precisione finale potrebbe risultare inferiore. In generale, è possibile migliorare l'energia consentendo un maggior numero di iterazioni di recupero della configurazione oppure aumentando il numero di campioni per lotto.\n",
        "\n",
        "Il secondo grafico mostra l'occupazione media di ciascun orbitale spaziale dopo l'iterazione finale. Possiamo notare che sia gli elettroni con spin-up che quelli con spin-down occupano i primi cinque orbitali con alta probabilità nelle nostre soluzioni.\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",
        "## Esempio di hardware su larga scala\n",
        "\n",
        "Ora eseguiamo un esempio più ampio su un hardware quantistico reale. In questa sede, ricaveremo uno spazio attivo per la molecola di azoto a partire dal set di basi \" cc-pVDZ \".\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "24ca3090-3f3b-4efb-a482-70b2e1b5d062",
      "metadata": {},
      "source": [
        "<span id=\"steps-1-4\" />\n",
        "\n",
        "### Passaggi da 1 a 4\n",
        "\n",
        "Qui riuniamo tutte le fasi in un unico flusso di lavoro su scala più ampia, che viene poi eseguito su hardware quantistico reale.\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",
        "## Passi successivi\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Suggerimenti\">\n",
        "  Se questo lavoro ti è sembrato interessante, potrebbero interessarti anche i seguenti contenuti:\n",
        "\n",
        "  * [Diagonalizzazione quantistica di Krylov basata su campioni di un modello a reticolo fermionico](/docs/tutorials/sample-based-krylov-quantum-diagonalization) – un tutorial correlato che utilizza circuiti di evoluzione temporale anziché un approccio variazionale\n",
        "  * [Scalare i flussi di lavoro chimici SQD con il risolutore Dice](/docs/addons/qiskit-addon-sqd/guides/integrate-dice-solver) : una pagina che illustra come utilizzare il software Dice, più efficiente, per la diagonalizzazione\n",
        "  * [Documentazione dell'API dell'add-on SQD](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) - guida di riferimento per la `diagonalize_fermionic_hamiltonian` funzione\n",
        "  * [*La chimica oltre i limiti della diagonalizzazione esatta su un supercomputer quantistico*](https://www.science.org/doi/10.1126/sciadv.adu9991) : l'articolo su cui si basa questo tutorial\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
}