{ "cells": [ { "cell_type": "markdown", "id": "ee388de5-2507-48ff-97a6-6777698d6256", "metadata": {}, "source": [ "---\n", "title: \"Diagonalización cuántica de Krylov basada en muestras de un modelo de red fermiónica\"\n", "description: \"Utiliza el algoritmo de diagonalización cuántica basado en muestras para simular el modelo de Anderson de impureza única utilizando hardware cuántico ruidoso.\"\n", "---\n", "\n", "{/* cspell:ignore fontdict fontsize nocc SQKD DMRG textrm varepsilon vecs pqrs ijkl */}\n", "\n", "\n", "\n", "# Diagonalización cuántica de Krylov basada en muestras de un modelo de red fermiónica\n", "\n", "*Estimación de uso: Nueve segundos 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": "a1b2c3d4-e5f6-7890-abcd-ef1234567890", "metadata": {}, "source": [ "\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](https://github.com/Qiskit/qiskit-addon-sqd) para calcular la energía del estado fundamental de un modelo de red mediante cadenas de bits obtenidas de una unidad de procesamiento cuántico (QPU).\n", "* Cómo utilizar [ffsim](https://github.com/qiskit-community/ffsim) para construir circuitos de evolución temporal para la simulación fermiónica.\n", "* Cómo combinar muestras de varios circuitos para su posprocesamiento con el algoritmo de diagonalización de Krylov basado en muestras (SKQD).\n", "\n", "\n", "\n", "## Requisitos previos\n", "\n", "Recomendamos a los usuarios que se familiaricen con los siguientes temas antes de seguir este tutorial:\n", "\n", "* [Diagonalización cuántica basada en muestras de un hamiltoniano químico](/docs/tutorials/sample-based-quantum-diagonalization)\n", "* [Diagonalización cuántica de Krylov de los hamiltonianos de red](/docs/tutorials/krylov-quantum-diagonalization)\n", "* [Qiskit primitives](/docs/guides/primitives)\n", "\n" ] }, { "cell_type": "markdown", "id": "dc5cc74e-06bf-45ac-a69b-81778138e08f", "metadata": {}, "source": [ "\n", "\n", "## En segundo plano\n", "\n", "Este tutorial muestra cómo utilizar la diagonalización cuántica basada en muestras (SQD) para estimar la energía del estado fundamental de un modelo fermiónico de celosía. En concreto, estudiamos el modelo unidimensional de una sola impureza de Anderson (SIAM), que se utiliza para describir las impurezas magnéticas incrustadas en los metales.\n", "\n", "Este tutorial sigue un flujo de trabajo similar al del tutorial relacionado [Diagonalización cuántica basada en muestras de un Hamiltoniano de química](/docs/tutorials/sample-based-quantum-diagonalization). Sin embargo, una diferencia clave radica en cómo se construyen los circuitos cuánticos. El otro tutorial utiliza un ansatz variacional heurístico, que resulta atractivo para los hamiltonianos químicos con millones de términos de interacción potenciales. Por otro lado, este tutorial utiliza circuitos que aproximan la evolución temporal por el Hamiltoniano. Tales circuitos pueden ser profundos, lo que hace que este enfoque sea mejor para aplicaciones a modelos reticulares. Los vectores de estado preparados por estos circuitos forman la base de un [subespacio de Krylov](https://en.wikipedia.org/wiki/Krylov_subspace) y, como resultado, el algoritmo converge de forma demostrable y eficiente al estado base, bajo supuestos adecuados.\n", "\n", "El enfoque utilizado en este tutorial puede verse como una combinación de las técnicas utilizadas en SQD y [la diagonalización cuántica de Krylov (KQD)](https://arxiv.org/abs/2407.14431). El enfoque combinado se denomina a veces diagonalización cuántica de Krylov basada en muestras (SQKD). Véase [Krylov quantum diagonalization of lattice Hamiltonians](/docs/tutorials/krylov-quantum-diagonalization) para un tutorial sobre el método KQD.\n", "\n", "Este tutorial se basa en el trabajo [\"Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization\"](https://arxiv.org/abs/2501.09702), al que puede consultarse para más detalles.\n", "\n", "\n", "\n", "### Modelo de Anderson de impureza única (SIAM)\n", "\n", "El Hamiltoniano SIAM unidimensional es una suma de tres términos:\n", "\n", "$$\n", "H = H_{\\textrm{imp}}+ H_\\textrm{bath} + H_\\textrm{hyb},\n", "$$\n", "\n", "donde\n", "\n", "$$\n", "\\begin{align*}\n", " H_\\textrm{imp} &= \\varepsilon \\left( \\hat{n}_{d\\uparrow} + \\hat{n}_{d\\downarrow} \\right) + U \\hat{n}_{d\\uparrow}\\hat{n}_{d\\downarrow}, \\\\\n", " H_\\textrm{bath} &= -t \\sum_{\\substack{\\mathbf{j} = 0\\\\ \\sigma\\in \\{\\uparrow, \\downarrow\\}}}^{L-1} \\left(\\hat{c}^\\dagger_{\\mathbf{j}, \\sigma}\\hat{c}_{\\mathbf{j}+1, \\sigma} + \\hat{c}^\\dagger_{\\mathbf{j}+1, \\sigma}\\hat{c}_{\\mathbf{j}, \\sigma} \\right), \\\\\n", " H_\\textrm{hyb} &= V\\sum_{\\sigma \\in \\{\\uparrow, \\downarrow \\}} \\left(\\hat{d}^\\dagger_\\sigma \\hat{c}_{0, \\sigma} + \\hat{c}^\\dagger_{0, \\sigma} \\hat{d}_{\\sigma} \\right).\n", "\\end{align*}\n", "$$\n", "\n", "Aquí, $c^\\dagger_{\\mathbf{j},\\sigma}/c_{\\mathbf{j},\\sigma}$ son los operadores fermiónicos de creación/aniquilación para el sitio de baño $\\mathbf{j}^{\\textrm{th}}$ con espín $\\sigma$, $\\hat{d}^\\dagger_{\\sigma}/\\hat{d}_{\\sigma}$ son operadores de creación/aniquilación para el modo de impureza, y $\\hat{n}_{d\\sigma} = \\hat{d}^\\dagger_{\\sigma} \\hat{d}_{\\sigma}$. $t$, $U$, y $V$ son números reales que describen las interacciones de salto, in situ, de hibridación, y $\\varepsilon$ es un número real que especifica el potencial químico.\n", "\n", "Nótese que el Hamiltoniano es una instancia específica del Hamiltoniano genérico de interacción-electrón,\n", "\n", "$$\n", "\\begin{align*}\n", " H &= \\sum_{\\substack{p, q \\\\ \\sigma}} h_{pq} \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{q\\sigma} + \\sum_{\\substack{p, q, r, s \\\\ \\sigma \\tau}} \\frac{h_{pqrs}}{2} \\hat{a}^\\dagger_{p\\sigma} \\hat{a}^\\dagger_{q\\tau} \\hat{a}_{s\\tau} \\hat{a}_{r\\sigma} \\\\\n", " &= H_1 + H_2,\n", "\\end{align*}\n", "$$\n", "\n", "donde $H_1$ consiste en términos de un cuerpo, que son cuadráticos en los operadores fermiónicos de creación y aniquilación, y $H_2$ consiste en términos de dos cuerpos, que son cuárticos. Para el SIAM,\n", "\n", "$$\n", "H_2 = U \\hat{n}_{d\\uparrow}\\hat{n}_{d\\downarrow}\n", "$$\n", "\n", "y $H_1$ contiene el resto de términos del Hamiltoniano. Para representar el Hamiltoniano programáticamente, almacenamos la matriz $h_{pq}$ y el tensor $h_{pqrs}$.\n", "\n", "\n", "\n", "### Bases de posición y momento\n", "\n", "Debido a la simetría traslacional aproximada en $H_\\textrm{bath}$, no esperamos que el estado base sea disperso en la base de posición (la base orbital en la que se especifica el Hamiltoniano más arriba). El rendimiento de SQD sólo está garantizado si el estado base es disperso, es decir, sólo tiene un peso significativo en un pequeño número de estados base computacionales. Para mejorar la dispersión del estado base, realizamos la simulación en la base orbital en la que $H_\\textrm{bath}$ es diagonal. Llamamos a esta base la *base del momento*. Dado que $H_\\textrm{bath}$ es un Hamiltoniano fermiónico cuadrático, puede diagonalizarse eficientemente mediante una rotación orbital.\n", "\n", "\n", "\n", "### Evolución temporal aproximada mediante el hamiltoniano\n", "\n", "Para aproximar la evolución temporal del hamiltoniano, utilizamos una descomposición de Trotter-Suzuki de segundo orden,\n", "\n", "$$\n", " e^{-i \\Delta t H} \\approx e^{-i\\frac{\\Delta t}{2} H_2} e^{-i\\Delta t H_1} e^{-i\\frac{\\Delta t}{2} H_2}.\n", "$$\n", "\n", "Bajo la [transformación Jordan-Wigner](https://en.wikipedia.org/wiki/Jordan%E2%80%93Wigner_transformation), la evolución temporal por $H_2$ equivale a una única puerta [CPhase](/docs/api/qiskit/qiskit.circuit.library.CPhaseGate) entre los orbitales spin-up y spin-down en el sitio de la impureza. Dado que $H_1$ es un Hamiltoniano fermiónico cuadrático, la evolución temporal por $H_1$ equivale a una rotación orbital.\n", "\n", "Los estados base de Krylov $\\{ |\\psi_k\\rangle \\}_{k=0}^{D-1}$, donde $D$ es la dimensión del subespacio de Krylov, se forman mediante la aplicación repetida de un único paso de Trotter, por lo que\n", "\n", "$$\n", " |\\psi_k\\rangle \\approx \\left[e^{-i\\frac{\\Delta t}{2} H_2} e^{-i\\Delta t H_1} e^{-i\\frac{\\Delta t}{2} H_2} \\right]^k\\ket{\\psi_0}.\n", "$$\n", "\n", "En el siguiente flujo de trabajo basado en SQD, tomaremos muestras de este conjunto de circuitos y posprocesaremos el conjunto combinado de cadenas de bits con SQD. Este enfoque contrasta con el utilizado en el tutorial relacionado [Sample-based quantum diagonalization of a chemistry Hamiltonian](/docs/tutorials/sample-based-quantum-diagonalization), donde las muestras se extrajeron de un único circuito variacional heurístico.\n", "\n" ] }, { "cell_type": "markdown", "id": "df66d697-102c-4a2c-80f3-f67fdda05573", "metadata": {}, "source": [ "\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.72 o posterior (`pip install ffsim`)\n", "\n" ] }, { "cell_type": "markdown", "id": "c3d4e5f6-a7b8-9012-cdef-123456789012", "metadata": {}, "source": [ "\n", "\n", "## Ejemplo de simulador a pequeña escala\n", "\n" ] }, { "cell_type": "markdown", "id": "8540487a-8033-49c2-9f30-022336105f64", "metadata": {}, "source": [ "\n", "\n", "### Paso 1: Asignar el problema a un circuito cuántico\n", "\n" ] }, { "cell_type": "markdown", "id": "4e6652bd-c97d-4f2f-98a7-d44857805cf1", "metadata": {}, "source": [ "En primer lugar, generamos el Hamiltoniano SIAM en la base de posición. El hamiltoniano está representado por la matriz $h_{pq}$ y el tensor $h_{pqrs}$. Luego, lo rotamos a la base del momento. En la base de posición, colocamos la impureza en el primer sitio. Sin embargo, cuando rotamos a la base de momento, movemos la impureza a un sitio central para facilitar las interacciones con otros orbitales.\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "b5cb9c28-4721-4141-8665-96885038e210", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pyscf.fci\n", "\n", "\n", "def siam_hamiltonian(\n", " norb: int,\n", " hopping: float,\n", " onsite: float,\n", " hybridization: float,\n", " chemical_potential: float,\n", ") -> tuple[np.ndarray, np.ndarray]:\n", " \"\"\"Hamiltonian for the single-impurity Anderson model.\"\"\"\n", " # Place the impurity on the first site\n", " impurity_orb = 0\n", "\n", " # One body matrix elements in the \"position\" basis\n", " h1e = np.zeros((norb, norb))\n", " np.fill_diagonal(h1e[:, 1:], -hopping)\n", " np.fill_diagonal(h1e[1:, :], -hopping)\n", " h1e[impurity_orb, impurity_orb + 1] = -hybridization\n", " h1e[impurity_orb + 1, impurity_orb] = -hybridization\n", " h1e[impurity_orb, impurity_orb] = chemical_potential\n", "\n", " # Two body matrix elements in the \"position\" basis\n", " h2e = np.zeros((norb, norb, norb, norb))\n", " h2e[impurity_orb, impurity_orb, impurity_orb, impurity_orb] = onsite\n", "\n", " return h1e, h2e\n", "\n", "\n", "def momentum_basis(norb: int) -> np.ndarray:\n", " \"\"\"Get the orbital rotation to change from the position to the momentum basis.\"\"\"\n", " n_bath = norb - 1\n", "\n", " # Orbital rotation that diagonalizes the bath (non-interacting system)\n", " hopping_matrix = np.zeros((n_bath, n_bath))\n", " np.fill_diagonal(hopping_matrix[:, 1:], -1)\n", " np.fill_diagonal(hopping_matrix[1:, :], -1)\n", " _, vecs = np.linalg.eigh(hopping_matrix)\n", "\n", " # Expand to include impurity\n", " orbital_rotation = np.zeros((norb, norb))\n", " # Impurity is on the first site\n", " orbital_rotation[0, 0] = 1\n", " orbital_rotation[1:, 1:] = vecs\n", "\n", " # Move the impurity to the center\n", " new_index = n_bath // 2\n", " perm = np.r_[1 : (new_index + 1), 0, (new_index + 1) : norb]\n", " orbital_rotation = orbital_rotation[:, perm]\n", "\n", " return orbital_rotation\n", "\n", "\n", "def rotated(\n", " h1e: np.ndarray, h2e: np.ndarray, orbital_rotation: np.ndarray\n", ") -> tuple[np.ndarray, np.ndarray]:\n", " \"\"\"Rotate the orbital basis of a Hamiltonian.\"\"\"\n", " h1e_rotated = np.einsum(\n", " \"ab,Aa,Bb->AB\",\n", " h1e,\n", " orbital_rotation,\n", " orbital_rotation.conj(),\n", " optimize=\"greedy\",\n", " )\n", " h2e_rotated = np.einsum(\n", " \"abcd,Aa,Bb,Cc,Dd->ABCD\",\n", " h2e,\n", " orbital_rotation,\n", " orbital_rotation.conj(),\n", " orbital_rotation,\n", " orbital_rotation.conj(),\n", " optimize=\"greedy\",\n", " )\n", " return h1e_rotated, h2e_rotated\n", "\n", "\n", "# Total number of spatial orbitals, including the bath sites and the impurity\n", "# This should be an even number\n", "norb = 8\n", "\n", "# System is half-filled\n", "nelec = (norb // 2, norb // 2)\n", "# One orbital is the impurity, the rest are bath sites\n", "n_bath = norb - 1\n", "\n", "# Hamiltonian parameters\n", "hybridization = 1.0\n", "hopping = 1.0\n", "onsite = 10.0\n", "chemical_potential = -0.5 * onsite\n", "\n", "# Generate Hamiltonian in position basis\n", "h1e, h2e = siam_hamiltonian(\n", " norb=norb,\n", " hopping=hopping,\n", " onsite=onsite,\n", " hybridization=hybridization,\n", " chemical_potential=chemical_potential,\n", ")\n", "\n", "# Rotate to momentum basis\n", "orbital_rotation = momentum_basis(norb)\n", "h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())\n", "# In the momentum basis, the impurity is placed in the center\n", "impurity_index = n_bath // 2\n", "\n", "# Use PySCF to compute the exact ground state energy\n", "reference_energy, _ = pyscf.fci.direct_spin1.kernel(h1e, h2e, norb, nelec)" ] }, { "cell_type": "markdown", "id": "e888edf2-7865-41a8-be57-6ddb72dd0cc7", "metadata": {}, "source": [ "A continuación, generamos los circuitos para producir los estados base de Krylov.\n", "Para cada especie de espín, el estado inicial $\\ket{\\psi_0}$ viene dado por la superposición de todas las posibles excitaciones de los tres electrones más cercanos al nivel de Fermi en los 4 modos vacíos más cercanos partiendo del estado $|00\\cdots 0011 \\cdots 11\\rangle$, y realizado mediante la aplicación de siete [XXPlusYYGates](/docs/api/qiskit/qiskit.circuit.library.XXPlusYYGate).\n", "Los estados evolucionados en el tiempo se producen mediante aplicaciones sucesivas de un escalón de Trotter de segundo orden.\n", "\n", "Para una descripción más detallada de este modelo y de cómo se diseñan los circuitos, consulte [\"Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization\"](https://arxiv.org/abs/2501.09702).\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "0f729f86-1814-4d5a-ae65-3f9614e103b3", "metadata": {}, "outputs": [], "source": [ "from typing import Sequence\n", "\n", "import ffsim\n", "import scipy\n", "from qiskit import QuantumCircuit, QuantumRegister\n", "from qiskit.circuit import CircuitInstruction, Qubit\n", "from qiskit.circuit.library import CPhaseGate, XGate, XXPlusYYGate\n", "\n", "\n", "def prepare_initial_state(qubits: Sequence[Qubit], norb: int, nocc: int):\n", " \"\"\"Prepare initial state.\"\"\"\n", " assert norb >= 8\n", " x_gate = XGate()\n", " rot = XXPlusYYGate(0.5 * np.pi, -0.5 * np.pi)\n", " for i in range(nocc):\n", " yield CircuitInstruction(x_gate, [qubits[i]])\n", " yield CircuitInstruction(x_gate, [qubits[norb + i]])\n", " for i in range(3):\n", " for j in range(nocc - i - 1, nocc + i, 2):\n", " yield CircuitInstruction(rot, [qubits[j], qubits[j + 1]])\n", " yield CircuitInstruction(\n", " rot, [qubits[norb + j], qubits[norb + j + 1]]\n", " )\n", " yield CircuitInstruction(rot, [qubits[j + 1], qubits[j + 2]])\n", " yield CircuitInstruction(\n", " rot, [qubits[norb + j + 1], qubits[norb + j + 2]]\n", " )\n", "\n", "\n", "def trotter_step(\n", " qubits: Sequence[Qubit],\n", " time_step: float,\n", " one_body_evolution: np.ndarray,\n", " h2e: np.ndarray,\n", " impurity_index: int,\n", " norb: int,\n", "):\n", " \"\"\"A Trotter step.\"\"\"\n", " # Assume the two-body interaction is just the on-site interaction of the impurity\n", " onsite = h2e[\n", " impurity_index, impurity_index, impurity_index, impurity_index\n", " ]\n", " # Two-body evolution for half the time\n", " yield CircuitInstruction(\n", " CPhaseGate(-0.5 * time_step * onsite),\n", " [qubits[impurity_index], qubits[norb + impurity_index]],\n", " )\n", " # One-body evolution for the full time\n", " yield CircuitInstruction(\n", " ffsim.qiskit.OrbitalRotationJW(norb, one_body_evolution), qubits\n", " )\n", " # Two-body evolution for half the time\n", " yield CircuitInstruction(\n", " CPhaseGate(-0.5 * time_step * onsite),\n", " [qubits[impurity_index], qubits[norb + impurity_index]],\n", " )\n", "\n", "\n", "# Time step\n", "time_step = 0.2\n", "# Number of Krylov basis states\n", "krylov_dim = 8\n", "\n", "# Initialize circuit\n", "qubits = QuantumRegister(2 * norb, name=\"q\")\n", "circuit = QuantumCircuit(qubits)\n", "\n", "# Generate initial state\n", "for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):\n", " circuit.append(instruction)\n", "circuit.measure_all()\n", "\n", "# Create list of circuits, starting with the initial state circuit\n", "circuits = [circuit.copy()]\n", "\n", "# Add time evolution circuits to the list\n", "one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)\n", "for i in range(krylov_dim - 1):\n", " # Remove measurements\n", " circuit.remove_final_measurements()\n", " # Append another Trotter step\n", " for instruction in trotter_step(\n", " qubits,\n", " time_step,\n", " one_body_evolution,\n", " h2e_momentum,\n", " impurity_index,\n", " norb,\n", " ):\n", " circuit.append(instruction)\n", " # Measure qubits\n", " circuit.measure_all()\n", " # Add a copy of the circuit to the list\n", " circuits.append(circuit.copy())" ] }, { "cell_type": "code", "execution_count": 3, "id": "9f2cc4d4-ecac-457a-bcae-558319668e1f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "circuits[0].draw(\"mpl\", scale=0.4, fold=-1)" ] }, { "cell_type": "code", "execution_count": 4, "id": "827976ec-4815-4707-80b1-e13fb2fef309", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "circuits[-1].draw(\"mpl\", scale=0.4, fold=-1)" ] }, { "cell_type": "markdown", "id": "b8ca6be4-61d9-47be-8099-8712c7ecc774", "metadata": {}, "source": [ "\n", "\n", "### Paso 2: Optimizar el problema para la ejecución cuántica\n", "\n" ] }, { "cell_type": "markdown", "id": "a3304e1b-9c7c-4212-8744-d1c62292eced", "metadata": {}, "source": [ "A continuación, optimizamos el circuito para un hardware específico. Por ahora, crearemos un backend genérico con un número determinado de qubits y un conjunto de puertas en el que los circuitos de evolución temporal se descomponen de forma natural.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "2d2fdbff-1e22-45af-a2eb-c334e4328c59", "metadata": {}, "outputs": [], "source": [ "from qiskit.providers.fake_provider import GenericBackendV2\n", "\n", "backend = GenericBackendV2(\n", " 2 * norb, basis_gates=[\"cp\", \"xx_plus_yy\", \"p\", \"x\"]\n", ")" ] }, { "cell_type": "markdown", "id": "ce3e6a52-7b99-49b6-8294-b005efa59cfc", "metadata": {}, "source": [ "Ahora, utilizamos Qiskit para transpilar los circuitos al backend de destino.\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "c8643533-9fec-40bf-a307-da8839b1e444", "metadata": {}, "outputs": [], "source": [ "from qiskit.transpiler import generate_preset_pass_manager\n", "\n", "pass_manager = generate_preset_pass_manager(\n", " optimization_level=3, backend=backend\n", ")\n", "isa_circuits = pass_manager.run(circuits)" ] }, { "cell_type": "markdown", "id": "6cfd3eea-e2d9-40a5-a449-1d3d790a5f2d", "metadata": {}, "source": [ "\n", "\n", "### Paso 3: Ejecutar con el comando « Qiskit primitives »\n", "\n" ] }, { "cell_type": "markdown", "id": "ad48e3e6-1013-45a1-942c-63766fed5819", "metadata": {}, "source": [ "Una vez optimizados los circuitos para su ejecución en hardware, estamos listos para ejecutarlos en el hardware de destino y recoger muestras para la estimación de la energía del estado fundamental. Después de utilizar la primitiva Sampler para muestrear las cadenas de bits de cada circuito, combinamos todos los resultados en un único diccionario de recuentos y trazamos las 20 cadenas de bits más muestreadas.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "80eee553-60d6-4258-88ab-d8d120418c36", "metadata": {}, "outputs": [], "source": [ "from qiskit.visualization import plot_histogram\n", "from qiskit.primitives import StatevectorSampler\n", "\n", "# Sample from the circuits\n", "sampler = StatevectorSampler()\n", "job = sampler.run(isa_circuits, shots=500)" ] }, { "cell_type": "code", "execution_count": 8, "id": "10af4663-7375-4b50-bae6-9f3d5106457b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from qiskit.primitives import BitArray\n", "\n", "# Combine the shots from the individual Trotter circuits\n", "bit_array = BitArray.concatenate_shots(\n", " [result.data.meas for result in job.result()]\n", ")\n", "\n", "plot_histogram(bit_array.get_counts(), number_to_keep=20)" ] }, { "cell_type": "markdown", "id": "2aa74455-d16b-4ac3-a354-54a79d5c5759", "metadata": {}, "source": [ "\n", "\n", "### Paso 4: Posprocesamiento y devolución del resultado al formato clásico deseado\n", "\n" ] }, { "cell_type": "markdown", "id": "d3f7713c-a4e6-407f-94de-0e5cfeb1134c", "metadata": {}, "source": [ "Ahora, ejecutamos el algoritmo SQD utilizando la función `diagonalize_fermionic_hamiltonian` . Consulte la [documentación de la API](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian) para obtener explicaciones sobre los argumentos de esta función.\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "7609d1e1-e8ef-48e1-a965-97927f403163", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration 1\n", "\tSubsample 0\n", "\t\tEnergy: -13.4222953188441\n", "\t\tSubspace dimension: 529\n", "\tSubsample 1\n", "\t\tEnergy: -13.42237556285828\n", "\t\tSubspace dimension: 784\n", "\tSubsample 2\n", "\t\tEnergy: -13.422045397387413\n", "\t\tSubspace dimension: 529\n", "Iteration 2\n", "\tSubsample 0\n", "\t\tEnergy: -13.422379583305478\n", "\t\tSubspace dimension: 900\n", "\tSubsample 1\n", "\t\tEnergy: -13.422376197704326\n", "\t\tSubspace dimension: 841\n", "\tSubsample 2\n", "\t\tEnergy: -13.422421162849295\n", "\t\tSubspace dimension: 1089\n", "Iteration 3\n", "\tSubsample 0\n", "\t\tEnergy: -13.422421164670345\n", "\t\tSubspace dimension: 1156\n", "\tSubsample 1\n", "\t\tEnergy: -13.422421492737689\n", "\t\tSubspace dimension: 1156\n", "\tSubsample 2\n", "\t\tEnergy: -13.422421205869572\n", "\t\tSubspace dimension: 1156\n", "Iteration 4\n", "\tSubsample 0\n", "\t\tEnergy: -13.422421494558726\n", "\t\tSubspace dimension: 1225\n", "\tSubsample 1\n", "\t\tEnergy: -13.422421492737689\n", "\t\tSubspace dimension: 1156\n", "\tSubsample 2\n", "\t\tEnergy: -13.422421492737689\n", "\t\tSubspace dimension: 1156\n" ] } ], "source": [ "from qiskit_addon_sqd.fermion import (\n", " SCIResult,\n", " diagonalize_fermionic_hamiltonian,\n", ")\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}\")\n", " print(\n", " f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n", " )\n", "\n", "\n", "rng = np.random.default_rng(24)\n", "result = diagonalize_fermionic_hamiltonian(\n", " h1e_momentum,\n", " h2e_momentum,\n", " bit_array,\n", " samples_per_batch=100,\n", " norb=norb,\n", " nelec=nelec,\n", " num_batches=3,\n", " max_iterations=5,\n", " symmetrize_spin=True,\n", " callback=callback,\n", " seed=rng,\n", ")" ] }, { "cell_type": "markdown", "id": "3dee9c61-fc42-48e9-8888-af8fc831cd5c", "metadata": {}, "source": [ "La siguiente celda de código muestra los resultados. El primer gráfico muestra la energía calculada en función del número de iteraciones de recuperación de la configuración, y el segundo gráfico muestra la ocupación media de cada orbital espacial tras la iteración final. Dado que se trata de un problema tan sencillo, la primera iteración ya nos acerca mucho a la energía exacta (fíjate en la escala del eje y).\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "b6879566-8bf5-4c28-bfb6-b2686692e3d3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Reference energy: -13.42249\n", "SQD energy: -13.42242\n", "Absolute error: 0.00007\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "min_es = [\n", " min(result, key=lambda res: res.energy).energy\n", " for result in result_history\n", "]\n", "min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])\n", "\n", "# Data for energies plot\n", "x1 = range(len(result_history))\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, min_es, label=\"energy\", marker=\"o\")\n", "axs[0].set_xticks(x1)\n", "axs[0].set_xticklabels(x1)\n", "axs[0].axhline(\n", " y=reference_energy,\n", " color=\"#BF5700\",\n", " linestyle=\"--\",\n", " label=\"reference energy\",\n", ")\n", "axs[0].set_title(\"Approximated Ground State Energy vs SQD Iterations\")\n", "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n", "axs[0].set_ylabel(\"Energy\", 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", "print(f\"Reference energy: {reference_energy:.5f}\")\n", "print(f\"SQD energy: {min_e:.5f}\")\n", "print(f\"Absolute error: {abs(min_e - reference_energy):.5f}\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "0c9f4976-d770-426a-822e-e6756c2cfbe7", "metadata": {}, "source": [ "\n", "\n", "### Comprueba el consumo energético\n", "\n", "Se garantiza que la energía devuelta por el SQD constituye un límite superior de la energía real del estado fundamental. El valor de la energía puede verificarse, ya que SQD también devuelve los coeficientes del vector de estado que aproxima el estado fundamental. Se puede calcular la energía a partir del vector de estado utilizando sus matrices de densidad reducidas de una y dos partículas, tal y como se muestra en la siguiente celda de código.\n", "\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "e2b9de72-61cf-49d3-a1b5-f043e4b16956", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Recomputed energy: -13.42242\n" ] } ], "source": [ "rdm1 = result.sci_state.rdm(rank=1, spin_summed=True)\n", "rdm2 = result.sci_state.rdm(rank=2, spin_summed=True)\n", "\n", "energy = np.sum(h1e_momentum * rdm1) + 0.5 * np.sum(h2e_momentum * rdm2)\n", "\n", "print(f\"Recomputed energy: {energy:.5f}\")" ] }, { "cell_type": "markdown", "id": "5221928a-79ff-4f54-90b4-1fe4d7739aae", "metadata": {}, "source": [ "\n", "\n", "## Ejemplo de hardware a gran escala\n", "\n", "Ahora ejecutamos un ejemplo más extenso en una QPU real.\n", "Para la energía de referencia, utilizamos los resultados de un cálculo [DMRG](https://en.wikipedia.org/wiki/Density_matrix_renormalization_group) que se realizó por separado.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "933037d8-847e-4986-80da-5ac8d677b2ff", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using backend ibm_boston\n", "Iteration 1\n", "\tSubsample 0\n", "\t\tEnergy: -28.63965951544449\n", "\t\tSubspace dimension: 9801\n", "\tSubsample 1\n", "\t\tEnergy: -28.625588929202006\n", "\t\tSubspace dimension: 9409\n", "\tSubsample 2\n", "\t\tEnergy: -28.647371834135498\n", "\t\tSubspace dimension: 8281\n", "Iteration 2\n", "\tSubsample 0\n", "\t\tEnergy: -28.67213260849567\n", "\t\tSubspace dimension: 29584\n", "\tSubsample 1\n", "\t\tEnergy: -28.670340686158816\n", "\t\tSubspace dimension: 27225\n", "\tSubsample 2\n", "\t\tEnergy: -28.669976379525988\n", "\t\tSubspace dimension: 31329\n", "Iteration 3\n", "\tSubsample 0\n", "\t\tEnergy: -28.68622875601382\n", "\t\tSubspace dimension: 36100\n", "\tSubsample 1\n", "\t\tEnergy: -28.698569623143126\n", "\t\tSubspace dimension: 34225\n", "\tSubsample 2\n", "\t\tEnergy: -28.694848533971882\n", "\t\tSubspace dimension: 33856\n", "Iteration 4\n", "\tSubsample 0\n", "\t\tEnergy: -28.69883392844593\n", "\t\tSubspace dimension: 42025\n", "\tSubsample 1\n", "\t\tEnergy: -28.701289495200996\n", "\t\tSubspace dimension: 38025\n", "\tSubsample 2\n", "\t\tEnergy: -28.699319594978245\n", "\t\tSubspace dimension: 45369\n", "Iteration 5\n", "\tSubsample 0\n", "\t\tEnergy: -28.701936886834154\n", "\t\tSubspace dimension: 51076\n", "\tSubsample 1\n", "\t\tEnergy: -28.702468711812013\n", "\t\tSubspace dimension: 53824\n", "\tSubsample 2\n", "\t\tEnergy: -28.702298147575938\n", "\t\tSubspace dimension: 52900\n", "Reference energy: -28.70660\n", "SQD energy: -28.70247\n", "Absolute error: 0.00413\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from qiskit_ibm_runtime import SamplerV2 as Sampler\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "\n", "# Model parameters\n", "norb = 20\n", "nelec = (norb // 2, norb // 2)\n", "n_bath = norb - 1\n", "hybridization = 1.0\n", "hopping = 1.0\n", "onsite = 10.0\n", "chemical_potential = -0.5 * onsite\n", "\n", "# Generate Hamiltonian and orbital rotation\n", "h1e, h2e = siam_hamiltonian(\n", " norb=norb,\n", " hopping=hopping,\n", " onsite=onsite,\n", " hybridization=hybridization,\n", " chemical_potential=chemical_potential,\n", ")\n", "orbital_rotation = momentum_basis(norb)\n", "h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())\n", "impurity_index = n_bath // 2\n", "\n", "# Set reference energy to DMRG value computed separately\n", "reference_energy = -28.70659686\n", "\n", "# Algorithm parameters\n", "time_step = 0.2\n", "krylov_dim = 8\n", "\n", "# Construct circuits\n", "qubits = QuantumRegister(2 * norb, name=\"q\")\n", "circuit = QuantumCircuit(qubits)\n", "for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):\n", " circuit.append(instruction)\n", "circuit.measure_all()\n", "circuits = [circuit.copy()]\n", "one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)\n", "for i in range(krylov_dim - 1):\n", " circuit.remove_final_measurements()\n", " for instruction in trotter_step(\n", " qubits,\n", " time_step,\n", " one_body_evolution,\n", " h2e_momentum,\n", " impurity_index,\n", " norb,\n", " ):\n", " circuit.append(instruction)\n", " circuit.measure_all()\n", " circuits.append(circuit.copy())\n", "\n", "# Initialize hardware backend\n", "service = QiskitRuntimeService()\n", "backend = service.least_busy(\n", " operational=True, simulator=False, min_num_qubits=127\n", ")\n", "print(f\"Using backend {backend.name}\")\n", "\n", "# Transpile to backend\n", "pass_manager = generate_preset_pass_manager(\n", " optimization_level=3, backend=backend\n", ")\n", "isa_circuits = pass_manager.run(circuits)\n", "\n", "# Sample from the circuits\n", "sampler = Sampler(backend)\n", "sampler.options.environment.job_tags = [\"TUT_SKQD\"]\n", "job = sampler.run(isa_circuits, shots=500)\n", "\n", "# Combine the shots from the individual Trotter circuits\n", "bit_array = BitArray.concatenate_shots(\n", " [result.data.meas for result in job.result()]\n", ")\n", "\n", "# Run configuration recovery and diagonalization\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}\")\n", " print(\n", " f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n", " )\n", "\n", "\n", "rng = np.random.default_rng(24)\n", "result = diagonalize_fermionic_hamiltonian(\n", " h1e_momentum,\n", " h2e_momentum,\n", " bit_array,\n", " samples_per_batch=100,\n", " norb=norb,\n", " nelec=nelec,\n", " num_batches=3,\n", " max_iterations=5,\n", " symmetrize_spin=True,\n", " callback=callback,\n", " seed=rng,\n", ")\n", "\n", "\n", "# Plot results\n", "min_es = [\n", " min(result, key=lambda res: res.energy).energy\n", " for result in result_history\n", "]\n", "min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])\n", "x1 = range(len(result_history))\n", "y2 = np.sum(result.orbital_occupancies, axis=0)\n", "x2 = range(len(y2))\n", "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n", "axs[0].plot(x1, min_es, label=\"energy\", marker=\"o\")\n", "axs[0].set_xticks(x1)\n", "axs[0].set_xticklabels(x1)\n", "axs[0].axhline(\n", " y=reference_energy,\n", " color=\"#BF5700\",\n", " linestyle=\"--\",\n", " label=\"reference energy\",\n", ")\n", "axs[0].set_title(\"Approximated Ground State Energy vs SQD Iterations\")\n", "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n", "axs[0].set_ylabel(\"Energy\", fontdict={\"fontsize\": 12})\n", "axs[0].legend()\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", "print(f\"Reference energy: {reference_energy:.5f}\")\n", "print(f\"SQD energy: {min_e:.5f}\")\n", "print(f\"Absolute error: {abs(min_e - reference_energy):.5f}\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "482ebea3-84b8-471b-bddc-b282e23162ad", "metadata": {}, "source": [ "\n", "\n", "## Próximos pasos\n", "\n", "\n", " Si te ha parecido interesante este trabajo, quizá te interese el siguiente material:\n", "\n", " * [Diagonalización cuántica basada en muestras de un hamiltoniano químico](/docs/tutorials/sample-based-quantum-diagonalization) : un tutorial relacionado que utiliza un enfoque variacional heurístico en lugar de circuitos de Trotter\n", " * [Diagonalización cuántica de Krylov de hamiltonianos de red](/docs/tutorials/krylov-quantum-diagonalization) : un tutorial sobre el método KQD\n", " * [Documentación de la API del complemento SQD](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian) : referencia de la `diagonalize_fermionic_hamiltonian` función\n", " * [*«Algoritmo centrado en el quantum para la diagonalización de Krylov basada en muestras*](https://arxiv.org/abs/2501.09702) »: el artículo en el que se basa este tutorial\n", "\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": 9 }, "nbformat": 4, "nbformat_minor": 4 }