{
"cells": [
{
"cell_type": "markdown",
"id": "090b6884",
"metadata": {},
"source": [
"---\n",
"title: \"Implementación de SQD\"\n",
"description: \"La diagonalización cuántica basada en muestras (SQD) se implementa en el contexto de la resolución del estado fundamental de una molécula de nitrógeno.\"\n",
"---\n",
"\n",
"{/* cspell:ignore fontdict fontsize milli hcore ccsd Motta */}\n",
"\n",
"\n",
"\n",
"# SQD para la estimación energética de un hamiltoniano químico\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ee7de1bc",
"metadata": {},
"source": [
"En esta lección, aplicaremos SQD para estimar la energía de estado fundamental de una molécula.\n",
"\n",
"En concreto, trataremos los siguientes temas utilizando el enfoque del patrón Qiskit $4$ -step:\n",
"\n",
"1. Paso 1: Asignar el problema a circuitos y operadores cuánticos\n",
" * Configurar el Hamiltoniano molecular para $N_2$.\n",
" * Explicar el clúster unitario local Jastrow (LUCJ), inspirado en la química y fácil de usar en hardware [\\[1\\]](#references)\n",
"2. Paso 2: Optimización para el hardware de destino\n",
" * Optimizar el número de puertas y el diseño del ansatz para su ejecución en hardware\n",
"3. Paso 3: Ejecutar en el hardware de destino\n",
" * Ejecuta el circuito optimizado en una QPU real para generar muestras del subespacio.\n",
"4. Paso 4: Tratamiento posterior de los resultados\n",
" * Introducir el bucle de recuperación de configuración autoconsistente [\\[2\\]](#references)\n",
" * Post-procesar el conjunto completo de muestras de cadenas de bits, utilizando el conocimiento previo del número de partículas y la ocupación orbital media calculada en la iteración más reciente.\n",
" * Crear probabilísticamente lotes de submuestras a partir de las cadenas de bits recuperadas.\n",
" * Proyectar y diagonalizar el Hamiltoniano molecular sobre cada subespacio muestreado.\n",
" * Guarda la energía mínima del estado base encontrada en todos los lotes y actualiza la ocupación orbital media.\n",
"\n",
"Utilizaremos varios paquetes de software a lo largo de la lección.\n",
"\n",
"* `PySCF` para definir la molécula y configurar el Hamiltoniano.\n",
"* `ffsim` para construir el ansatz LUCJ.\n",
"* `Qiskit` para transpilar el ansatz para su ejecución en hardware.\n",
"* `Qiskit IBM Runtime` para ejecutar el circuito en una QPU y recoger muestras.\n",
"* `Qiskit addon SQD` recuperación de la configuración y estimación de la energía del estado base mediante proyección de subespacios y diagonalización de matrices.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "719a9c0e-8c00-4fab-ba23-f2c8a7ebb573",
"metadata": {},
"source": [
"\n",
"\n",
"## 1. Asignar el problema a circuitos y operadores cuánticos\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d02f97af",
"metadata": {},
"source": [
"\n",
"\n",
"### Hamiltoniano molecular\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "5dac27ce-56b1-4f7f-ab83-ac153c004e82",
"metadata": {},
"source": [
"Un Hamiltoniano molecular toma la forma genérica:\n",
"\n",
"$$\n",
"\\hat{H} = \\sum_{ \\substack{pr\\\\\\sigma} } h_{pr} \\, \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{r\\sigma}\n",
"+\n",
"\\sum_{ \\substack{prqs\\\\\\sigma\\tau} }\n",
"\\frac{(pr|qs)}{2} \\,\n",
"\\hat{a}^\\dagger_{p\\sigma}\n",
"\\hat{a}^\\dagger_{q\\tau}\n",
"\\hat{a}_{s\\tau}\n",
"\\hat{a}_{r\\sigma}\n",
"$$\n",
"\n",
"$\\hat{a}^\\dagger_{p\\sigma}$ / $\\hat{a}_{p\\sigma}$ son los operadores fermiónicos de creación/aniquilación asociados al $p$ -ésimo elemento del conjunto de bases y al espín $\\sigma$. $h_{pr}$ y $(pr|qs)$ son las integrales electrónicas de uno y dos cuerpos. Utilizando pySCF, definiremos la molécula y calcularemos las integrales de uno y dos cuerpos del Hamiltoniano para el conjunto de bases `6-31g`.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cc987c08-7261-4c4c-a06b-609d7003efe9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"converged SCF energy = -108.835236570774\n",
"CASCI E = -109.046671778080 E(CI) = -32.8155692383188 S^2 = 0.0000000\n"
]
}
],
"source": [
"import warnings\n",
"import pyscf\n",
"import pyscf.cc\n",
"import pyscf.mcscf\n",
"\n",
"warnings.filterwarnings(\"ignore\")\n",
"\n",
"# Specify molecule properties\n",
"open_shell = False\n",
"spin_sq = 0\n",
"\n",
"# Build N2 molecule\n",
"mol = pyscf.gto.Mole()\n",
"mol.build(\n",
" atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]], # Two N atoms 1 angstrom apart\n",
" basis=\"6-31g\",\n",
" symmetry=\"Dooh\",\n",
")\n",
"\n",
"# Define active space\n",
"n_frozen = 2\n",
"active_space = range(n_frozen, mol.nao_nr())\n",
"\n",
"# Get molecular integrals\n",
"scf = pyscf.scf.RHF(mol).run()\n",
"num_orbitals = len(active_space)\n",
"n_electrons = int(sum(scf.mo_occ[active_space]))\n",
"num_elec_a = (n_electrons + mol.spin) // 2\n",
"num_elec_b = (n_electrons - mol.spin) // 2\n",
"cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))\n",
"mo = cas.sort_mo(active_space, base=0)\n",
"hcore, nuclear_repulsion_energy = cas.get_h1cas(mo) # hcore: one-body integrals\n",
"eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals) # eri: two-body integrals\n",
"\n",
"# Compute exact energy for comparison\n",
"exact_energy = cas.run().e_tot"
]
},
{
"cell_type": "markdown",
"id": "3aa4e0f0",
"metadata": {},
"source": [
"En esta lección, utilizaremos la transformación Jordan-Wigner (JW) para mapear una función de onda fermiónica a una función de onda qubit de forma que pueda ser preparada utilizando un circuito cuántico. La transformación JW mapea el espacio de Fock de fermiones en M orbitales espaciales en el espacio de Hilbert de 2M qubits, es decir, un orbital espacial se divide en dos *orbitales de espín*, uno asociado con un electrón de espín arriba ( $\\alpha$ ) y otro con espín abajo ( $\\beta$ ). Un orbital de espín puede estar ocupado o desocupado. Normalmente, cuando nos referimos al número de orbitales, estaremos utilizando el número de orbitales *espaciales*. El número de orbitales de espín será el doble. En los circuitos cuánticos, representaremos cada orbital de espín con un qubit. Así, un conjunto de qubits representará los orbitales de espín hacia arriba o $\\alpha$, y otro conjunto representará los orbitales de espín hacia abajo o $\\beta$. Por ejemplo, la molécula $N_2$ para el conjunto de bases `6-31g` tiene orbitales espaciales $16$ (es decir, $16$ $\\alpha$ + $16$ $\\beta$ = $32$ orbitales de espín). Por lo tanto, necesitaremos un circuito cuántico $32$ -qubit (es posible que necesitemos qubits ancilla adicionales, como se discutirá más adelante). Los qubits se miden en base computacional para generar cadenas de bits, que representan configuraciones electrónicas o determinantes (de Slater). A lo largo de esta lección, utilizaremos indistintamente los términos cadenas de bits, configuraciones y determinantes. Las cadenas de bits nos indican la ocupación de electrones en orbitales de espín: un $1$ en una posición de bit significa que el orbital de espín correspondiente está ocupado, mientras que un $0$ significa que el orbital de espín está vacío. Dado que los problemas de estructura electrónica preservan las partículas, sólo debe ocuparse un número fijo de orbitales de espín. La molécula $N_2$ tiene $5$ electrones de espín ascendente ( $\\alpha$ ) y $5$ electrones de espín descendente ( $\\beta$ ). Así, cualquier cadena de bits que represente los orbitales $\\alpha$ y $\\beta$ debe tener cinco $1\\text{s}$ cada uno para la molécula $N_2$.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "f6a89c89-85e1-4c0d-abeb-8ab8b154a7ba",
"metadata": {},
"source": [
"\n",
"\n",
"### 1.1 Circuito cuántico para la generación de muestras: el enfoque LUCJ\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "a464c865-1528-45c2-8ede-325583b15976",
"metadata": {},
"source": [
"En esta lección, utilizaremos el ansatz local unitario acoplado de Jastrow (LUCJ) \\ \\[[1\\]](#references) para la preparación del estado cuántico y el posterior muestreo. En primer lugar, explicaremos los distintos componentes del ansatz UCJ completo y las aproximaciones realizadas en la versión local del mismo. A continuación, utilizando el paquete ffsim, construiremos el ansatz LUCJ y lo optimizaremos utilizando el transpilador Qiskit para su ejecución en hardware.\n",
"\n",
"El ansatz UCJ tiene la siguiente forma (para un producto de $L$ capas o repeticiones del operador UCJ.)\n",
"\n",
"$$\n",
"|\\psi\\rangle = \\prod_{\\mu=1}^{L}{(e^{K^{\\mu}} \\times {e^{iJ^{\\mu}}} \\times {e^{-K^{\\mu}}})} |\\Phi_{0}\\rangle\n",
"$$\n",
"\n",
"donde, $\\vert \\Phi_{0} \\rangle$ es un estado de referencia, típicamente tomado como el estado Hartree-Fock (HF). Como el estado Hartree-Fock se define por tener ocupados los orbitales de menor número, la preparación del estado HF implicará aplicar puertas X para poner a uno los qubits correspondientes a los orbitales ocupados. Por ejemplo, el bloque de preparación del estado HF para 4 orbitales espaciales y 2 up- y 2 down-spin puede tener el siguiente aspecto:\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "0f0bb614",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "4a32d88b",
"metadata": {},
"source": [
"Una sola repetición del operador UCJ ${(e^{K^{(\\mu)}} \\times {e^{iJ^{(\\mu)}}} \\times {e^{-K^{(\\mu)}}})}$ consiste en una evolución diagonal de Coulomb ( $e^{iJ^{(\\mu)}}$ ) intercalada por rotaciones orbitales ( $e^{K^{(\\mu)}}$ y $e^{-K^{(\\mu)}}$ ).\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "963e8386-39b6-40b7-9740-fffcf1573fe6",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "271b1825-bc78-4957-bbf8-21a714e45ed3",
"metadata": {},
"source": [
"Los bloques de rotación orbital funcionan con una sola especie de espín ( $\\alpha$ (up-spin)/ $\\beta$ (down-spin)). Para cada especie de electrón, la rotación orbital consiste en una capa de compuertas de un solo qubit $R_{z}$ seguidas de una secuencia de compuertas de rotación de 2 qubits Given ( $XX + YY$ gates).\n",
"\n",
"Las puertas de 2 qubits actúan sobre orbitales de espín adyacentes (qubits vecinos más cercanos) y, por tanto, pueden implementarse en las QPU de IBM® sin necesidad de puertas SWAP.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "2e8ac1d2-8f04-4591-921b-8ba0174e4ad0",
"metadata": {},
"source": [
"Esquema \n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ebc9b48c-e77d-4af5-b8cf-460daeeadcdb",
"metadata": {},
"source": [
"El $e^{iJ^{(\\mu)}}$, también conocido como operador diagonal de Coulomb, consta de tres bloques. Dos de ellos funcionan en los mismos sectores de espín ( $e^{iJ_{\\alpha \\alpha}^{(\\mu)}}$ y $e^{iJ_{\\beta \\beta}^{(\\mu)}}$ ), y uno funciona entre dos sectores de espín ( $e^{iJ_{\\alpha \\beta}^{(\\mu)}}$ ).\n",
"\n",
"Todos los bloques en $e^{iJ^{(\\mu)}}$ consisten en puertas número-número $U_{nn}(\\phi)$ [\\[1\\]](#references). Una puerta $U_{nn}(\\phi)$ puede descomponerse a su vez en una puerta $R_{ZZ}(\\frac{\\phi}{2})$ seguida de dos puertas $Rz(-\\frac{\\phi}{2})$ de un solo qubit que actúan sobre dos qubits separados.\n",
"\n",
"Los componentes del mismo espín ( $J_{\\alpha \\alpha}$ y $J_{\\beta \\beta}$ ) tienen $U_{nn}$ puertas entre todos los pares posibles de qubits. Sin embargo, como las QPU superconductoras tienen una conectividad restrictiva, los qubits deben intercambiarse para realizar puertas entre qubits no adyacentes.\n",
"\n",
"Por ejemplo, considere el siguiente bloque $e^{iJ_{\\alpha \\alpha}^{(\\mu)}}$ (o $e^{iJ_{\\beta \\beta}^{(\\mu)}}$ ) para orbitales espaciales $N = 4$. Para una conectividad lineal de qubits, las tres últimas puertas no son directamente implementables, ya que funcionan entre qubits no adyacentes (por ejemplo, Q0 y Q2 no están directamente conectados). Por lo tanto, necesitamos puertas SWAP para hacerlas adyacentes (la siguiente figura muestra un ejemplo con $3$ puertas SWAP).\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d3e24a20-1c86-4aea-8300-13bd2ead4e8f",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ef3f1d36-e8ab-46f0-bddb-fa8fb884e531",
"metadata": {},
"source": [
"A continuación, $J_{\\alpha \\beta}$ implementa puertas entre los mismos orbitales indexados de diferentes sectores de espín (por ejemplo, entre $0\\alpha$ y $0\\beta$ ). Del mismo modo, si los qubits no son físicamente adyacentes en una QPU, estas puertas también requerirán SWAPs.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9afe0036-318c-43b9-9b2c-39a808bda82c",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "6219aa8e-8a63-4bfc-b52b-64507292471d",
"metadata": {},
"source": [
"A partir de la discusión anterior, el ansatz UCJ se enfrenta a algunos obstáculos para la ejecución HW, ya que necesita puertas SWAP debido a las interacciones qubit no adyacentes. La variante local del ansatz UCJ, LUCJ, aborda este reto eliminando algunos $U_{nn}$ del operador de Coulomb diagonal.\n",
"\n",
"En los mismos bloques de especies de electrones, $J_{\\alpha \\alpha}$ y $J_{\\beta \\beta}$ ), sólo mantenemos las puertas $U_{nn}$ compatibles con la conectividad de vecino más cercano y eliminamos las puertas entre qubits no adyacentes en la versión LUCJ. La siguiente figura muestra el bloque LUCJ después de eliminar las puertas no adyacentes.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "24f54a55-4a09-4508-997a-96306835c7e6",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "165d0a19-0366-4af8-8ea0-20df351f49f5",
"metadata": {},
"source": [
"A continuación, la versión LUCJ del bloque $J_{\\alpha \\beta}$ que funciona entre diferentes especies de electrones puede adoptar diferentes formas en función de la topología del dispositivo.\n",
"\n",
"También en este caso, la versión LUCJ se deshace de las puertas no compatibles. La siguiente figura muestra variantes del bloque $J_{\\alpha \\beta}$ para diferentes topologías de qubits, incluyendo rejilla, hexagonal, heavy-hex y lineal.\n",
"\n",
"* **Cuadrado** : podemos tener puertas $U_{nn}$ entre todos los orbitales $\\alpha$ y $\\beta$ sin ningún SWAP, y por lo tanto, no necesitamos eliminar ninguna puerta $U_{nn}$.\n",
"* **Heavy-hex** : Las interacciones $\\alpha$ - $\\beta$ se mantienen entre cada $4$ -ésimo orbital de espín indexado (como el 0º, 4º y 8º) y están mediadas por *ancilla*, es decir, necesitamos qubits ancilla entre las cadenas lineales que representan los orbitales $\\alpha$ y $\\beta$. Este acuerdo necesita un número limitado de SWAPs.\n",
"* **Hexagonal** : Todos los demás orbitales, como los orbitales 0, 2 y 4 indexados, se convierten en vecinos más cercanos cuando $\\alpha$ y $\\beta$ se disponen en dos cadenas lineales adyacentes.\n",
"* **Lineal** : Sólo se conectan un orbital $\\alpha$ y otro $\\beta$, lo que significa que el bloque $J_{\\alpha \\beta}$ sólo tendrá una puerta.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8ec1e433-bc00-42bc-bdbb-288c56c32f9d",
"metadata": {},
"source": [
"Diagramas de \n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "3e070615-2b4a-4866-8c6c-d2d038447f45",
"metadata": {},
"source": [
"Aunque la eliminación de puertas del ansatz UCJ para construir la versión LUCJ lo hace más compatible con HW, el ansatz pierde algo de expresividad. Por lo tanto, pueden ser necesarias más repeticiones ( $L$ ) del operador UCJ modificado cuando se utiliza el ansatz LUCJ.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "04367dac",
"metadata": {},
"source": [
"\n",
"\n",
"### 1.2 Inicialización del ansatz LUCJ\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "f5eae55f",
"metadata": {},
"source": [
"El LUCJ es un ansatz parametrizado, y necesitamos inicializar los parámetros antes de la ejecución del hardware. Una forma de inicializar el ansatz es utilizando las amplitudes `t1` y `t2` del método clásico de clústeres simples y dobles acoplados (CCSD), donde las amplitudes `t1` son el coeficiente de los operadores de excitación simple y las amplitudes `t2` son para los operadores de excitación doble.\n",
"\n",
"Obsérvese que aunque la inicialización del ansatz LUCJ con `t1` y `t2` amplitudes genera resultados decentes, los parámetros del ansatz pueden necesitar una mayor optimización.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "1835e74e-3354-425f-8596-c574f03e7a6e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"E(CCSD) = -109.0398256929733 E_corr = -0.20458912219883\n"
]
}
],
"source": [
"# Get CCSD t2 amplitudes for initializing the ansatz\n",
"ccsd = pyscf.cc.CCSD(\n",
" scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]\n",
")\n",
"ccsd.run()\n",
"\n",
"t1 = ccsd.t1\n",
"t2 = ccsd.t2"
]
},
{
"cell_type": "markdown",
"id": "9f0842fa",
"metadata": {},
"source": [
"\n",
"\n",
"### 1.3 Construcción del enfoque LUCJ utilizando `ffsim`\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "e5f7f7e8-30e3-40e8-9b85-bf057b35b766",
"metadata": {},
"source": [
"Utilizaremos el paquete [ffsim](https://github.com/qiskit-community/ffsim/tree/main) para crear e inicializar el ansatz con `t1` y `t2` amplitudes calculadas anteriormente. Dado que nuestra molécula tiene un estado Hartree-Fock de cáscara cerrada, utilizaremos la variante de espín equilibrado del ansatz UCJ, [UCJOpSpinBalanced](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced).\n",
"\n",
"Como el hardware de IBM tiene una topología heavy-hex, adoptaremos el patrón *zig-zag* utilizado en [\\[1\\]](#references) y explicado anteriormente para las interacciones qubit. En este patrón, los orbitales (qubits) con el mismo espín están conectados con una topología de línea (círculos rojos y azules). Debido a la topología hexagonal pesada, los orbitales para diferentes espines tienen conexiones entre cada 4 orbitales, es decir, el 0, el 4, el 8, etc. (círculos morados).\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "6799421a-4425-404a-9994-47215957054d",
"metadata": {},
"outputs": [],
"source": [
"import ffsim\n",
"from qiskit import QuantumCircuit, QuantumRegister\n",
"\n",
"n_reps = 2\n",
"alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]\n",
"alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]\n",
"\n",
"ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
" t2=t2,\n",
" t1=t1,\n",
" n_reps=n_reps,\n",
" interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),\n",
")\n",
"\n",
"nelec = (num_elec_a, num_elec_b)\n",
"\n",
"# create an empty quantum circuit\n",
"qubits = QuantumRegister(2 * num_orbitals, name=\"q\")\n",
"circuit = QuantumCircuit(qubits)\n",
"\n",
"# prepare Hartree-Fock state as the reference state and append it to the quantum circuit\n",
"circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)\n",
"\n",
"# apply the UCJ operator to the reference state\n",
"circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n",
"circuit.measure_all()\n",
"# circuit.decompose().draw(\"mpl\", scale=0.5, fold=-1)"
]
},
{
"cell_type": "markdown",
"id": "3b2de316",
"metadata": {},
"source": [
"El ansatz LUCJ con capas repetidas puede optimizarse fusionando algunos bloques adyacentes. Consideremos un caso para `n_reps=2`. Los dos bloques de rotación orbital del centro pueden fusionarse en un único bloque de rotación orbital. El paquete `ffsim` dispone de un gestor de pasadas denominado `ffsim.qiskit.PRE_INIT` para optimizar el circuito fusionando dichos bloques adyacentes.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "7cb99cd9",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "e9ad460d",
"metadata": {},
"source": [
"\n",
"\n",
"## 2. Optimizar para el hardware de destino\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "cecf2994",
"metadata": {},
"source": [
"En primer lugar, buscamos un backend de nuestra elección. Optimizaremos nuestro circuito para el backend, y luego ejecutaremos el circuito optimizado en el mismo backend para generar muestras para el subespacio.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4c8dbbef-378e-4ae3-9290-bde58b72024d",
"metadata": {},
"outputs": [],
"source": [
"from qiskit_ibm_runtime import QiskitRuntimeService\n",
"\n",
"service = QiskitRuntimeService()\n",
"# Use the least-busy backend or specify a quantum computer using the syntax commented out below.\n",
"backend = service.least_busy(operational=True, simulator=False)\n",
"# backend = service.backend(\"ibm_brisbane\")"
]
},
{
"cell_type": "markdown",
"id": "1b8ba388-6904-4725-ad97-2a31c0adaa07",
"metadata": {},
"source": [
"A continuación, recomendamos los siguientes pasos para optimizar el ansatz y hacerlo compatible con el hardware.\n",
"\n",
"* Seleccione qubits físicos (`initial_layout`) del hardware de destino que se adhieran al patrón en zig-zag (dos cadenas lineales con un qubit ancilla entre ellas) descrito anteriormente. La disposición de los qubits según este patrón da lugar a un circuito eficiente compatible con hardware y con menos puertas.\n",
"* Genere un gestor de pases por etapas utilizando la función [`generate_preset_pass_manager`](/docs/api/qiskit/qiskit.transpiler.generate_preset_pass_manager) de Qiskit con su elección de `backend` y `initial_layout`.\n",
"* Establezca la etapa `pre_init` de su gestor de pases escalonados en `ffsim.qiskit.PRE_INIT`. `ffsim.qiskit.PRE_INIT` incluye pases de transpilador Qiskit que descomponen las puertas en rotaciones orbitales y luego fusionan las rotaciones orbitales, lo que resulta en menos puertas en el circuito final.\n",
"* Ejecuta el gestor de pases en tu circuito.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "ecd2675c-a302-43fb-80f4-d658d56360d5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Gate counts (w/o pre-init passes): OrderedDict({'rz': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})\n",
"Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})\n"
]
}
],
"source": [
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
"\n",
"spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]\n",
"spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]\n",
"\n",
"initial_layout = spin_a_layout + spin_b_layout\n",
"\n",
"pass_manager = generate_preset_pass_manager(\n",
" optimization_level=3, backend=backend, initial_layout=initial_layout\n",
")\n",
"\n",
"# without PRE_INIT passes\n",
"isa_circuit = pass_manager.run(circuit)\n",
"print(f\"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}\")\n",
"\n",
"# with PRE_INIT passes\n",
"# We will use the circuit generated by this pass manager for hardware execution\n",
"pass_manager.pre_init = ffsim.qiskit.PRE_INIT\n",
"isa_circuit = pass_manager.run(circuit)\n",
"print(f\"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}\")"
]
},
{
"cell_type": "markdown",
"id": "72a780f9",
"metadata": {},
"source": [
"\n",
"\n",
"## 3. Ejecutar en el hardware de destino\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "4cd9164c-697a-4aa6-b71f-86720e3d5b66",
"metadata": {},
"source": [
"Tras optimizar el circuito para su ejecución en hardware, estamos listos para ejecutarlo en el hardware de destino y recopilar muestras para la estimación de la energía del estado fundamental. Como solo tenemos un circuito, utilizaremos el y `qiskit-ibm-runtime`[Modo de ejecución de tareas](/docs/guides/execution-modes) ejecutaremos nuestro circuito.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "9f542448-5767-4e12-85c8-dd8584545dbb",
"metadata": {},
"outputs": [],
"source": [
"from qiskit_ibm_runtime import SamplerV2 as Sampler\n",
"\n",
"sampler = Sampler(mode=backend)\n",
"sampler.options.dynamical_decoupling.enable = True\n",
"\n",
"job = sampler.run([isa_circuit], shots=10_000) # Takes approximately 5sec of QPU time"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8c5eb3e9-2a5a-423b-ac18-7bba7f69f18f",
"metadata": {},
"outputs": [],
"source": [
"# Run cell after IQX job completion\n",
"primitive_result = job.result()\n",
"pub_result = primitive_result[0]\n",
"counts = pub_result.data.meas.get_counts()"
]
},
{
"cell_type": "markdown",
"id": "e3ebf77a-8efd-43d6-bd19-e26423921c6f",
"metadata": {},
"source": [
"\n",
"\n",
"## 4. Resultados posteriores al procesamiento\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "1e7f0ecd",
"metadata": {},
"source": [
"La parte de posprocesamiento del flujo de trabajo de SQD puede resumirse mediante el siguiente diagrama.\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8c51a527",
"metadata": {},
"source": [
"El muestreo del ansatz LUCJ en la base de cálculo genera un conjunto de configuraciones ruidosas $\\tilde{\\mathcal{\\chi}}$, que se utilizan en la rutina de posprocesamiento. Se trata de un método denominado *recuperación de configuraciones* (que se explicará más adelante) para corregir probabilísticamente las configuraciones con números de electrones incorrectos. A continuación, las configuraciones que sólo tienen los números de electrones correctos $\\tilde{\\mathcal{\\chi}}_{R}$ se submuestrean y se distribuyen en varios lotes en función de la frecuencia de aparición de cada configuración única. Cada lote de muestras define un subespacio ( $\\mathcal{S^{(k)}}$ ). A continuación, el Hamiltoniano molecular, $H$, se proyecta en subespacios:\n",
"\n",
"$$\n",
"H_{\\mathcal{S}^{(k)}} = P_{\\mathcal{S}^{(k)}} H _{\\mathcal{S}^{(k)}} \\text{ with } P_{\\mathcal{S}^{(k)}} = \\sum_{x \\in \\mathcal{S}^{(k)}} \\vert x \\rangle \\langle x \\vert\n",
"$$\n",
"\n",
"Cada Hamiltoniano proyectado $H_{\\mathcal{S}^{(k)}}$ se introduce en un Eigensolver, donde se diagonaliza para calcular los valores propios y los vectores propios para reconstruir un estado propio. En esta lección, proyectamos y diagonalizamos el Hamiltoniano utilizando el paquete `qiskit-addon-sqd` que utiliza el método de Davidson de PySCF para la diagonalización.\n",
"\n",
"$$\n",
"H_{\\mathcal{S}^{(k)}} \\vert \\psi^{(k)} \\rangle = E^{(k)} \\vert \\psi^{(k)} \\rangle\n",
"$$\n",
"\n",
"A continuación, recogemos el valor propio más bajo (energía) de los lotes y también calculamos la ocupación orbital media, $\\text{n}$. La información sobre la ocupación media se utiliza en el paso de recuperación de la configuración para corregir probabilísticamente las configuraciones con ruido.\n",
"\n",
"A continuación, explicamos en detalle el bucle de recuperación de la configuración autoconsistente y mostramos ejemplos concretos de código para implementar los pasos mencionados con el fin de estimar la energía del estado fundamental del Hamiltoniano de $N_2$.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "39997b5d-8bd3-4bc8-9ade-11fa5cfb34f9",
"metadata": {},
"source": [
"\n",
"\n",
"### 4.1 Recuperación de la configuración: descripción general\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "b05881ac-80ce-47a6-ac28-c1237f85bc0a",
"metadata": {},
"source": [
"Cada bit de una cadena de bits (determinante de Slater) representa un orbital de espín. La mitad derecha de una cadena de bits representa orbitales de espín ascendente y la mitad izquierda orbitales de espín descendente. Un `1` significa que el orbital está ocupado por un electrón, y un `0` significa que el orbital está vacío. Conocemos a priori el número correcto de partículas (electrón de espín ascendente y electrón de espín descendente). Supongamos que tenemos un determinante $x$ con $N_x$ electrones (es decir, hay $N_x$ números de $1$ s en la cadena de bits) en él. El número correcto de partículas es $N$. Si $N_x \\neq N$, entonces sabemos que la cadena de bits está corrompida por ruido. La rutina de configuración autoconsistente intenta corregir la cadena de bits volteando probabilísticamente $|N_x - N|$ bits aprovechando la información de ocupación orbital media. La ocupación orbital media ( $n$ ) nos indica la probabilidad de que un orbital esté ocupado por un electrón. Si $N_x < N$, tenemos menos electrones y necesitamos voltear algunos $0$ s a $1$ s y viceversa.\n",
"\n",
"La probabilidad de flipping puede ser $|x[i] - avg\\_occupancy[i]|$ para `i`-th spin orbital. En [\\[2\\]](#references), los autores utilizaron una probabilidad ponderada de volteo utilizando la función ReLU modificada.\n",
"\n",
"$$\n",
"\\begin{align}\n",
" w(y) = \\begin{cases}\n",
"\n",
" \\delta \\frac{y}{h} & \\text{if } y \\leq h\\\\ \\nonumber\n",
"\n",
" \\delta + (1 - \\delta) \\frac{y - h}{1 - h} & \\text{if } y > h\n",
"\n",
"\\end{cases}\n",
"\\end{align}\n",
"$$\n",
"\n",
"Aquí $h$ define la ubicación de la \"esquina\" de la función ReLU, y el parámetro $\\delta$ define el valor de la función ReLU en la esquina. Para $\\delta = 0$, $w$ se convierte en la función ReLU verdadera, y para $\\delta >0$, se convierte en ReLU\\* modificada\\*. En el artículo, los autores utilizaron $\\delta = 0.01$ y $h =$ número de partículas alfa (o beta)/número de orbitales de espín alfa (o beta) $= N/M$ (factor de llenado).\n",
"\n",
"La ocupación orbital media ( $n$ ) no se conoce a priori. La primera iteración de la estimación del estado básico comienza con configuraciones que sólo tienen números de partículas correctos en ambas especies de espín. Después de la primera iteración, tenemos una estimación del estado base, y utilizando la estimación, podemos construir la primera conjetura de $n$. Esta conjetura de $n$ se utiliza para recuperar configuraciones, ejecutar la siguiente iteración de estimación del estado base, y auto-consistentemente refinar la conjetura de $n$. El proceso se repite hasta que se cumple un criterio de parada.\n",
"\n",
"Considere el siguiente ejemplo para $N = 2$ y $x = |1000\\rangle$ ( $N_x = 1$ ). Tenemos que dar la vuelta a uno de los 0s a 1 para corregirlo para los números de partículas, y las opciones son `1100`, `1010`, y `1001`. En función de la probabilidad de volteo, se seleccionará una de las opciones como *configuración recuperada* (o la cadena de bits con el número correcto de partículas).\n",
"\n",
"Supongamos que en la primera iteración ejecutamos dos lotes, y los estados básicos estimados a partir de ellos son:\n",
"\n",
"$$\n",
"\\begin{align}\\nonumber\n",
" \\text{Batch0: } \\vert \\psi \\rangle &= 0.8 \\times \\vert 1001 \\rangle + 0.6 \\times \\vert 0110 \\rangle \\\\ \\nonumber\n",
" \\text{Batch1: } \\vert \\psi \\rangle &= \\frac{1}{\\sqrt{3}} \\left( \\vert 1001 \\rangle + \\vert 0101 \\rangle + \\vert 0110 \\rangle \\right) \\nonumber\n",
"\\end{align}\n",
"$$\n",
"\n",
"Utilizando los estados base computacionales y sus amplitudes, podemos calcular la probabilidad de las ocupaciones de electrones (en breve *ocupaciones* ) por orbital de espín (qubit) (nótese que probabilidad = |amplitud| $^2$ ). A continuación tabulamos las ocupaciones por qubit para cada cadena de bits que aparece en el estado básico estimado y calculamos la ocupación orbital total para un lote. Tenga en cuenta que, según la convención de ordenación de Qiskit, el bit situado más a la derecha representa qubit-0 ( Q0 ), y el situado más a la izquierda representa Q3.\n",
"\n",
"Ocupación ( Batch0 ):\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :--------------: | :------: | :------: | :------: | :------: |\n",
"| 1001 | 0.64 | 0.0 | 0.0 | 0.64 |\n",
"| 0110 | 0.0 | 0.36 | 0.36 | 0.0 |\n",
"| **n** *(Batch0)* | **0.64** | **0.36** | **0.36** | **0.64** |\n",
"\n",
"Ocupación ( Batch1 )\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :--------------: | :------: | :------: | :------: | :------: |\n",
"| 1001 | 0.33 | 0.00 | 0.00 | 0.33 |\n",
"| 0101 | 0.0 | 0.33 | 0.00 | 0.33 |\n",
"| 0110 | 0.0 | 0.33 | 0.33 | 0.00 |\n",
"| **n** *(Batch1)* | **0.33** | **0.66** | **0.33** | **0.66** |\n",
"\n",
"Ocupación (media de los lotes)\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :--------------: | :------: | :------: | :------: | :------: |\n",
"| **n** *(Batch0)* | 0.64 | 0.36 | 0.36 | 0.64 |\n",
"| **n** *(Batch1)* | 0.33 | 0.66 | 0.33 | 0.66 |\n",
"| **n** *(media)* | **0.49** | **0.51** | **0.35** | **0.65** |\n",
"\n",
"Utilizando la ocupación orbital media calculada anteriormente, podemos hallar las probabilidades de volteo para diferentes orbitales en la configuración $x = \\vert 1000 \\rangle$. Como el orbital representado por Q3 ya está ocupado y no es necesario voltearlo, fijamos su p(flip) en $0$. Para el resto de orbitales, que están desocupados, la probabilidad de volteo es $\\vert x[i] - \\text{n}[i] \\vert$ para cada uno. Junto con p(flip), también calculamos el peso de probabilidad asociado a flipping utilizando la función ReLU modificada descrita anteriormente.\n",
"\n",
"Probabilidad de flip ( $x = \\vert 1000 \\rangle$, $\\delta = 0.01$, $h = N/M = 2/4 = 0.50$ )\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :------------------------------------------: | :-: | :--: | :---: | :--: |\n",
"| p(flip) ( $\\vert x[i] - \\text{n}[i] \\vert$ ) | 0 | 0.51 | 0.35 | 0.65 |\n",
"| w(p(flip)) | 0 | 0.03 | 0.007 | 0.31 |\n",
"\n",
"Finalmente, utilizando las probabilidades ponderadas anteriores, podemos voltear uno de los orbitales desocupados Q2, Q1, y Q0. Basándose en los valores anteriores, lo más probable es que se invierta Q0, y una posible configuración recuperada puede ser $\\vert \\text{1001} \\rangle$.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "c1940235",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "5c614290",
"metadata": {},
"source": [
"El proceso completo de recuperación de configuraciones autoconsistentes puede resumirse como sigue:\n",
"\n",
"**Primera iteración:** Supongamos que las cadenas de bits (configuraciones o determinantes de Slater) generadas por el ordenador cuántico forman un conjunto $\\widetilde{\\chi}$, que incluye tanto configuraciones con número correcto ( $\\widetilde{\\chi}_{correct}$ ) como incorrecto ( $\\widetilde{\\chi}_{incorrect}$ ) de partículas en cada sector de espín.\n",
"\n",
"1. Las configuraciones de ( $\\widetilde{\\chi}_{correct}$ ) se muestrean aleatoriamente para crear lotes $(\\mathcal{S}^{(1)}, \\cdots, \\mathcal{S}^{(K)})$ de vectores para la proyección del subespacio. El número de lotes y las muestras de cada lote son parámetros definidos por el usuario. Cuanto mayor sea el número de muestras de cada lote, mayor será la dimensión del subespacio y más exigente será la diagonalización desde el punto de vista informático. Por otra parte, un número demasiado pequeño de muestras puede pasar por alto los vectores de apoyo del estado básico y dar lugar a una estimación incorrecta.\n",
"2. Ejecute el solucionador de estados propios (es decir, la proyección sobre el subespacio y la diagonalización) en los lotes y obtenga los estados propios aproximados. $|\\psi^{(1)}\\rangle, \\cdots, |\\psi^{(K)}\\rangle$.\n",
"3. A partir de los estados propios aproximados, construya la primera conjetura para $n$.\n",
"\n",
"**Iteraciones posteriores:**\n",
"\n",
"1. Utilizando $n$ corregimos las configuraciones con número de partícula erróneo en $\\widetilde{\\chi}_{incorrect}$. Supongamos que las llamamos $\\widetilde{\\chi}_{correct\\_new}$. Entonces, $\\widetilde{\\chi}_{recovered} (\\widetilde{\\chi}_{R}) = \\widetilde{\\chi}_{correct} \\cup \\widetilde{\\chi}_{correct\\_new}$ forma el nuevo conjunto de configuraciones con número de partículas correcto.\n",
"2. $\\widetilde{\\chi}_{R}$ se muestrea para crear los lotes $\\mathcal{S}^{(1)}, \\cdots, \\mathcal{S}^{(K)}$.\n",
"3. El solucionador de estados propios se ejecuta con nuevos lotes y genera nuevas estimaciones de los estados básicos $|\\psi^{(1)}\\rangle, \\cdots, |\\psi^{(K)}\\rangle$.\n",
"4. A partir de los estados propios aproximados se construyen conjeturas refinadas para $n$.\n",
"5. Si no se cumple el criterio de parada, vuelva al paso `2.1`.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "5eea548a",
"metadata": {},
"source": [
"\n",
"\n",
"### 4.2 Estimación del estado fundamental\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "31dfe00e",
"metadata": {},
"source": [
"En primer lugar, transformaremos los recuentos en una matriz de cadenas de bits y una matriz de probabilidades para su postprocesamiento.\n",
"\n",
"Cada fila de la matriz representa una única cadena de bits. Dado que los qubits se indexan desde la derecha de una cadena de bits en Qiskit, la columna `0` representa el qubit `N-1`, y la columna `N-1` representa el qubit `0`, donde `N` es el número de qubits.\n",
"\n",
"Los orbitales alfa se representan en el rango de índice de columna `(N, N/2]` (mitad derecha), y los orbitales beta en el rango de columna `(N/2, 0]` (mitad izquierda).\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "71550274",
"metadata": {},
"outputs": [],
"source": [
"from qiskit_addon_sqd.counts import counts_to_arrays\n",
"\n",
"# Convert counts into bitstring and probability arrays\n",
"bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)"
]
},
{
"cell_type": "markdown",
"id": "acbdbdc5",
"metadata": {},
"source": [
"Hay algunas opciones controladas por el usuario que son importantes para esta técnica:\n",
"\n",
"* `iterations`: Número de iteraciones de recuperación de configuración autoconsistente\n",
"* `n_batches`: Número de lotes de configuraciones utilizados por las diferentes llamadas al solucionador de estados propios\n",
"* `samples_per_batch`: Número de configuraciones únicas a incluir en cada lote\n",
"* `max_davidson_cycles`: Número máximo de ciclos Davidson ejecutados por cada eigensolver\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a2d22e0b-8f51-42a0-858f-ad0297cb0bae",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting configuration recovery iteration 0\n",
" Batch 0 subspace dimension: 21609\n",
" Batch 1 subspace dimension: 21609\n",
" Batch 2 subspace dimension: 21609\n",
" Batch 3 subspace dimension: 21609\n",
" Batch 4 subspace dimension: 21609\n",
"Starting configuration recovery iteration 1\n",
" Batch 0 subspace dimension: 609961\n",
" Batch 1 subspace dimension: 616225\n",
" Batch 2 subspace dimension: 627264\n",
" Batch 3 subspace dimension: 633616\n",
" Batch 4 subspace dimension: 624100\n",
"Starting configuration recovery iteration 2\n",
" Batch 0 subspace dimension: 564001\n",
" Batch 1 subspace dimension: 605284\n",
" Batch 2 subspace dimension: 582169\n",
" Batch 3 subspace dimension: 559504\n",
" Batch 4 subspace dimension: 591361\n",
"Starting configuration recovery iteration 3\n",
" Batch 0 subspace dimension: 550564\n",
" Batch 1 subspace dimension: 549081\n",
" Batch 2 subspace dimension: 531441\n",
" Batch 3 subspace dimension: 527076\n",
" Batch 4 subspace dimension: 531441\n",
"Starting configuration recovery iteration 4\n",
" Batch 0 subspace dimension: 544644\n",
" Batch 1 subspace dimension: 580644\n",
" Batch 2 subspace dimension: 527076\n",
" Batch 3 subspace dimension: 531441\n",
" Batch 4 subspace dimension: 537289\n"
]
}
],
"source": [
"import numpy as np\n",
"from qiskit_addon_sqd.configuration_recovery import recover_configurations\n",
"from qiskit_addon_sqd.fermion import (\n",
" bitstring_matrix_to_ci_strs,\n",
" solve_fermion,\n",
")\n",
"from qiskit_addon_sqd.subsampling import postselect_and_subsample\n",
"\n",
"rng = np.random.default_rng(24)\n",
"# SQD options\n",
"iterations = 5\n",
"\n",
"# Eigenstate solver options\n",
"n_batches = 5\n",
"samples_per_batch = 500\n",
"max_davidson_cycles = 300\n",
"\n",
"# Self-consistent configuration recovery loop\n",
"e_hist = np.zeros((iterations, n_batches)) # energy history\n",
"s_hist = np.zeros((iterations, n_batches)) # spin history\n",
"occupancy_hist = []\n",
"avg_occupancy = None\n",
"for i in range(iterations):\n",
" print(f\"Starting configuration recovery iteration {i}\")\n",
" # On the first iteration, we have no orbital occupancy information from the\n",
" # solver, so we begin with the full set of noisy configurations.\n",
" if avg_occupancy is None:\n",
" bs_mat_tmp = bitstring_matrix_full\n",
" probs_arr_tmp = probs_arr_full\n",
"\n",
" # If we have average orbital occupancy information, we use it to refine\n",
" # the full set of noisy configurations.\n",
" else:\n",
" bs_mat_tmp, probs_arr_tmp = recover_configurations(\n",
" bitstring_matrix_full,\n",
" probs_arr_full,\n",
" avg_occupancy,\n",
" num_elec_a,\n",
" num_elec_b,\n",
" rand_seed=rng,\n",
" )\n",
"\n",
" # Create batches of subsamples. We postselect here to remove configurations\n",
" # with incorrect hamming weight during iteration 0, since no config recovery was performed.\n",
" batches = postselect_and_subsample(\n",
" bs_mat_tmp,\n",
" probs_arr_tmp,\n",
" hamming_right=num_elec_a,\n",
" hamming_left=num_elec_b,\n",
" samples_per_batch=samples_per_batch,\n",
" num_batches=n_batches,\n",
" rand_seed=rng,\n",
" )\n",
"\n",
" # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.\n",
" e_tmp = np.zeros(n_batches)\n",
" s_tmp = np.zeros(n_batches)\n",
" occs_tmp = []\n",
" coeffs = []\n",
" for j in range(n_batches):\n",
" strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])\n",
" print(f\" Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}\")\n",
" energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(\n",
" batches[j],\n",
" hcore,\n",
" eri,\n",
" open_shell=open_shell,\n",
" spin_sq=spin_sq,\n",
" max_cycle=max_davidson_cycles,\n",
" )\n",
" energy_sci += nuclear_repulsion_energy\n",
" e_tmp[j] = energy_sci\n",
" s_tmp[j] = spin\n",
" occs_tmp.append(avg_occs)\n",
" coeffs.append(coeffs_sci)\n",
"\n",
" # Combine batch results\n",
" avg_occupancy = tuple(np.mean(occs_tmp, axis=0))\n",
"\n",
" # Track optimization history\n",
" e_hist[i, :] = e_tmp\n",
" s_hist[i, :] = s_tmp\n",
" occupancy_hist.append(avg_occupancy)"
]
},
{
"cell_type": "markdown",
"id": "e4f6ec2c-032d-4e18-8ea4-cf867cba5054",
"metadata": {},
"source": [
"\n",
"\n",
"### 4.3 Discusión de los resultados\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "592eabb7-18f3-4622-8710-bfc43ad6cdec",
"metadata": {},
"source": [
"El primer gráfico muestra que, tras unas pocas iteraciones, estimamos la energía del estado básico con un margen de \\~24 mH (la precisión química suele aceptarse en 1 kcal/mol $\\approx$ 1.6 mH ). El segundo gráfico muestra la ocupación media de cada orbital espacial tras la iteración final. Podemos ver que tanto los electrones de espín arriba como los de espín abajo ocupan los cinco primeros orbitales con alta probabilidad en nuestras soluciones.\n",
"\n",
"Aunque la energía estimada del estado básico es decente, no está dentro del límite de precisión química ( $\\pm \\approx 1.6$ mH ). Este desfase puede atribuirse a la pequeña dimensión del subespacio que utilizamos anteriormente para la proyección y la diagonalización. Como utilizamos `samples_per_batch=500`, el subespacio está abarcado por el máximo de vectores de $500$, que son los vectores que faltan en el soporte del estado base. Aumentar el parámetro `samples_per_batch` debería mejorar la precisión a costa de más recursos informáticos clásicos y tiempo de ejecución.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "61959e69-a182-4636-abcb-a32349fc9076",
"metadata": {},
"outputs": [],
"source": [
"# Data for energies plot\n",
"x1 = range(iterations)\n",
"min_e = [np.min(e) for e in e_hist]\n",
"e_diff = [abs(e - exact_energy) for e in min_e]\n",
"yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]\n",
"\n",
"# Chemical accuracy (+/- 1 milli-Hartree)\n",
"chem_accuracy = 0.001\n",
"\n",
"# Data for avg spatial orbital occupancy\n",
"y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]\n",
"x2 = range(len(y2))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "8cd90034-6ef3-41bd-a847-c115cade82f7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact energy: -109.04667 Ha\n",
"SQD energy: -109.02234 Ha\n",
"Absolute error: 0.02434 Ha\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n",
"\n",
"# Plot energies\n",
"axs[0].plot(x1, e_diff, label=\"energy error\", marker=\"o\")\n",
"axs[0].set_xticks(x1)\n",
"axs[0].set_xticklabels(x1)\n",
"axs[0].set_yticks(yt1)\n",
"axs[0].set_yticklabels(yt1)\n",
"axs[0].set_yscale(\"log\")\n",
"axs[0].set_ylim(1e-6)\n",
"axs[0].axhline(\n",
" y=chem_accuracy, color=\"#BF5700\", linestyle=\"--\", label=\"chemical accuracy\"\n",
")\n",
"axs[0].set_title(\"Approximated Ground State Energy Error vs SQD Iterations\")\n",
"axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n",
"axs[0].set_ylabel(\"Energy Error (Ha)\", fontdict={\"fontsize\": 12})\n",
"axs[0].legend()\n",
"\n",
"# Plot orbital occupancy\n",
"axs[1].bar(x2, y2, width=0.8)\n",
"axs[1].set_xticks(x2)\n",
"axs[1].set_xticklabels(x2)\n",
"axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n",
"axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n",
"axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n",
"\n",
"print(f\"Exact energy: {exact_energy:.5f} Ha\")\n",
"print(f\"SQD energy: {min_e[-1]:.5f} Ha\")\n",
"print(f\"Absolute error: {e_diff[-1]:.5f} Ha\")\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "e1530472",
"metadata": {},
"source": [
"\n",
"\n",
"#### Ejercicio para el lector\n",
"\n",
"Aumente progresivamente el parámetro `samples_per_batch` (por ejemplo, de $1000$ a $10000$ a un paso de $1000$; permitido mi memoria de su ordenador) y compare las energías estimadas del estado básico.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d2b2241f",
"metadata": {},
"source": [
"\n",
"\n",
"## Referencias\n",
"\n",
"\\[1] M. Motta et al., \"Bridging physical intuition and hardware efficiency for correlated electronic states: the local unitary cluster Jastrow ansatz for electronic structure\" (2023). [Química. Sci., 2023, 14, 11213](https://pubs.rsc.org/en/content/articlehtml/2023/sc/d3sc02516k).\n",
"\n",
"\\[2] J. Robledo-Moreno et al., \"Química más allá de las soluciones exactas en un superordenador centrado en la cuántica\" (2024). [arXiv:quant-ph/2405.05068](https://arxiv.org/abs/2405.05068).\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"
}
},
"nbformat": 4,
"nbformat_minor": 5
}