{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "9e40af77-7f0f-4dd6-ab0a-420cf396050e",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Diagonalización cuántica basada en muestras de un hamiltoniano químico\"\n",
        "description: \"Utiliza el algoritmo de diagonalización cuántica basado en muestras para simular una molécula de nitrógeno utilizando hardware cuántico ruidoso.\"\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",
        "# Diagonalización cuántica basada en muestras de un hamiltoniano químico\n",
        "\n",
        "*Estimación de uso: menos de un minuto en un procesador Heron r2 (NOTA: Esto es sólo una estimación. Su tiempo de ejecución puede variar)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9701917b",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## Resultados del aprendizaje\n",
        "\n",
        "Una vez completado este tutorial, los usuarios deberían comprender:\n",
        "\n",
        "* Cómo utilizar el [complemento SQD Qiskit](/docs/addons/qiskit-addon-sqd) para aproximar la energía del estado fundamental de un sistema molecular utilizando cadenas de bits muestreadas a partir de una unidad de procesamiento cuántico (QPU).\n",
        "* Cómo utilizar [ffsim](https://github.com/qiskit-community/ffsim) para construir un circuito Jastrow de clúster unitario local (LUCJ) para la simulación de química cuántica.\n",
        "\n",
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## Requisitos previos\n",
        "\n",
        "Recomendamos a los usuarios que se familiaricen con los siguientes temas antes de seguir este tutorial:\n",
        "\n",
        "* Química cuántica y segunda cuantización\n",
        "* Uso de la primitiva «Sampler» para obtener muestras de circuitos cuánticos\n",
        "\n",
        "<span id=\"background\" />\n",
        "\n",
        "## En segundo plano\n",
        "\n",
        "En este tutorial, mostramos cómo realizar el posprocesamiento de muestras cuánticas con ruido para aproximar el estado fundamental de la molécula de nitrógeno $\\text{N}_2$ a la longitud de enlace de equilibrio, utilizando el [complemento SQD de Qiskit](https://github.com/Qiskit/qiskit-addon-sqd) para implementar el [algoritmo de diagonalización cuántica basada en muestras (SQD)](https://arxiv.org/abs/2405.05068). En la [documentación](/docs/addons/qiskit-addon-sqd) correspondiente se puede encontrar más información sobre el software, incluido un [ejemplo sencillo](/docs/addons/qiskit-addon-sqd/guides/quickstart) para empezar.\n",
        "\n",
        "Este tutorial está recomendado para usuarios que estén familiarizados con la química cuántica; en concreto, con el cálculo de las energías del estado fundamental de una molécula. Para obtener una guía detallada sobre el flujo de trabajo, consulta el [curso](/learning/courses/quantum-diagonalization-algorithms) sobre el algoritmo de diagonalización cuántica.\n",
        "\n",
        "El SQD es una técnica para hallar los valores propios y los vectores propios de operadores cuánticos, como el hamiltoniano de un sistema cuántico, mediante la combinación de la computación cuántica y la computación clásica distribuida. La computación distribuida clásica se utiliza para procesar muestras obtenidas de un procesador cuántico, así como para proyectar y diagonalizar un hamiltoniano objetivo en un subespacio que estas abarcan. Un flujo de trabajo basado en SQD consta de los siguientes pasos:\n",
        "\n",
        "1. Elija un ansatz de circuito y aplíquelo en un ordenador cuántico a un estado de referencia (en este caso, el estado [Hartree-Fock](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method) ).\n",
        "2. Muestras de cadenas de bits del estado cuántico resultante.\n",
        "3. Ejecuta el procedimiento *de recuperación de la configuración autoconsistente* en las cadenas de bits para obtener la aproximación del estado fundamental.\n",
        "\n",
        "Se sabe que SQD funciona bien cuando el estado propio objetivo es disperso: la función de onda está soportada en un conjunto de estados base $\\mathcal{S} = \\{|x\\rangle \\}$ cuyo tamaño no aumenta exponencialmente con el tamaño del problema.\n",
        "\n",
        "<span id=\"quantum-chemistry\" />\n",
        "\n",
        "### Química cuántica\n",
        "\n",
        "El Hamiltoniano de un sistema molecular puede escribirse como\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",
        "donde $h_{pr}$ y $h_{prqs}$ son números complejos denominados integrales moleculares que pueden calcularse a partir de la especificación de la molécula utilizando un programa informático. En este tutorial, calculamos las integrales utilizando el paquete de software [PySCF](https://pyscf.org/) paquete de software.\n",
        "\n",
        "Para más detalles sobre cómo se obtiene el hamiltoniano molecular, consulte un libro de texto sobre química cuántica (por ejemplo, *Modern Quantum Chemistry* de Szabo y Ostlund). Para una explicación de alto nivel de cómo los problemas de química cuántica se trasladan a los ordenadores cuánticos, consulte la conferencia [*Mapping Problems to Qubits*](https://youtube.com/watch?v=TyFU6r8uEsE\\&t=900) de la Escuela Mundial de Verano Qiskit 2024.\n",
        "\n",
        "<span id=\"local-unitary-cluster-jastrow-lucj-ansatz\" />\n",
        "\n",
        "### Enfoque del clúster unitario local de Jastrow (LUCJ)\n",
        "\n",
        "El SQD requiere un ansatz de circuito cuántico del que extraer muestras. En este tutorial utilizaremos el enfoque [del clúster unitario local de Jastrow (LUCJ)](https://pubs.rsc.org/en/content/articlelanding/2023/sc/d3sc02516k), debido a su combinación de fundamentación física y facilidad de implementación en hardware. Utilizaremos [ffsim](https://qiskit-community.github.io/ffsim/) para construir el circuito de aproximación.\n",
        "\n",
        "El enfoque LUCJ se adapta a las QPU con conectividad de qubits limitada. Los orbitales de espín se asignan a qubits de tal manera que el ansatz no requiere el uso de puertas SWAP. IBM® El hardware presenta una topología de qubits de red hexagonal densa, en cuyo caso podemos adoptar un patrón «en zigzag», tal y como se muestra a continuación. En este esquema, los orbitales con el mismo espín se asignan a qubits con una topología lineal (círculos rojos y azules), y cada cuatro orbitales espaciales hay una conexión entre orbitales de espín diferente, facilitada por un qubit auxiliar (círculos morados).\n",
        "\n",
        "![Diagrama de mapeo Qubit para el ansatz LUCJ en una red hexagonal pesada](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",
        "### Recuperación de configuración autoconsistente\n",
        "\n",
        "El procedimiento de recuperación de la configuración autoconsistente está diseñado para extraer la mayor cantidad de señal posible de muestras cuánticas ruidosas. Dado que el Hamiltoniano molecular conserva el número de partícula y el espín Z, tiene sentido elegir un ansatz de circuito que también conserve estas simetrías. Cuando se aplica al estado Hartree-Fock, el estado resultante tiene un número de partícula y un espín Z fijos en la configuración sin ruido. Por lo tanto, las mitades spin- $\\alpha$ y spin- $\\beta$ de cualquier cadena de bits muestreada desde este estado deben tener el mismo [peso Hamming](https://en.wikipedia.org/wiki/Hamming_weight) que en el estado Hartree-Fock. Debido a la presencia de ruido en los procesadores cuánticos actuales, algunas cadenas de bits medidas violarán esta propiedad. Una forma sencilla de postselección descartaría estas cadenas de bits, pero esto es un desperdicio porque esas cadenas de bits aún podrían contener alguna señal. El procedimiento de recuperación autoconsistente intenta recuperar parte de esa señal en el postprocesado. El procedimiento es iterativo y requiere como entrada una estimación de las ocupaciones medias de cada orbital en el estado fundamental, que se calcula primero a partir de las muestras brutas. El procedimiento se ejecuta en un bucle, y cada iteración tiene los siguientes pasos:\n",
        "\n",
        "1. Para cada cadena de bits que viole las simetrías especificadas, voltea sus bits con un procedimiento probabilístico diseñado para acercar la cadena de bits a la estimación actual de las ocupaciones orbitales medias, para obtener una nueva cadena de bits.\n",
        "2. Recoge todas las cadenas de bits antiguas y nuevas que cumplan las simetrías, y submuestrea subconjuntos de un tamaño fijo, elegido de antemano.\n",
        "3. Para cada subconjunto de cadenas de bits, proyecte el Hamiltoniano en el subespacio abarcado por los vectores base correspondientes (véase la [sección anterior](#quantum-chemistry) para una descripción de estos vectores base), y calcule una estimación del estado base del Hamiltoniano proyectado en un ordenador clásico.\n",
        "4. Actualizar la estimación de las ocupaciones orbitales medias con la estimación del estado fundamental con la energía más baja.\n",
        "\n",
        "<span id=\"sqd-workflow-diagram\" />\n",
        "\n",
        "### Diagrama del flujo de trabajo SQD\n",
        "\n",
        "El flujo de trabajo de SQD se representa en el siguiente diagrama:\n",
        "\n",
        "![Diagrama de flujo de trabajo del 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",
        "## Requisitos\n",
        "\n",
        "Antes de empezar este tutorial, asegúrate de que tienes instalado lo siguiente:\n",
        "\n",
        "* Qiskit SDK v1.0 o posterior, con soporte [de visualización](/docs/api/qiskit/visualization)\n",
        "* Qiskit Runtime v0.22 o posterior (`pip install qiskit-ibm-runtime`)\n",
        "* Complemento SQD Qiskit v0.11 o posterior (`pip install qiskit-addon-sqd`)\n",
        "* ffsim v0.0.75 o posterior (`pip install ffsim`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c6e44a31",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## Configuración\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",
        "## Ejemplo de simulador a pequeña escala\n",
        "\n",
        "En este tutorial, hallaremos una aproximación al estado fundamental de una molécula de nitrógeno cerca de su distancia de enlace de equilibrio. En primer lugar, utilizamos un pequeño conjunto de bases « STO-6G » para poder simular el experimento y asegurarnos de que funciona.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "afeb054c",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### Paso 1: Asignar entradas clásicas a un problema cuántico\n",
        "\n",
        "En primer lugar, especificamos la molécula y sus propiedades.\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": [
        "Antes de construir el circuito LUCJ ansatz, primero realizamos un cálculo CCSD en la siguiente celda de código. Las [amplitudes $t_1$ y $t_2$](https://en.wikipedia.org/wiki/Coupled_cluster#Cluster_operator) de este cálculo se utilizarán para inicializar los parámetros del 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": [
        "Ahora utilizamos [ffsim](https://github.com/qiskit-community/ffsim) para crear el circuito de aproximación. Dado que nuestra molécula presenta un estado de Hartree-Fock de capa cerrada, utilizamos la variante de equilibrio de espín del enfoque UCJ, [UCJOpSpinBalanced](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced). En el `from_t_amplitudes` método hemos establecido `optimize=True` para habilitar la doble factorización «comprimida» de las amplitudes de $t_2$ (véase [el enfoque del clúster unitario local de Jastrow (LUCJ)](https://qiskit-community.github.io/ffsim/explanations/lucj.html#Parameter-initialization-from-CCSD) en la documentación de ffsim para más detalles).\n",
        "\n",
        "Dado que el ansatz LUCJ se adapta a la conectividad disponible de la QPU, debemos inicializar el backend de la QPU antes de crear el ansatz. Por ahora, crearemos un backend genérico con un mapa de acoplamiento hexagonal fuerte y un conjunto de puertas en el que el enfoque LUCJ se descompone de forma natural. A continuación, utilizaremos `ffsim.qiskit.generate_lucj_pass_manager` para crear un gestor de pases especializado en la transpilación del enfoque LUCJ al backend especificado, siguiendo el esquema «en zigzag» descrito en la [sección de antecedentes sobre el enfoque LUCJ](#local-unitary-cluster-jastrow-lucj-ansatz). Esta función utiliza una heurística de puntuación para minimizar los errores asociados al diseño seleccionado, lo cual es importante si tu backend es una QPU real o un simulador con un modelo de ruido. Además de devolver el gestor de pases, esta función también devuelve los pares de acoplamiento alfa-beta que se pueden implementar en el hardware. Si no es posible implementar todos los pares, se emite una advertencia.\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",
        "### Paso 2: Optimizar para la ejecución en hardware cuántico\n",
        "\n",
        "A continuación, optimizamos el circuito para un hardware específico. Por lo general, este paso consiste en inicializar el backend de hardware y un gestor de pases para dicho backend. Sin embargo, dado que el enfoque LUCJ se adapta a la conectividad del hardware, ya hemos llevado a cabo estas acciones en el paso anterior. Lo único que queda por hacer es ejecutar el gestor de pases en el circuito para transpilarlo a un circuito ISA que pueda ejecutarse directamente en la 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",
        "### Paso 3: Ejecutar utilizando Qiskit primitives\n",
        "\n",
        "Después de optimizar el circuito para su ejecución en hardware, estamos listos para ejecutarlo en el hardware de destino y recoger muestras para la estimación de la energía del estado de tierra. Como sólo tenemos un circuito, utilizaremos el [modo de ejecución Job](/docs/guides/execution-modes) de Qiskit Runtime y ejecutaremos nuestro 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",
        "### Paso 4: Procesamiento posterior y devolución del resultado en el formato clásico deseado\n",
        "\n",
        "Una métrica útil para evaluar la calidad de la salida de la QPU es el número de configuraciones válidas devueltas. Una configuración válida tiene el número correcto de partículas y el espín Z, lo que significa que la mitad derecha de la cadena de bits tiene un peso de Hamming igual al número de electrones con espín hacia arriba, y la mitad izquierda tiene un peso de Hamming igual al número de electrones con espín hacia abajo. La siguiente celda calcula la fracción de configuraciones muestreadas que son válidas.\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": [
        "Todas las cadenas de bits son válidas porque estamos realizando un muestreo del circuito en un simulador sin ruido. Cuando se ejecuta en una QPU ruidosa, la fracción será inferior a uno, pero es de esperar que sea mayor que la fracción que cabría esperar si las cadenas de bits se muestrearan de forma aleatoria y uniforme, tal y como se calcula en la siguiente celda.\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": [
        "Ahora, estimamos la energía del estado fundamental del Hamiltoniano utilizando la función `diagonalize_fermionic_hamiltonian` . Esta función realiza el procedimiento de recuperación de la configuración autoconsistente para refinar iterativamente las muestras cuánticas ruidosas con el fin de mejorar la estimación de la energía. Pasamos una función callback para poder guardar los resultados intermedios para su posterior análisis. Consulte [la documentación de la API](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) para obtener explicaciones sobre los argumentos de `diagonalize_fermionic_hamiltonian`.\n",
        "\n",
        "Aquí, utilizamos el `initial_occupancies` argumento para `diagonalize_fermionic_hamiltonian` especificar la configuración de Hartree-Fock como la estimación inicial para las ocupaciones orbitales en el estado fundamental. Este enfoque es adecuado para sistemas en los que el estado fundamental tiene un apoyo significativo en la configuración de Hartree-Fock, pero puede que no sea apropiado en otras situaciones, aunque métodos computacionales más avanzados podrían proporcionar mejores estimaciones iniciales en esos casos. Especificar `initial_occupancies` también permite que la recuperación de la configuración se ejecute incluso si no se han muestreado configuraciones válidas, como puede ocurrir al muestrear un circuito grande en una QPU ruidosa. Sin este argumento, la recuperación de la configuración fallaría y generaría un error si no se proporcionaran configuraciones válidas.\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",
        "#### Visualizar los resultados\n",
        "\n",
        "El primer gráfico muestra que, en esta simulación, ya nos encontramos muy `1 mH` cerca de la respuesta exacta tras la primera iteración (por lo general, se acepta que la precisión química es de `1 kcal/mol`$\\approx$`1.6 mH`). Sin embargo, se trata de un sistema pequeño y, dado que las muestras no contienen ruido, no es necesario recuperar la configuración. En un sistema más grande que se ejecute en una QPU ruidosa, es posible que se necesiten varias iteraciones de recuperación de la configuración y que la precisión final sea menor. Por lo general, se puede mejorar la energía permitiendo más iteraciones de recuperación de la configuración o aumentando el número de muestras por lote.\n",
        "\n",
        "El segundo gráfico muestra la ocupación media de cada orbital espacial tras la iteración final. Podemos ver que tanto los electrones de espín arriba como los de espín abajo ocupan los cinco primeros orbitales con alta probabilidad en nuestras soluciones.\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",
        "## Ejemplo de hardware a gran escala\n",
        "\n",
        "Ahora ejecutamos un ejemplo más amplio en hardware cuántico real. En este caso, derivaremos un espacio activo para la molécula de nitrógeno a partir del conjunto de bases « cc-pVDZ ».\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "24ca3090-3f3b-4efb-a482-70b2e1b5d062",
      "metadata": {},
      "source": [
        "<span id=\"steps-1-4\" />\n",
        "\n",
        "### Pasos 1 a 4\n",
        "\n",
        "Aquí reunimos todos los pasos en un único flujo de trabajo a mayor escala, que luego se ejecuta en hardware cuántico real.\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",
        "## Próximos pasos\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Recomendaciones\">\n",
        "  Si te ha parecido interesante este trabajo, quizá te interese el siguiente material:\n",
        "\n",
        "  * [Diagonalización cuántica de Krylov basada en muestras de un modelo de red fermiónica](/docs/tutorials/sample-based-krylov-quantum-diagonalization) : tutorial relacionado que utiliza circuitos de evolución temporal en lugar de un enfoque variacional\n",
        "  * [Amplía los flujos de trabajo químicos de SQD con el solucionador Dice](/docs/addons/qiskit-addon-sqd/guides/integrate-dice-solver) : una página que explica cómo utilizar el software Dice, más eficiente, para la diagonalización\n",
        "  * [Documentación de la API del complemento SQD](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) : referencia de la `diagonalize_fermionic_hamiltonian` función\n",
        "  * [*La química más allá de la escala de la diagonalización exacta en un superordenador cuántico*](https://www.science.org/doi/10.1126/sciadv.adu9991) : el artículo en el que se basa este 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
}