{
"cells": [
{
"cell_type": "markdown",
"id": "8e82ead7",
"metadata": {},
"source": [
"---\n",
"title: \"Simulación de la dispersión de neutrones en materiales cuánticos mediante circuitos cuánticos\"\n",
"description: \"Calcular el factor de estructura dinámica de un imán cuántico utilizando circuitos de Trotter y simulación MPS.\"\n",
"---\n",
"\n",
"{/* cspell:ignore DSFs spinon DMRG viridis fontsize vmax vmin cadetblue antiferromagnet COBYQA */}\n",
"\n",
"\n",
"\n",
"# Simulación de la dispersión de neutrones en materiales cuánticos mediante circuitos cuánticos\n",
"\n",
"*Estimación de tiempo de ejecución: 13 minutos en un procesador Heron r2 (NOTA: Se trata únicamente de una estimación. (El tiempo de ejecución puede variar.)*\n",
"\n",
"\n",
" Utiliza este tutorial para aprender la implementación paso a paso, con el cálculo clásico ejecutándose localmente en tu ordenador portátil y la ejecución en hardware en una QPU de IBM Quantum®. El ejemplo a gran escala requiere una cantidad considerable de memoria, y la compilación cuántica aproximada (AQC) puede tardar varias horas en un ordenador portátil estándar. Para delegar la compresión a los recursos de computación y memori Qiskit Serverless, utiliza el [tutorial](/docs/tutorials/simulate-neutron-scattering-with-a-serverless-workflow) complementario sobre «Serverless».\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "11033666",
"metadata": {},
"source": [
"\n",
"\n",
"## Resultados del aprendizaje\n",
"\n",
"* Cómo se relacionan los espectros de dispersión inelástica de neutrones (INS) con los factores de estructura dinámicos (DSF) de los modelos cuánticos de espín.\n",
"* Cómo preparar un estado fundamental, aplicar una perturbación local y llevar a cabo la evolución temporal de Trotter en un circuito cuántico.\n",
"* Cómo utilizar la compilación cuántica aproximada (AQC) con `qiskit-addon-aqc-tensor` para comprimir circuitos de Trotter profundos para su ejecución en hardware.\n",
"* Cómo extraer la función de Green retardada (RGF) a partir de los valores esperados de los qubits y aplicarles la transformada de Fourier para obtener una DSF.\n",
"\n",
"\n",
"\n",
"## Requisitos previos\n",
"\n",
"* [Fundamentos de la información cuántica](/learning/courses/basics-of-quantum-information)\n",
"* [Diseño de algoritmos variacionales](/learning/courses/variational-algorithm-design)\n",
"* [Introducción a « Qiskit primitives » (Estimador y muestreador)](/docs/guides/qiskit-runtime-primitives)\n",
"\n",
"\n",
"\n",
"## En segundo plano\n",
"\n",
"En este tutorial, reproducimos los resultados de [Lee et al., arXiv:2603.15608](https://arxiv.org/abs/2603.15608).\n",
"\n",
"La dispersión de neutrones ayuda a los investigadores a caracterizar materiales relevantes para la sostenibilidad, incluidos los electrodos de las baterías. Este tutorial analiza cómo los circuitos cuánticos pueden simular excitaciones magnéticas y relacionar los modelos teóricos con las mediciones de dispersión. Su ejemplo de imán cuántico desarrolla métodos para estudiar el comportamiento de los materiales, lo que contribuye a ampliar el conjunto de herramientas de investigación para el descubrimiento de nuevos materiales.\n",
"\n",
"\n",
"\n",
"### La dispersión inelástica de neutrones y el factor de estructura dinámico\n",
"\n",
"La dispersión inelástica de neutrones (INS) es una de las técnicas experimentales más potentes para estudiar las excitaciones magnéticas en los materiales cuánticos. Cuando un haz de neutrones térmicos o fríos incide sobre un cristal, los neutrones individuales intercambian tanto momento $\\mathbf{q}$ como energía $\\omega$ con el subsistema magnético. La intensidad de dispersión medida es proporcional al factor de estructura dinámica (DSF),\n",
"\n",
"$S^{\\alpha\\beta}(q,\\omega) = \\sum_{j} e^{-iq\\,j}\\int_{-\\infty}^{\\infty} dt\\; e^{i\\omega t}\\, \\langle S_0^{\\alpha}(0)\\, S_j^{\\beta}(t)\\rangle,$\n",
"\n",
"que codifica todas las correlaciones espacio-temporales de los grados de libertad de espín.\n",
"\n",
"\n",
"\n",
"### KCuF$_3$: un imán canónico de tipo líquido de Luttinger\n",
"\n",
"El fluoruro de potasio y cobre ( KCuF$_3$ ) es un antiferromagnético cuasiunidimensional en el que las cadenas de iones de espín- $\\frac{1}{2}$ Cu $^{2+}$ interactúan a través de un intercambio de Heisenberg entre vecinos más cercanos $J$, mientras que el acoplamiento entre cadenas es solo $\\sim 2.7\\%$ de $J$. En $T = 6\\;\\mathrm{K}$, donde se dispone de datos de espectroscopia de resonancia de espín in situ (INS), el espectro está dominado por excitaciones de espínones fraccionadas, características de un líquido de Tomonaga-Luttinger. Dado que la dinámica intracadena viene bien descrita por el hamiltoniano unidimensional de tipo «spin- $\\frac{1}{2}$ XXZ» en el punto isotrópico ( $\\epsilon = 1$ ),\n",
"\n",
"$H = J\\sum_{i}\\left[S_i^Z S_{i+1}^Z + \\epsilon\\left(S_i^X S_{i+1}^X + S_i^Y S_{i+1}^Y\\right)\\right],$\n",
"\n",
"KCuF$_3$ constituye un punto de referencia ideal para la simulación cuántica: el hamiltoniano es lo suficientemente sencillo como para implementarlo en un procesador cuántico, pero el estado fundamental presenta un fuerte entrelazamiento y el espectro de excitación presenta un amplio continuo de dos espinones.\n",
"\n",
"Nota: este tutorial establece $J = 1$ como unidad de energía y adopta la normalización $H = J\\sum_i[\\ldots]$, que se corresponde con el hamiltoniano local $H_\\mathrm{loc} = J(S^X S^X + S^Y S^Y + S^Z S^Z)$ utilizado en la implementación del circuito del artículo (fig. S3 (del anexo). El hamiltoniano completo del artículo (ecuación 3) conlleva un factor global adicional de 2, por lo que el « $J$ » del artículo es el doble del « $J$ » utilizado aquí.\n",
"\n",
"\n",
"\n",
"### Lo que simulamos y medimos\n",
"\n",
"La magnitud física que calculamos es la **función de Green retardada** (RGF), definida como la función de correlación espín-espín dependiente del tiempo\n",
"\n",
"$G^R_{\\alpha,\\beta}(j, j_c, t) = -\\frac{i}{2}\\,\\langle\\psi_{\\mathrm{GS}}|\\,S_j^\\alpha(t)\\,S_{j_c}^\\beta(0) - S_{j_c}^\\beta(0)\\,S_j^\\alpha(t)\\,|\\psi_{\\mathrm{GS}}\\rangle,$\n",
"\n",
"donde $j_c$ es un sitio de referencia (el centro de la cadena) y $S_j^\\alpha(t) = e^{iHt}S_j^\\alpha e^{-iHt}$ es el operador de espín según la interpretación de Heisenberg. En este tutorial nos centramos en el componente « $zz$ » ( $\\alpha = \\beta = z$ ). En un ordenador cuántico, se accede al RGF preparando el estado fundamental, aplicando una perturbación local en $j_c$, haciendo evolucionar en el tiempo el estado perturbado y midiendo el valor esperado de un solo qubit $\\langle\\sigma_j^z\\rangle$ en cada sitio $j$ para cada paso temporal. La idea clave es que cada $\\langle\\sigma_j^z\\rangle$ representa la diferencia con respecto a la magnetización del estado fundamental; dado que el antiferromagneto isotrópico de Heisenberg tiene una magnetización neta por sitio igual a cero ( $\\langle\\sigma_j^z\\rangle_\\mathrm{GS} = 0$ ), el valor bruto medido da directamente como resultado $G^R(j, j_c, t)$ sin necesidad de realizar ninguna resta explícita.\n",
"\n",
"Al recopilar $G^R(j, j_c, t)$ en todos los puntos y pasos temporales, obtenemos un conjunto de datos bidimensional al que luego se le aplica una transformada de Fourier tanto en el espacio como en el tiempo para obtener el **factor de estructura dinámico** $S(q,\\omega)$. El DSF es la magnitud que se mide directamente en un experimento INS: nos indica qué excitaciones magnéticas existen en cada momento $q$ y energía $\\omega$. Para la cadena de Heisenberg isotrópica, el espectro de excitación exacto es un **continuo de dos espinones**, una banda ancha de intensidad de dispersión cuya forma sirve como un riguroso punto de referencia de extremo a extremo para la simulación cuántica: valida a la vez la preparación del estado fundamental, la perturbación, la evolución temporal de Trotter y el protocolo de medición.\n",
"\n",
"\n",
"\n",
"### Flujo de trabajo de simulación cuántica\n",
"\n",
"El flujo de trabajo reproduce la física de un evento del INS. (1) Preparamos el estado fundamental de muchos cuerpos $|\\psi_{\\mathrm{GS}}\\rangle$ en $n$ qubits; (2) aplicamos una perturbación local de inversión de espín $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$ en el centro de la cadena para simular la transferencia de espín del neutrón; (3) hacemos evolucionar el sistema bajo $H$ durante pasos de tiempo discretos utilizando la trotterización de segundo orden; y (4) medimos $\\langle\\sigma_i^z\\rangle$ en cada qubit en cada paso para obtener la RGF. A continuación, una transformada de Fourier discreta bidimensional da como resultado el DSF $S(q,\\omega)$.\n",
"\n",
"La magnitud observable que medimos en cada paso temporal es $\\sigma_i^z$ en cada qubit $i$. En Qiskit, esto se representa como una lista de `SparsePauliOp` operadores: un operador de un solo qubit $Z$ integrado en la cadena de identidad de $n$ -qubit para cada sitio. Estas variables observables se construyen una vez por cada instancia del problema durante la fase de mapeo del problema (paso 1) y se reutilizan para todos los circuitos a esa escala.\n",
"\n",
"\n",
"\n",
"### Compilación cuántica aproximada (AQC)\n",
"\n",
"Los circuitos Deep Trotter pueden comprimirse mediante **la compilación cuántica aproximada (AQC)**, que sustituye las primeras capas de Trotter por un ansatz parametrizado más corto, cuyos parámetros se optimizan de forma clásica para maximizar la fidelidad a nivel de MPS con respecto al circuito profundo original. Los pasos restantes de Trotter se añaden tal cual, lo que da lugar a un circuito «mixto» de AQC + Trotter con un número considerablemente menor de puertas de dos qubits.\n",
"\n",
"\n",
"\n",
"### Simulación MPS\n",
"\n",
"En el caso de un sistema unidimensional, los métodos de estado de producto matricial (MPS) permiten simular de forma eficaz tanto la preparación del estado fundamental (mediante el grupo de renormalización de la matriz de densidad, o DMRG) como la evolución temporal a nivel de circuito. Al controlar la dimensión del enlace $\\chi$, se establece un equilibrio entre la precisión y el coste computacional. En este tutorial utilizamos la simulación MPS para `qiskit-addon-aqc-tensor` calcular ansätze de AQC de alta fidelidad que comprimen los circuitos de Trotter profundos para su ejecución en hardware.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "de1ad5a2",
"metadata": {},
"source": [
"\n",
"\n",
"## Requisitos\n",
"\n",
"Antes de empezar este tutorial, asegúrate de tener instalado lo siguiente:\n",
"\n",
"* Qiskit SDK con funciones [de visualización](/docs/api/qiskit/visualization)\n",
"* Qiskit Runtime (`pip install qiskit-ibm-runtime`)\n",
"* `qiskit-addon-aqc-tensor` con `quimb` y JAX extras (`pip install 'qiskit-addon-aqc-tensor[quimb-jax]'`)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9938e4bd",
"metadata": {},
"source": [
"\n",
"\n",
"## Configuración\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "71835ff6",
"metadata": {},
"outputs": [],
"source": [
"import timeit\n",
"import warnings\n",
"from collections.abc import Iterator, Sequence\n",
"from functools import partial\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import quimb.tensor as qtn\n",
"import scipy.optimize\n",
"from numpy.typing import NDArray\n",
"from qiskit import QuantumCircuit\n",
"from qiskit.circuit import CircuitInstruction, ParameterVector, Qubit\n",
"from qiskit.circuit.library import PauliEvolutionGate\n",
"from qiskit.primitives import StatevectorEstimator\n",
"from qiskit.quantum_info import SparsePauliOp\n",
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
"from qiskit_addon_aqc_tensor import generate_ansatz_from_circuit\n",
"from qiskit_addon_aqc_tensor.objective import MaximizeStateFidelity\n",
"from qiskit_addon_aqc_tensor.simulation import tensornetwork_from_circuit\n",
"from qiskit_addon_aqc_tensor.simulation.quimb import QuimbSimulator\n",
"from qiskit_ibm_runtime import EstimatorV2 as Estimator\n",
"from qiskit_ibm_runtime import QiskitRuntimeService\n",
"from qiskit_quimb import quimb_circuit\n",
"from scipy.sparse import SparseEfficiencyWarning\n",
"\n",
"# scipy.sparse.linalg.expm (reached via PauliEvolutionGate.to_matrix) prefers\n",
"# CSC and warns when handed the CSR matrix from SparsePauliOp.to_matrix; the\n",
"# result is unaffected, so silence the cosmetic warning.\n",
"warnings.filterwarnings(\"ignore\", category=SparseEfficiencyWarning)\n",
"\n",
"\n",
"def xxz_hamiltonian_mpo(\n",
" n_qubits: int, interaction: float = 1.0, anisotropy: float = 1.0\n",
") -> qtn.MatrixProductOperator:\n",
" \"\"\"1D XXZ Hamiltonian as a quimb MPO.\n",
"\n",
" Builds the Hamiltonian using ``qtn.SpinHam1D``.\n",
"\n",
" Args:\n",
" n_qubits: Number of sites.\n",
" interaction: Overall interaction strength (J in the paper).\n",
" anisotropy: XY/Z anisotropy (ε in the paper). ``ε = 1`` is the\n",
" isotropic Heisenberg point; ``ε = 0`` is the Ising limit.\n",
"\n",
" Returns:\n",
" The Hamiltonian as a matrix product operator.\n",
" \"\"\"\n",
" builder = qtn.SpinHam1D(S=1 / 2)\n",
" builder += interaction * anisotropy * 0.5, \"+\", \"-\"\n",
" builder += interaction * anisotropy * 0.5, \"-\", \"+\"\n",
" builder += interaction, \"Z\", \"Z\"\n",
" return builder.build_mpo(L=n_qubits)\n",
"\n",
"\n",
"def build_ground_state_ansatz(n_qubits: int, n_layers: int) -> QuantumCircuit:\n",
" \"\"\"Hamiltonian variational ansatz (HVA) circuit for ground-state preparation.\n",
"\n",
" Starts from a product of singlet pairs and applies alternating\n",
" odd/even layers of parameterized XXZ pair-evolution gates.\n",
"\n",
" The returned circuit is parameterized: it carries a ``ParameterVector``\n",
" named ``\"theta\"`` of length ``2 * n_layers`` whose values must be\n",
" assigned (e.g. via ``circuit.assign_parameters``) before simulation.\n",
" Element ``2 * r`` is the odd-layer angle and element ``2 * r + 1`` is the\n",
" even-layer angle of layer ``r``.\n",
"\n",
" Args:\n",
" n_qubits: Number of qubits (must be even).\n",
" n_layers: Number of HVA layers.\n",
"\n",
" Returns:\n",
" The parameterized HVA preparation circuit.\n",
" \"\"\"\n",
" theta = ParameterVector(\"theta\", 2 * n_layers)\n",
" circuit = QuantumCircuit(n_qubits)\n",
" # Initial singlet product state\n",
" for i in range(n_qubits // 2):\n",
" circuit.x(2 * i)\n",
" circuit.x(2 * i + 1)\n",
" circuit.h(2 * i + 1)\n",
" circuit.cx(2 * i + 1, 2 * i)\n",
" # Variational layers: each pair gate is exp(-i theta (XX+YY+ZZ)/2)\n",
" pair_ham = SparsePauliOp(\n",
" [\"XX\", \"YY\", \"ZZ\"], coeffs=[0.5, 0.5, 0.5]\n",
" ) # H_pair (HVA form)\n",
" for r in range(n_layers):\n",
" for i in range(1, (n_qubits + 1) // 2): # odd layer\n",
" circuit.append(\n",
" PauliEvolutionGate(pair_ham, time=theta[2 * r]),\n",
" [2 * i - 1, 2 * i],\n",
" )\n",
" for i in range(n_qubits // 2): # even layer\n",
" circuit.append(\n",
" PauliEvolutionGate(pair_ham, time=theta[2 * r + 1]),\n",
" [2 * i, 2 * i + 1],\n",
" )\n",
" return circuit\n",
"\n",
"\n",
"def optimize_ground_state_ansatz(\n",
" ansatz: QuantumCircuit,\n",
" x0: NDArray[np.floating],\n",
" target_mps: qtn.MatrixProductState,\n",
" *,\n",
" max_bond: int | None = None,\n",
" cutoff: float = 1e-10,\n",
" method: str = \"COBYQA\",\n",
" options: dict | None = None,\n",
") -> scipy.optimize.OptimizeResult:\n",
" \"\"\"Optimize HVA parameters by maximizing fidelity with a target MPS.\n",
"\n",
" The HVA circuit is simulated as a matrix product state with the given\n",
" bond-dimension truncation, and the parameters are optimized to maximize\n",
" the state fidelity ``||**2`` with the DMRG ground\n",
" state. Both states are normalized, so the minimized objective is the\n",
" infidelity ``1 - ||**2``.\n",
"\n",
" Args:\n",
" ansatz: Parameterized HVA circuit from ``build_ground_state_ansatz``.\n",
" The length of ``x0`` must equal ``ansatz.num_parameters``.\n",
" x0: Initial parameters.\n",
" target_mps: Target MPS (DMRG ground state) to maximize fidelity with.\n",
" max_bond: Maximum MPS bond dimension during gate application.\n",
" cutoff: Singular-value cutoff during gate application.\n",
" method: ``scipy.optimize.minimize`` method.\n",
" options: Options dict forwarded to ``scipy.optimize.minimize``.\n",
"\n",
" Returns:\n",
" The Scipy OptimizeResult.\n",
" \"\"\"\n",
"\n",
" def infidelity(params: NDArray[np.floating]) -> float:\n",
" circuit = ansatz.assign_parameters(params)\n",
" circuit_mps = quimb_circuit(\n",
" circuit.decompose([\"PauliEvolution\"]),\n",
" quimb_circuit_class=qtn.CircuitMPS,\n",
" max_bond=max_bond,\n",
" cutoff=cutoff,\n",
" )\n",
" return 1 - abs(circuit_mps.psi.H @ target_mps) ** 2\n",
"\n",
" return scipy.optimize.minimize(\n",
" infidelity, np.asarray(x0), method=method, options=options\n",
" )\n",
"\n",
"\n",
"def trotter_evolution(\n",
" qubits: Sequence[Qubit],\n",
" interaction: float,\n",
" anisotropy: float,\n",
" time_step: float,\n",
" n_steps: int,\n",
") -> Iterator[CircuitInstruction]:\n",
" \"\"\"Second-order Trotter steps of the XXZ pair Hamiltonian.\n",
"\n",
" While the paper used a hand-optimized circuit for the Trotter steps, we use\n",
" PauliEvolutionGate here for simplicity and generality. The final two-qubit gate\n",
" count and gate depth are equivalent when transpiled with ``optimization_level=3``.\n",
"\n",
" Args:\n",
" qubits: Qubits to act on (length ``n_qubits``).\n",
" interaction: Overall interaction strength (J in the paper).\n",
" anisotropy: XY/Z anisotropy (ε in the paper). ``ε = 1`` is the\n",
" isotropic Heisenberg point; ``ε = 0`` is the Ising limit.\n",
" time_step: Per-step Trotter time.\n",
" n_steps: Number of Trotter steps.\n",
"\n",
" Yields:\n",
" ``CircuitInstruction``s implementing the Trotter steps.\n",
" \"\"\"\n",
" if n_steps == 0:\n",
" return\n",
" n_qubits = len(qubits)\n",
" pair_ham = SparsePauliOp(\n",
" [\"XX\", \"YY\", \"ZZ\"],\n",
" coeffs=[\n",
" 0.25 * interaction * anisotropy,\n",
" 0.25 * interaction * anisotropy,\n",
" 0.25 * interaction,\n",
" ],\n",
" )\n",
" half_evo = PauliEvolutionGate(pair_ham, time=time_step / 2)\n",
" full_evo = PauliEvolutionGate(pair_ham, time=time_step)\n",
" for i in range(n_qubits // 2): # half even layer\n",
" yield CircuitInstruction(half_evo, (qubits[2 * i], qubits[2 * i + 1]))\n",
" for i in range(n_qubits // 2 - 1): # full odd layer\n",
" yield CircuitInstruction(\n",
" full_evo, (qubits[2 * i + 1], qubits[2 * i + 2])\n",
" )\n",
" for _ in range(n_steps - 1): # interior steps\n",
" for i in range(n_qubits // 2):\n",
" yield CircuitInstruction(\n",
" full_evo, (qubits[2 * i], qubits[2 * i + 1])\n",
" )\n",
" for i in range(n_qubits // 2 - 1):\n",
" yield CircuitInstruction(\n",
" full_evo, (qubits[2 * i + 1], qubits[2 * i + 2])\n",
" )\n",
" for i in range(n_qubits // 2): # half even layer\n",
" yield CircuitInstruction(half_evo, (qubits[2 * i], qubits[2 * i + 1]))\n",
"\n",
"\n",
"def get_dsf(\n",
" n_qubits: int,\n",
" rgf_mat: NDArray[np.floating],\n",
" time_step: float,\n",
" n_steps: int,\n",
" n_points_momentum: int,\n",
" n_points_frequency: int,\n",
") -> NDArray[np.floating]:\n",
" \"\"\"Compute the dynamical structure factor from the retarded Green's function.\n",
"\n",
" Uses the center-site approximation and a discrete Fourier transform.\n",
" The result is symmetrized about the momentum axis and clipped to\n",
" non-negative values, ready for plotting.\n",
"\n",
" Args:\n",
" n_qubits: Number of qubits (sites).\n",
" rgf_mat: RGF matrix of shape ``(n_steps, n_qubits)``.\n",
" time_step: Trotter time-step size.\n",
" n_steps: Number of time steps.\n",
" n_points_momentum: Number of momentum points.\n",
" n_points_frequency: Number of frequency points.\n",
"\n",
" Returns:\n",
" DSF array of shape ``(n_points_frequency, n_points_momentum)``,\n",
" symmetrized about the momentum axis and clipped to non-negative values.\n",
" \"\"\"\n",
" max_frequency = np.pi / time_step\n",
" momentum_range = np.linspace(0, 2 * np.pi, n_points_momentum)\n",
" frequency_range = np.linspace(0, max_frequency, n_points_frequency)\n",
" result = np.zeros((frequency_range.shape[0], momentum_range.shape[0]))\n",
" center = n_qubits // 2 - 1\n",
" for iw, w in enumerate(frequency_range):\n",
" exponent = np.exp(1j * w * time_step * np.arange(1, n_steps + 1))\n",
" # S = sigma/2, so two spin operators contribute a factor of (1/2)(1/2) = 1/4.\n",
" rgf_omega = (\n",
" np.dot(rgf_mat.T, exponent) * time_step / 4\n",
" ) # S(omega): time Fourier slice of the Green's function\n",
" for iq, q in enumerate(momentum_range):\n",
" momentum_phases = np.exp(\n",
" -1j * q * np.arange(-center, center + 2, 1)\n",
" )\n",
" result[iw, iq] = np.imag(np.dot(rgf_omega, momentum_phases))\n",
" result = -(result + result[:, ::-1]) / 2\n",
" result = np.clip(result, a_min=0, a_max=None)\n",
" return result\n",
"\n",
"\n",
"def plot_dsf(\n",
" dsf: NDArray[np.floating],\n",
" time_step: float,\n",
" n_points_momentum: int,\n",
" n_points_frequency: int,\n",
" title: str | None = None,\n",
") -> None:\n",
" \"\"\"Heat-map of the dynamical structure factor.\n",
"\n",
" Args:\n",
" dsf: DSF array of shape ``(n_points_frequency, n_points_momentum)``.\n",
" time_step: Trotter time-step size.\n",
" n_points_momentum: Number of momentum points.\n",
" n_points_frequency: Number of frequency points.\n",
" title: Optional plot title.\n",
" \"\"\"\n",
" max_frequency = np.pi / time_step\n",
" momentum_range = np.linspace(0, 2 * np.pi, n_points_momentum)\n",
" frequency_range = np.linspace(0, max_frequency, n_points_frequency)\n",
" x, y = np.meshgrid(momentum_range, frequency_range)\n",
" fig, ax = plt.subplots(figsize=(8, 5))\n",
" c = ax.pcolormesh(x, y, dsf / np.max(dsf), cmap=\"viridis\", shading=\"auto\")\n",
" fig.colorbar(c, ax=ax, label=\"Normalized intensity\")\n",
" ax.set_ylim(0, 3.6)\n",
" ax.set_xlim(0, 2 * np.pi)\n",
" ax.set_xlabel(r\"$q$\", fontsize=16)\n",
" ax.set_ylabel(r\"$\\tilde{\\omega} = \\omega / J$\", fontsize=16)\n",
" ax.set_xticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi])\n",
" ax.set_xticklabels([\"0\", r\"$\\pi/2$\", r\"$\\pi$\", r\"$3\\pi/2$\", r\"$2\\pi$\"])\n",
" if title:\n",
" ax.set_title(title, fontsize=14)\n",
" plt.tight_layout()\n",
" plt.show()\n",
"\n",
"\n",
"def plot_rgf(\n",
" n_qubits: int,\n",
" rgf_mat: NDArray[np.floating],\n",
" time_step: float,\n",
" n_steps: int,\n",
" title: str | None = None,\n",
") -> None:\n",
" \"\"\"Heat-map of the retarded Green's function in real space and time.\n",
"\n",
" Args:\n",
" n_qubits: Number of qubits (sites).\n",
" rgf_mat: RGF matrix of shape ``(n_steps, n_qubits)``.\n",
" time_step: Trotter time-step size.\n",
" n_steps: Number of time steps.\n",
" title: Optional plot title.\n",
" \"\"\"\n",
" fig, ax = plt.subplots(figsize=(8, 6))\n",
" qubit_axis = np.arange(n_qubits)\n",
" t_axis = np.arange(1, n_steps + 1) * time_step\n",
" x, y = np.meshgrid(qubit_axis, t_axis)\n",
" c = ax.pcolormesh(\n",
" x,\n",
" y,\n",
" np.real(rgf_mat),\n",
" cmap=\"RdBu\",\n",
" vmax=0.5,\n",
" vmin=-0.5,\n",
" shading=\"auto\",\n",
" )\n",
" fig.colorbar(c, ax=ax, label=r\"Re $G^R(j, j_c, t)$\")\n",
" ax.set_xlabel(\"Qubit\", fontsize=16)\n",
" ax.xaxis.set_major_locator(\n",
" plt.matplotlib.ticker.MaxNLocator(integer=True)\n",
" )\n",
" ax.set_ylabel(r\"Time\", fontsize=16)\n",
" if title:\n",
" ax.set_title(title, fontsize=14)\n",
" plt.tight_layout()\n",
" plt.show()\n",
"\n",
"\n",
"def uniform_2q_depth(circuit: QuantumCircuit) -> int:\n",
" \"\"\"Two-qubit gate depth in a standardized basis.\"\"\"\n",
" pass_manager = generate_preset_pass_manager(\n",
" optimization_level=0, basis_gates=[\"cz\", \"id\", \"rz\", \"sx\", \"x\"]\n",
" )\n",
" return pass_manager.run(circuit).depth(\n",
" lambda inst: inst.operation.num_qubits == 2\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "d28d01a6",
"metadata": {},
"source": [
"\n",
"\n",
"## Ejemplo de simulador a pequeña escala\n",
"\n",
"En primer lugar, demostramos el flujo de trabajo completo con **10 qubits**, optimizando el ansatz del estado fundamental HVA mediante simulación MPS y utilizando el simulador de vectores de estado de Qiskit para la evolución temporal. El DMRG proporciona una energía de referencia del estado fundamental y un MPS de referencia. Este ejemplo a pequeña escala nos permite validar cada paso antes de ampliarlo.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "3818ed8f",
"metadata": {},
"source": [
"\n",
"\n",
"### Paso 1: Asignar entradas clásicas a un problema cuántico\n",
"\n",
"Comenzamos definiendo el modelo físico y construyendo los circuitos cuánticos.\n",
"\n",
"**Hamiltoniano.** KCuF$_3$ se modela mediante el hamiltoniano XXZ de 1D en el punto isotrópico ( $J = 1$, $\\epsilon = 1$, tomando $J = 1$ como unidad de energía).\n",
"\n",
"**Estado fundamental.** Utilizamos un circuito de enfoque variacional `build_ground_state_ansatz` hamiltoniano (HVA) como circuito de preparación del estado fundamental. `CircuitMPS`Los parámetros del HVA se optimizan de forma clásica maximizando la fidelidad del estado $|\\langle\\psi_{\\mathrm{HVA}}(\\theta)|\\psi_{\\mathrm{DMRG}}\\rangle|^2$ con el MPS del estado fundamental del DMRG, donde el estado del HVA se evalúa simulando el circuito como un estado de producto matricial con quimb's. Utilizamos `scipy.optimize.minimize` para la optimización. También se calcula, a modo de referencia, la energía $\\langle H \\rangle$ del ansatz optimizado.\n",
"\n",
"**Puertas para trotones.** Cada término de interacción de «vecino más cercano» $e^{-i\\Delta t\\, H_{\\mathrm{pair}}}$, donde $H_{\\mathrm{pair}} = (J/4)(XX + YY + ZZ)$, se construye con `PauliEvolutionGate(H_pair, time=time_step)`. Qiskit sintetiza esto en la descomposición óptima de tres CNOT durante la transpilación.\n",
"\n",
"**Perturbación.** Una puerta « $R_z(\\pi/2)$ » aplicada al qubit central implementa un « $U_{j_c} = \\frac{1}{\\sqrt{2}}(I - i\\sigma^z_{j_c})$ », imitando el «spin-flip» local producido por un neutrón dispersado.\n",
"\n",
"**Observables.** Construimos un observable de tipo « $\\sigma^z$ » para cada sitio de qubit. Estos `SparsePauliOp` objetos se pasan a la primitiva «Estimator» en el paso 3 para extraer un valor de « $\\langle\\sigma_i^z\\rangle$ » en cada paso temporal.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "58305fbc",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ground-state energy (DMRG): -4.258035\n",
"Optimizing ground state ansatz...\n",
"Finished optimizing ground state ansatz in 8.850687621976249 seconds.\n",
"Ground state ansatz fidelity: 0.984277\n",
"Ground state ansatz energy: -4.232565\n",
"Built 10 circuits, deepest 2q depth (uniform basis) = 163\n",
"Defined 10 Z observables.\n"
]
}
],
"source": [
"# -- Physical parameters --\n",
"n_qubits = 10\n",
"interaction = 1.0 # J\n",
"anisotropy = 1.0 # ε (isotropic point)\n",
"time_step = 0.6\n",
"n_steps = 10\n",
"mps_max_bond = 32\n",
"mps_cutoff = 1e-8\n",
"center = n_qubits // 2 - 1\n",
"\n",
"# -- Hamiltonian MPO --\n",
"ham_mpo = xxz_hamiltonian_mpo(n_qubits, interaction, anisotropy)\n",
"\n",
"# -- Reference ground-state energy via DMRG --\n",
"dmrg = qtn.DMRG2(ham_mpo)\n",
"dmrg.solve(tol=1e-8)\n",
"print(f\"Ground-state energy (DMRG): {dmrg.energy:.6f}\")\n",
"\n",
"# -- Build ground state ansatz circuit --\n",
"gs_n_layers = 3\n",
"gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)\n",
"\n",
"# -- Optimize ground state ansatz parameters to maximize fidelity with the DMRG MPS --\n",
"# Initialize odd-layer angles near 0 (where the inter-pair gate is the\n",
"# identity) and even-layer angles near pi/2 (where the intra-pair gate\n",
"# is a SWAP, since 0.5*(XX + YY + ZZ) = SWAP - I/2).\n",
"rng = np.random.default_rng(12345)\n",
"x0 = np.tile([0, np.pi / 2], gs_n_layers) + rng.normal(\n",
" scale=0.1, size=2 * gs_n_layers\n",
")\n",
"print(\"Optimizing ground state ansatz...\")\n",
"t0 = timeit.default_timer()\n",
"result = optimize_ground_state_ansatz(\n",
" gs_ansatz,\n",
" x0,\n",
" dmrg.state,\n",
" max_bond=mps_max_bond,\n",
" cutoff=mps_cutoff,\n",
" options=dict(maxiter=100),\n",
")\n",
"t1 = timeit.default_timer()\n",
"print(f\"Finished optimizing ground state ansatz in {t1 - t0} seconds.\")\n",
"print(f\"Ground state ansatz fidelity: {1 - result.fun:.6f}\")\n",
"\n",
"gs_circuit = gs_ansatz.assign_parameters(result.x)\n",
"gs_circuit_mps = quimb_circuit(\n",
" gs_circuit.decompose([\"PauliEvolution\"]),\n",
" quimb_circuit_class=qtn.CircuitMPS,\n",
" max_bond=mps_max_bond,\n",
" cutoff=mps_cutoff,\n",
")\n",
"gs_ansatz_energy = qtn.expec_TN_1D(\n",
" gs_circuit_mps.psi.H, ham_mpo, gs_circuit_mps.psi\n",
")\n",
"print(f\"Ground state ansatz energy: {gs_ansatz_energy:.6f}\")\n",
"\n",
"# -- Build circuits for each time step --\n",
"perturbed = gs_circuit.copy()\n",
"perturbed.rz(np.pi / 2, center)\n",
"\n",
"circuits = []\n",
"for t in range(1, n_steps + 1):\n",
" circuit = perturbed.copy()\n",
" for instr in trotter_evolution(\n",
" circuit.qubits, interaction, anisotropy, time_step, t\n",
" ):\n",
" circuit.append(instr)\n",
" circuits.append(circuit)\n",
"print(\n",
" f\"Built {len(circuits)} circuits, deepest 2q depth (uniform basis) = \"\n",
" f\"{uniform_2q_depth(circuits[-1])}\"\n",
")\n",
"\n",
"# -- Observables: Z on each qubit site --\n",
"observables = [\n",
" SparsePauliOp.from_sparse_list([(\"Z\", [i], 1)], num_qubits=n_qubits)\n",
" for i in range(n_qubits)\n",
"]\n",
"print(f\"Defined {len(observables)} Z observables.\")"
]
},
{
"cell_type": "markdown",
"id": "39191552",
"metadata": {},
"source": [
"\n",
"\n",
"### Paso 2: Optimizar el problema para su ejecución en hardware cuántico\n",
"\n",
"En el caso de un hardware real, los circuitos de Trotter mencionados anteriormente serían demasiado complejos. **La compilación cuántica aproximada (AQC)** resuelve este problema sustituyendo las primeras capas de Trotter de « $k$ » (incluido el circuito del estado fundamental) por un ansatz parametrizado más corto, optimizado para maximizar la fidelidad a nivel de MPS con respecto al circuito profundo original. Los pasos restantes de Trotter se añaden tal cual, lo que da lugar a un circuito «AQC + Trotter» menos profundo.\n",
"\n",
"Para equilibrar la expresividad y la profundidad del circuito, se utilizan dos enfoques: un **enfoque de una sola capa** (generado a partir de un único paso de Trotter) comprime los primeros pasos temporales, y un **enfoque** más profundo de dos capas (generado a partir de dos pasos de Trotter) comprime los siguientes pasos, en los que se requiere una mayor fidelidad.\n",
"\n",
"El proceso de AQC consta de cuatro pasos secundarios:\n",
"\n",
"1. **Construir circuitos objetivo** : los primeros circuitos « $k_1 + k_2$ » del paso 1 sirven directamente como objetivos de AQC.\n",
"2. **Calcular el MPS de destino** : simular cada circuito de destino como un estado de producto matricial utilizando `quimb.tensor.CircuitMPS`.\n",
"3. **Generar y optimizar ansätze** : `generate_ansatz_from_circuit` crea un ansatz parametrizado de una capa y otro de dos capas; los parámetros se optimizan mediante L-BFGS-B con gradientes acelerados por JAX para minimizar $1 - |\\langle\\psi_{\\mathrm{ansatz}}|\\psi_{\\mathrm{target}}\\rangle|^2$. Los primeros pasos $k_1$ utilizan el ansatz de una capa y los siguientes pasos $k_2$ utilizan el ansatz de dos capas; dentro de cada etapa, cada paso parte de los parámetros optimizados del paso anterior, y los parámetros se restablecen a los valores predeterminados de la etapa al final de esta.\n",
"4. **Montar circuitos mixtos** : para los pasos temporales posteriores al punto de control del AQC, añadir capas exactas de Trotter al circuito AQC optimizado de dos capas utilizando `trotter_evolution`.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "e78f8d9a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Step 2b — target MPS:\n",
" k=1: max bond = 22\n",
" k=2: max bond = 22\n",
" k=3: max bond = 26\n",
" k=4: max bond = 27\n",
" k=5: max bond = 30\n",
"\n",
"Step 2c — one-layer ansatz: 389 parameters, 2q depth (uniform basis) = 27\n",
"Step 2c — two-layer ansatz: 470 parameters, 2q depth (uniform basis) = 33\n",
" k=1: fidelity = 1.0000, 2q depth (uniform basis) = 27, 11.8s\n",
" k=2: fidelity = 0.9980, 2q depth (uniform basis) = 27, 13.2s\n",
" k=3: fidelity = 0.9890, 2q depth (uniform basis) = 27, 13.7s\n",
" k=4: fidelity = 0.9983, 2q depth (uniform basis) = 33, 16.4s\n",
" k=5: fidelity = 0.9957, 2q depth (uniform basis) = 33, 16.4s\n",
"\n",
"Step 2d — assembled 10 circuits\n",
" At step 10 (uniform basis): Trotter 2q depth = 163, AQC+Trotter 2q depth = 99; 2q gates: Trotter = 371, AQC+Trotter = 200\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Number of time steps to compress into an AQC ansatz with one layer\n",
"aqc_n_steps_1 = 3\n",
"# Number of time steps to compress into an AQC ansatz with two layers\n",
"aqc_n_steps_2 = 2\n",
"aqc_n_steps_total = aqc_n_steps_1 + aqc_n_steps_2\n",
"\n",
"# ── Step 2a: Target circuits (first aqc_n_steps_total circuits from Step 1) ──\n",
"target_circuits = {\n",
" k: circuits[k - 1] for k in range(1, aqc_n_steps_total + 1)\n",
"}\n",
"# ── Step 2b: Compute target MPS ──\n",
"# Decompose PauliEvolutionGate to RXX/RYY/RZZ before passing to the AQC MPS\n",
"# backend, which does not understand PauliEvolutionGate natively.\n",
"aqc_sim = QuimbSimulator(\n",
" quimb_circuit_factory=partial(\n",
" qtn.CircuitMPS,\n",
" gate_opts=dict(max_bond=mps_max_bond, cutoff=mps_cutoff),\n",
" ),\n",
" autodiff_backend=\"jax\",\n",
")\n",
"print(\"\\nStep 2b — target MPS:\")\n",
"target_mps = {}\n",
"for k in range(1, aqc_n_steps_total + 1):\n",
" target_mps[k] = tensornetwork_from_circuit(\n",
" target_circuits[k].decompose([\"PauliEvolution\"]), aqc_sim\n",
" )\n",
" print(f\" k={k}: max bond = {target_mps[k].psi.max_bond()}\")\n",
"\n",
"# ── Step 2c: Generate ansätze (one-layer and two-layer) and optimize parameters ──\n",
"ansatz_1, initial_params_1 = generate_ansatz_from_circuit(\n",
" target_circuits[1].decompose([\"PauliEvolution\"]),\n",
" qubits_initially_zero=True,\n",
")\n",
"initial_params_1 = np.array(initial_params_1)\n",
"ansatz_2, initial_params_2 = generate_ansatz_from_circuit(\n",
" target_circuits[2].decompose([\"PauliEvolution\"]),\n",
" qubits_initially_zero=True,\n",
")\n",
"initial_params_2 = np.array(initial_params_2)\n",
"print(\n",
" f\"\\nStep 2c — one-layer ansatz: {ansatz_1.num_parameters} parameters, \"\n",
" f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_1)}\"\n",
")\n",
"print(\n",
" f\"Step 2c — two-layer ansatz: {ansatz_2.num_parameters} parameters, \"\n",
" f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_2)}\"\n",
")\n",
"\n",
"aqc_circuits = {}\n",
"aqc_params = {}\n",
"for k in range(1, aqc_n_steps_total + 1):\n",
" if k <= aqc_n_steps_1:\n",
" ansatz, base_params = ansatz_1, initial_params_1\n",
" else:\n",
" ansatz, base_params = ansatz_2, initial_params_2\n",
" # Warm-start from the previous step only within the same stage\n",
" same_stage = (k - 1 >= 1) and (\n",
" (k - 1 <= aqc_n_steps_1) == (k <= aqc_n_steps_1)\n",
" )\n",
" x0 = aqc_params[k - 1] if same_stage else base_params\n",
" obj = MaximizeStateFidelity(target_mps[k], ansatz, aqc_sim)\n",
" t0 = timeit.default_timer()\n",
" result = scipy.optimize.minimize(\n",
" obj.loss_function,\n",
" x0,\n",
" method=\"L-BFGS-B\",\n",
" jac=True,\n",
" options=dict(maxiter=100),\n",
" )\n",
" elapsed = timeit.default_timer() - t0\n",
" aqc_params[k] = result.x\n",
" aqc_circuits[k] = ansatz.assign_parameters(result.x)\n",
" print(\n",
" f\" k={k}: fidelity = {1 - result.fun:.4f}, \"\n",
" f\"2q depth (uniform basis) = \"\n",
" f\"{uniform_2q_depth(aqc_circuits[k])}, \"\n",
" f\"{elapsed:.1f}s\"\n",
" )\n",
"\n",
"# ── Step 2d: Assemble full circuit set (AQC + Trotter) ──\n",
"all_circuits = []\n",
"for k in range(1, aqc_n_steps_total + 1):\n",
" all_circuits.append(aqc_circuits[k])\n",
"base = aqc_circuits[aqc_n_steps_total]\n",
"for k in range(1, n_steps - aqc_n_steps_total + 1):\n",
" circuit = base.copy()\n",
" for instr in trotter_evolution(\n",
" circuit.qubits, interaction, anisotropy, time_step, k\n",
" ):\n",
" circuit.append(instr)\n",
" all_circuits.append(circuit)\n",
"\n",
"full_depths = [\n",
" uniform_2q_depth(circuits[k - 1]) for k in range(1, n_steps + 1)\n",
"]\n",
"aqc_depths = [uniform_2q_depth(circuit) for circuit in all_circuits]\n",
"full_2q = [\n",
" circuits[k - 1].decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n",
" for k in range(1, n_steps + 1)\n",
"]\n",
"aqc_2q = [\n",
" circuit.decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n",
" for circuit in all_circuits\n",
"]\n",
"print(f\"\\nStep 2d — assembled {len(all_circuits)} circuits\")\n",
"print(\n",
" f\" At step {n_steps} (uniform basis): \"\n",
" f\"Trotter 2q depth = {full_depths[-1]}, AQC+Trotter 2q depth = {aqc_depths[-1]}; \"\n",
" f\"2q gates: Trotter = {full_2q[-1]}, AQC+Trotter = {aqc_2q[-1]}\"\n",
")\n",
"\n",
"steps_axis = np.arange(1, n_steps + 1)\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.plot(steps_axis, full_depths, \"-o\", color=\"black\", label=\"Full Trotter\")\n",
"ax.plot(\n",
" steps_axis, aqc_depths, \"-o\", color=\"cadetblue\", label=\"AQC + Trotter\"\n",
")\n",
"ax.set_xlabel(\"Trotter steps\", fontsize=13)\n",
"ax.xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=True))\n",
"ax.set_ylabel(\"2q gate depth (uniform basis)\", fontsize=13)\n",
"ax.set_title(f\"AQC circuit-depth reduction ({n_qubits} qubits)\", fontsize=13)\n",
"ax.legend(fontsize=11)\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "43ad6048",
"metadata": {},
"source": [
"\n",
"\n",
"### Paso 3: Ejecutar con el comando « Qiskit primitives »\n",
"\n",
"Simulamos cada circuito compilado con AQC utilizando la `StatevectorEstimator` primitiva.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "b8f9e41c",
"metadata": {},
"outputs": [],
"source": [
"estimator = StatevectorEstimator()\n",
"pubs = [(circuit, observables) for circuit in all_circuits]\n",
"job = estimator.run(pubs)\n",
"result = job.result()"
]
},
{
"cell_type": "markdown",
"id": "cb12662c",
"metadata": {},
"source": [
"\n",
"\n",
"### Paso 4: Realizar el posprocesamiento y obtener el resultado en el formato clásico deseado\n",
"\n",
"A continuación, extraemos el valor esperado $\\langle\\sigma_i^z\\rangle$ para cada qubit $i$ en cada paso temporal $t$. Estos valores forman la matriz de la función de Green retardada $G^R(j, j_c, t)$. A continuación, aplicamos la transformada de Fourier a la función de Green retardada (RGF) para obtener el factor de estructura dinámico $S(q, \\omega)$ y representamos gráficamente tanto la RGF como el DSF. `get_dsf` aplica la simetría especular y recorta los valores negativos internamente.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "4585d139",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rgf_mat = np.stack([pub_result.data.evs for pub_result in result])\n",
"\n",
"# -- Compute DSF --\n",
"n_points_momentum, n_points_frequency = 100, 100\n",
"spectrum = get_dsf(\n",
" n_qubits,\n",
" rgf_mat,\n",
" time_step,\n",
" n_steps,\n",
" n_points_momentum,\n",
" n_points_frequency,\n",
")\n",
"\n",
"# -- Plot retarded Green's function --\n",
"plot_rgf(\n",
" n_qubits,\n",
" rgf_mat,\n",
" time_step,\n",
" n_steps,\n",
" title=f\"Retarded Green's function — {n_qubits} qubits (simulation)\",\n",
")\n",
"\n",
"# -- Plot DSF --\n",
"plot_dsf(\n",
" spectrum,\n",
" time_step,\n",
" n_points_momentum,\n",
" n_points_frequency,\n",
" title=f\"Dynamical structure factor — {n_qubits} qubits (simulation)\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "ls-intro",
"metadata": {},
"source": [
"\n",
"\n",
"## Ejecución de hardware a gran escala\n",
"\n",
"Ahora ampliamos la escala hasta **los 50 qubits**. A esta escala, las optimizaciones alcanzan fidelidades inferiores a las del ejemplo a pequeña escala: la fidelidad del ansatz del estado fundamental desciende hasta aproximadamente 0.65, y las fidelidades de AQC disminuyen hasta aproximadamente 0.7 en los últimos puntos de control. Esto es de esperar, y se puede mejorar la precisión aumentando el número de capas del ansatz del estado fundamental (`gs_n_layers`) o el número de iteraciones de optimización (`maxiter`) con un coste clásico adicional. Cabe señalar también que las fidelidades del AQC de dos capas resultan inferiores a las del de una sola capa. No se trata de una regresión: los pasos temporales posteriores generan más entrelazamiento y son simplemente más difíciles de comprimir, por lo que para ellos se utiliza el ansatz de dos capas, que es más expresivo.\n",
"\n",
"La optimización de AQC también puede requerir varias horas de tiempo de cálculo clásico (aproximadamente seis horas en la ejecución que se muestra aquí, la mayor parte de las cuales se dedica a los puntos de control de dos capas del paso 2c ). Para reducir el tiempo real, plantéate ejecutar este cuaderno en un hardware clásico más potente, como un sistema de computación de alto rendimiento (HPC). Como alternativa, puedes reducir la escala a una instancia del problema más pequeña (por ejemplo, con menos qubits o menos pasos temporales); en ese caso, tus resultados diferirán de los que se muestran aquí.\n",
"\n",
"El código que aparece a continuación sigue la misma estructura de cuatro pasos que el ejemplo a pequeña escala. Los parámetros del estado fundamental de la HVA se optimizan de nuevo mediante una simulación MPS. En la QPU, activamos el desacoplamiento dinámico (DD), el «Pauli twirling» y la extinción de errores de lectura con «twirling» (TREX) para suprimir y mitigar los errores. La siguiente tabla resume las diferencias entre el experimento a gran escala y el de pequeña escala:\n",
"\n",
"| | PEQUEÑA ESCALA | GRAN ESCALA |\n",
"| ------------------------------------------- | ---------------------- | --------------------------------- |\n",
"| Qubits | 10 | 50 |\n",
"| Intervalos de tiempo | 10 | 20 |\n",
"| Puntos de control de AQC (1 capa + 2 capas) | 3 + 2 = 5 | 6 + 4 = 10 |\n",
"| Capas del ansatz del estado fundamental | 3 | 5 |\n",
"| Dimensión máxima de la unión MPS | 32 | 128 |\n",
"| Estimador | `StatevectorEstimator` | QPU con DD, giros de Pauli y TREX |\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9e85fbd1",
"metadata": {},
"source": [
"\n",
" Durante la optimización de AQC (paso « 2c », que se describe a continuación), es posible que aparezcan `stderr` mensajes como los siguientes:\n",
"\n",
" ```\n",
" [Compiling module jit_MakeArrayFn for CPU] Very slow compile? ...\n",
" The operation took 2m14s\n",
" ```\n",
"\n",
" Se trata de mensajes de diagnóstico benignos de XLA, el compilador que sustenta `qiskit-addon-aqc-tensor`la función «autodiff» de JAX. Con 50 qubits y una dimensión de enlace MPS de 128, XLA tarda un par de minutos en compilar la función de gradiente la primera vez que se traza. La compilación se realiza correctamente y los resultados de la optimización no se ven afectados.\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "fbe57e1a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ground-state energy (DMRG): -21.972109\n",
"Optimizing ground state ansatz...\n",
"Finished optimizing ground state ansatz in 132.2076231740648 seconds.\n",
"Ground state ansatz fidelity: 0.645956\n",
"Ground state ansatz energy: -21.616744\n",
"Built 20 circuits, deepest 2q depth (uniform basis) = 307\n",
"Defined 50 Z observables.\n",
"\n",
"Step 2b — target MPS:\n",
" k=1: max bond = 44\n",
" k=2: max bond = 46\n",
" k=3: max bond = 53\n",
" k=4: max bond = 62\n",
" k=5: max bond = 75\n",
" k=6: max bond = 96\n",
" k=7: max bond = 118\n",
" k=8: max bond = 128\n",
" k=9: max bond = 128\n",
" k=10: max bond = 128\n",
"\n",
"Step 2c — one-layer ansatz: 2971 parameters, 2q depth (uniform basis) = 39\n",
"Step 2c — two-layer ansatz: 3412 parameters, 2q depth (uniform basis) = 45\n",
" k=1: fidelity = 1.0000, 2q depth (uniform basis) = 39, 215.6s\n",
" k=2: fidelity = 0.9676, 2q depth (uniform basis) = 39, 269.8s\n",
" k=3: fidelity = 0.9014, 2q depth (uniform basis) = 39, 293.5s\n",
" k=4: fidelity = 0.8755, 2q depth (uniform basis) = 39, 279.6s\n",
" k=5: fidelity = 0.8450, 2q depth (uniform basis) = 39, 266.4s\n",
" k=6: fidelity = 0.8146, 2q depth (uniform basis) = 39, 371.0s\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"E0626 02:59:07.229070 1508496 slow_operation_alarm.cc:73] \n",
"********************************\n",
"[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.\n",
"********************************\n",
"E0626 02:59:21.713796 1508477 slow_operation_alarm.cc:140] The operation took 2m14.484932074s\n",
"\n",
"********************************\n",
"[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.\n",
"********************************\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" k=7: fidelity = 0.7312, 2q depth (uniform basis) = 45, 10865.3s\n",
" k=8: fidelity = 0.7720, 2q depth (uniform basis) = 45, 1450.3s\n",
" k=9: fidelity = 0.7316, 2q depth (uniform basis) = 45, 3388.6s\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"E0626 07:21:05.148724 1508496 slow_operation_alarm.cc:73] \n",
"********************************\n",
"[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.\n",
"********************************\n",
"E0626 07:21:24.286095 1508477 slow_operation_alarm.cc:140] The operation took 2m19.137566498s\n",
"\n",
"********************************\n",
"[Compiling module jit_MakeArrayFn for CPU] Very slow compile? If you want to file a bug, run with envvar XLA_FLAGS=--xla_dump_to=/tmp/foo and attach the results.\n",
"********************************\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" k=10: fidelity = 0.6726, 2q depth (uniform basis) = 45, 3396.1s\n",
"\n",
"Step 2d — assembled 20 circuits\n",
" At step 20 (uniform basis): Trotter 2q depth = 307, AQC+Trotter 2q depth = 171; 2q gates: Trotter = 3775, AQC+Trotter = 1913\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Backend: ibm_fez\n",
"Transpiled 2q depth (deepest, ISA on ibm_fez): 105 (2q gates: 2574)\n",
"Job ID: d8v39vhropqc738biotg\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# ── Parameters ──────────────────────────────────────────────────────────────\n",
"n_qubits = 50 # 10 → 50\n",
"interaction = 1.0 # J\n",
"anisotropy = 1.0 # ε (isotropic point)\n",
"time_step = 0.6\n",
"n_steps = 20 # 10 → 20\n",
"aqc_n_steps_1 = 6 # 3 → 6\n",
"aqc_n_steps_2 = 4 # 2 → 4\n",
"aqc_n_steps_total = aqc_n_steps_1 + aqc_n_steps_2\n",
"gs_n_layers = 5 # 3 → 5\n",
"mps_max_bond = 128 # 32 -> 128\n",
"mps_cutoff = 1e-8\n",
"center = n_qubits // 2 - 1\n",
"\n",
"# ── Step 1: Map ──────────────────────────────────────────────────────────────\n",
"ham_mpo = xxz_hamiltonian_mpo(n_qubits, interaction, anisotropy)\n",
"\n",
"dmrg = qtn.DMRG2(ham_mpo)\n",
"dmrg.solve(tol=1e-8)\n",
"print(f\"Ground-state energy (DMRG): {dmrg.energy:.6f}\")\n",
"\n",
"gs_ansatz = build_ground_state_ansatz(n_qubits, gs_n_layers)\n",
"\n",
"rng = np.random.default_rng(12345)\n",
"x0 = np.tile([0.0, np.pi / 2], gs_n_layers) + rng.normal(\n",
" scale=0.1, size=2 * gs_n_layers\n",
")\n",
"print(\"Optimizing ground state ansatz...\")\n",
"t0 = timeit.default_timer()\n",
"result = optimize_ground_state_ansatz(\n",
" gs_ansatz,\n",
" x0,\n",
" dmrg.state,\n",
" max_bond=mps_max_bond,\n",
" cutoff=mps_cutoff,\n",
" options=dict(maxiter=100),\n",
")\n",
"t1 = timeit.default_timer()\n",
"print(f\"Finished optimizing ground state ansatz in {t1 - t0} seconds.\")\n",
"print(f\"Ground state ansatz fidelity: {1 - result.fun:.6f}\")\n",
"\n",
"gs_circuit = gs_ansatz.assign_parameters(result.x)\n",
"gs_circuit_mps = quimb_circuit(\n",
" gs_circuit.decompose([\"PauliEvolution\"]),\n",
" quimb_circuit_class=qtn.CircuitMPS,\n",
" max_bond=mps_max_bond,\n",
" cutoff=mps_cutoff,\n",
")\n",
"gs_ansatz_energy = qtn.expec_TN_1D(\n",
" gs_circuit_mps.psi.H, ham_mpo, gs_circuit_mps.psi\n",
")\n",
"print(f\"Ground state ansatz energy: {gs_ansatz_energy:.6f}\")\n",
"\n",
"perturbed = gs_circuit.copy()\n",
"perturbed.rz(np.pi / 2, center)\n",
"\n",
"circuits = []\n",
"for t in range(1, n_steps + 1):\n",
" circuit = perturbed.copy()\n",
" for instr in trotter_evolution(\n",
" circuit.qubits, interaction, anisotropy, time_step, t\n",
" ):\n",
" circuit.append(instr)\n",
" circuits.append(circuit)\n",
"print(\n",
" f\"Built {len(circuits)} circuits, deepest 2q depth (uniform basis) = \"\n",
" f\"{uniform_2q_depth(circuits[-1])}\"\n",
")\n",
"\n",
"observables = [\n",
" SparsePauliOp.from_sparse_list([(\"Z\", [i], 1)], num_qubits=n_qubits)\n",
" for i in range(n_qubits)\n",
"]\n",
"print(f\"Defined {len(observables)} Z observables.\")\n",
"\n",
"# ── Step 2: AQC ──────────────────────────────────────────────────────────────\n",
"target_circuits = {\n",
" k: circuits[k - 1] for k in range(1, aqc_n_steps_total + 1)\n",
"}\n",
"\n",
"aqc_sim = QuimbSimulator(\n",
" quimb_circuit_factory=partial(\n",
" qtn.CircuitMPS,\n",
" gate_opts=dict(max_bond=mps_max_bond, cutoff=mps_cutoff),\n",
" ),\n",
" autodiff_backend=\"jax\",\n",
")\n",
"print(\"\\nStep 2b — target MPS:\")\n",
"target_mps = {}\n",
"for k in range(1, aqc_n_steps_total + 1):\n",
" target_mps[k] = tensornetwork_from_circuit(\n",
" target_circuits[k].decompose([\"PauliEvolution\"]), aqc_sim\n",
" )\n",
" print(f\" k={k}: max bond = {target_mps[k].psi.max_bond()}\")\n",
"\n",
"ansatz_1, initial_params_1 = generate_ansatz_from_circuit(\n",
" target_circuits[1].decompose([\"PauliEvolution\"]),\n",
" qubits_initially_zero=True,\n",
")\n",
"initial_params_1 = np.array(initial_params_1)\n",
"ansatz_2, initial_params_2 = generate_ansatz_from_circuit(\n",
" target_circuits[2].decompose([\"PauliEvolution\"]),\n",
" qubits_initially_zero=True,\n",
")\n",
"initial_params_2 = np.array(initial_params_2)\n",
"print(\n",
" f\"\\nStep 2c — one-layer ansatz: {ansatz_1.num_parameters} parameters, \"\n",
" f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_1)}\"\n",
")\n",
"print(\n",
" f\"Step 2c — two-layer ansatz: {ansatz_2.num_parameters} parameters, \"\n",
" f\"2q depth (uniform basis) = {uniform_2q_depth(ansatz_2)}\"\n",
")\n",
"\n",
"aqc_circuits = {}\n",
"aqc_params = {}\n",
"for k in range(1, aqc_n_steps_total + 1):\n",
" if k <= aqc_n_steps_1:\n",
" ansatz, base_params = ansatz_1, initial_params_1\n",
" else:\n",
" ansatz, base_params = ansatz_2, initial_params_2\n",
" # Warm-start from the previous step only within the same stage\n",
" same_stage = (k - 1 >= 1) and (\n",
" (k - 1 <= aqc_n_steps_1) == (k <= aqc_n_steps_1)\n",
" )\n",
" x0 = aqc_params[k - 1] if same_stage else base_params\n",
" obj = MaximizeStateFidelity(target_mps[k], ansatz, aqc_sim)\n",
" t0 = timeit.default_timer()\n",
" result = scipy.optimize.minimize(\n",
" obj.loss_function,\n",
" x0,\n",
" method=\"L-BFGS-B\",\n",
" jac=True,\n",
" options=dict(maxiter=100),\n",
" )\n",
" elapsed = timeit.default_timer() - t0\n",
" aqc_params[k] = result.x\n",
" aqc_circuits[k] = ansatz.assign_parameters(result.x)\n",
" print(\n",
" f\" k={k}: fidelity = {1 - result.fun:.4f}, \"\n",
" f\"2q depth (uniform basis) = \"\n",
" f\"{uniform_2q_depth(aqc_circuits[k])}, \"\n",
" f\"{elapsed:.1f}s\"\n",
" )\n",
"\n",
"all_circuits = []\n",
"for k in range(1, aqc_n_steps_total + 1):\n",
" all_circuits.append(aqc_circuits[k])\n",
"base = aqc_circuits[aqc_n_steps_total]\n",
"for k in range(1, n_steps - aqc_n_steps_total + 1):\n",
" circuit = base.copy()\n",
" for instr in trotter_evolution(\n",
" circuit.qubits, interaction, anisotropy, time_step, k\n",
" ):\n",
" circuit.append(instr)\n",
" all_circuits.append(circuit)\n",
"\n",
"full_depths = [\n",
" uniform_2q_depth(circuits[k - 1]) for k in range(1, n_steps + 1)\n",
"]\n",
"aqc_depths = [uniform_2q_depth(circuit) for circuit in all_circuits]\n",
"full_2q = [\n",
" circuits[k - 1].decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n",
" for k in range(1, n_steps + 1)\n",
"]\n",
"aqc_2q = [\n",
" circuit.decompose([\"PauliEvolution\"]).num_nonlocal_gates()\n",
" for circuit in all_circuits\n",
"]\n",
"print(f\"\\nStep 2d — assembled {len(all_circuits)} circuits\")\n",
"print(\n",
" f\" At step {n_steps} (uniform basis): \"\n",
" f\"Trotter 2q depth = {full_depths[-1]}, AQC+Trotter 2q depth = {aqc_depths[-1]}; \"\n",
" f\"2q gates: Trotter = {full_2q[-1]}, AQC+Trotter = {aqc_2q[-1]}\"\n",
")\n",
"\n",
"steps_axis = np.arange(1, n_steps + 1)\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.plot(steps_axis, full_depths, \"-o\", color=\"black\", label=\"Full Trotter\")\n",
"ax.plot(\n",
" steps_axis, aqc_depths, \"-o\", color=\"cadetblue\", label=\"AQC + Trotter\"\n",
")\n",
"ax.set_xlabel(\"Trotter steps\", fontsize=13)\n",
"ax.xaxis.set_major_locator(plt.matplotlib.ticker.MaxNLocator(integer=True))\n",
"ax.set_ylabel(\"2q gate depth (uniform basis)\", fontsize=13)\n",
"ax.set_title(f\"AQC circuit-depth reduction ({n_qubits} qubits)\", fontsize=13)\n",
"ax.legend(fontsize=11)\n",
"plt.tight_layout()\n",
"plt.show()\n",
"\n",
"# ── Step 3: Execute on IBM Quantum hardware ───────────────────────────────────\n",
"# (replaces StatevectorEstimator)\n",
"service = QiskitRuntimeService()\n",
"backend = service.least_busy(\n",
" min_num_qubits=n_qubits,\n",
" operational=True,\n",
" simulator=False,\n",
" filters=lambda x: x.configuration().processor_type[\"family\"] == \"Heron\",\n",
")\n",
"print(f\"Backend: {backend.name}\")\n",
"\n",
"pm = generate_preset_pass_manager(optimization_level=3, backend=backend)\n",
"isa_circuits = pm.run(all_circuits, num_processes=1)\n",
"isa_2q_depths = [\n",
" isa_circuit.depth(lambda inst: inst.operation.num_qubits == 2)\n",
" for isa_circuit in isa_circuits\n",
"]\n",
"print(\n",
" f\"Transpiled 2q depth (deepest, ISA on {backend.name}): \"\n",
" f\"{max(isa_2q_depths)} \"\n",
" f\"(2q gates: {max(isa_circuit.num_nonlocal_gates() for isa_circuit in isa_circuits)})\"\n",
")\n",
"\n",
"estimator = Estimator(backend)\n",
"estimator.options.environment.job_tags = [\"TUT_SNS\"]\n",
"estimator.options.dynamical_decoupling.enable = True\n",
"estimator.options.dynamical_decoupling.sequence_type = \"XY4\"\n",
"estimator.options.twirling.enable_gates = True\n",
"estimator.options.twirling.num_randomizations = 1000\n",
"estimator.options.twirling.shots_per_randomization = 128\n",
"estimator.options.resilience.measure_mitigation = True\n",
"estimator.options.resilience.measure_noise_learning.num_randomizations = 32\n",
"estimator.options.resilience.measure_noise_learning.shots_per_randomization = 100\n",
"\n",
"pubs = [\n",
" (\n",
" isa_circuit,\n",
" [obs.apply_layout(isa_circuit.layout) for obs in observables],\n",
" )\n",
" for isa_circuit in isa_circuits\n",
"]\n",
"job = estimator.run(pubs)\n",
"print(f\"Job ID: {job.job_id()}\")\n",
"\n",
"result = job.result()\n",
"\n",
"# ── Step 4: Post-process ──────────────────────────────────────────────────────\n",
"rgf_mat = np.stack([pub_result.data.evs for pub_result in result])\n",
"\n",
"n_points_momentum, n_points_frequency = 100, 100\n",
"spectrum = get_dsf(\n",
" n_qubits,\n",
" rgf_mat,\n",
" time_step,\n",
" n_steps,\n",
" n_points_momentum,\n",
" n_points_frequency,\n",
")\n",
"\n",
"plot_rgf(\n",
" n_qubits,\n",
" rgf_mat,\n",
" time_step,\n",
" n_steps,\n",
" title=f\"Retarded Green's function — {n_qubits} qubits (QPU)\",\n",
")\n",
"plot_dsf(\n",
" spectrum,\n",
" time_step,\n",
" n_points_momentum,\n",
" n_points_frequency,\n",
" title=rf\"KCuF$_3$ DSF — {n_qubits} qubits (QPU)\"\n",
" \"\\n(AQC + DD + Pauli twirling + TREX)\",\n",
")"
]
},
{
"cell_type": "markdown",
"id": "ls-discussion",
"metadata": {},
"source": [
"Los resultados del hardware reproducen las características clave del continuo de dos espinones: la intensidad de dispersión se concentra cerca del vector de onda antiferromagnético $q = \\pi$ a baja energía y está limitada por debajo por la dispersión sinusoidal de los espinones, con un amplio continuo de peso espectral por encima de ella, en lugar de un único modo bien definido. Se trata de la misma estructura medida mediante dispersión inelástica de neutrones en KCuF$_3$, lo que valida el flujo de trabajo completo —que incluye la preparación del estado fundamental, la perturbación, la evolución de Trotter comprimida mediante AQC y la medición con mitigación de errores— a escala de 50 qubits.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ls-nextsteps",
"metadata": {},
"source": [
"\n",
"\n",
"## Próximos pasos\n",
"\n",
"Si este trabajo te ha parecido interesante, quizá te interese el siguiente material:\n",
"\n",
"\n",
" * [Simulación de la dispersión de neutrones con un flujo de trabajo sin servidor basado en AQC y dinámica de Trotter](/docs/tutorials/simulate-neutron-scattering-with-a-serverless-workflow) : un tutorial complementario que ejecuta este mismo experimento mediante una plantilla de función de Qiskit Serverless ya implementada\n",
" * [Lee et al., «Evaluación comparativa de la simulación cuántica mediante experimentos de dispersión de neutrones» ( arXiv:2603.15608 )](https://arxiv.org/abs/2603.15608) : el artículo de referencia en el que se basa este tutorial\n",
" * [Técnicas de mitigación y supresión de errores](/docs/guides/error-mitigation-and-suppression-techniques) — DD, «Pauli twirling» y TREX utilizadas en los experimentos de hardware\n",
" * [Compilación cuántica aproximada para circuitos de evolución temporal](/docs/tutorials/approximate-quantum-compilation-for-time-evolution) : tutorial sobre AQC-Tensor\n",
" * [Documentación de AQC-Tensor](https://qiskit.github.io/qiskit-addon-aqc-tensor/)\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": 780
},
"nbformat": 4,
"nbformat_minor": 5
}