{
"cells": [
{
"cell_type": "markdown",
"id": "frontmatter",
"metadata": {},
"source": [
"---\n",
"title: \"Mejorar una estimación del SQD mediante la optimización orbital\"\n",
"description: \"Mejora de una estimación SQD mediante optimización orbital para la última versión de la diagonalización cuántica basada en muestras (SQD)\"\n",
"---\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "bb5a576d",
"metadata": {},
"source": [
"\n",
"\n",
"# Mejorar una estimación del SQD mediante la optimización orbital\n",
"\n",
"La diagonalización cuántica basada en muestras (SQD) aproxima la energía del estado fundamental mediante\n",
"la diagonalización del hamiltoniano en un subespacio fijo de configuraciones electrónicas. Esa\n",
"estimación depende de la base orbital en la que se expresa el hamiltoniano, y\n",
"*la optimización orbital* (OO) aprovecha esta libertad para reducir la energía sin ampliar\n",
"el subespacio.\n",
"\n",
"Esta guía ejecuta SQD sobre una molécula « $N_2$ » y, a continuación, mejora el resultado mediante la optimización orbital,\n",
"utilizando [`ffsim`](https://qiskit-community.github.io/ffsim/) para representar\n",
"el hamiltoniano y hallar la rotación orbital que minimiza la energía.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "77ab953b",
"metadata": {},
"source": [
"\n",
"\n",
"## Ejecutar SQD\n",
"\n",
"Construimos las integrales moleculares para $N_2$ en la base de orbitales moleculares (MO), generamos\n",
"muestras aleatorias uniformes y ejecutamos SQD para obtener una aproximación del estado fundamental.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "b8d5618e",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-16T01:36:48.334261Z",
"iopub.status.busy": "2026-07-16T01:36:48.334063Z",
"iopub.status.idle": "2026-07-16T01:37:42.585853Z",
"shell.execute_reply": "2026-07-16T01:37:42.584303Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"converged SCF energy = -108.835236570775\n",
"CASCI E = -109.046671778080 E(CI) = -32.8155692383187 S^2 = 0.0000000\n"
]
}
],
"source": [
"import numpy as np\n",
"import pyscf\n",
"import pyscf.cc\n",
"import pyscf.mcscf\n",
"from qiskit_addon_sqd.counts import generate_bit_array_uniform\n",
"from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian\n",
"\n",
"# Specify molecule properties\n",
"num_orbitals = 16\n",
"num_elec_a = num_elec_b = 5\n",
"spin_sq = 0\n",
"\n",
"# Build N2 molecule\n",
"mol = pyscf.gto.Mole()\n",
"mol.build(\n",
" atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]],\n",
" basis=\"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)\n",
"eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)\n",
"\n",
"# Compute exact energy\n",
"exact_energy = cas.run().e_tot\n",
"\n",
"# Create a seed to control randomness throughout this workflow\n",
"rng = np.random.default_rng(24)\n",
"\n",
"\n",
"# Generate random samples\n",
"bit_array = generate_bit_array_uniform(\n",
" 10_000, num_orbitals * 2, rand_seed=rng\n",
")\n",
"\n",
"# Run SQD\n",
"result = diagonalize_fermionic_hamiltonian(\n",
" hcore,\n",
" eri,\n",
" bit_array,\n",
" samples_per_batch=100,\n",
" norb=num_orbitals,\n",
" nelec=(num_elec_a, num_elec_b),\n",
" num_batches=1,\n",
" max_iterations=5,\n",
" symmetrize_spin=True,\n",
" seed=rng,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "0ca70028",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-16T01:37:42.590790Z",
"iopub.status.busy": "2026-07-16T01:37:42.589477Z",
"iopub.status.idle": "2026-07-16T01:37:42.595563Z",
"shell.execute_reply": "2026-07-16T01:37:42.595133Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact energy: -109.04667178\n",
"SQD energy: -108.98469255\n"
]
}
],
"source": [
"sqd_energy = result.energy + nuclear_repulsion_energy\n",
"print(f\"Exact energy: {exact_energy:.8f}\")\n",
"print(f\"SQD energy: {sqd_energy:.8f}\")"
]
},
{
"cell_type": "markdown",
"id": "9160d09c",
"metadata": {},
"source": [
"\n",
"\n",
"## Optimizar los orbitales\n",
"\n",
"La optimización orbital busca una rotación orbital que reduzca la energía variacional\n",
"\n",
"$$\n",
"E = \\langle \\psi | \\mathcal{U}^\\dagger\\, H\\, \\mathcal{U} | \\psi \\rangle\n",
"$$\n",
"\n",
"de la aproximación del estado fundamental SQD $|\\psi\\rangle$. Una rotación orbital viene dada por\n",
"una matriz unitaria de tipo « $N \\times N$ » $\\mathbf{U}$ (donde $N$ es el número de orbitales espaciales), que\n",
"actúa sobre el estado de muchos cuerpos a través del operador\n",
"\n",
"$$\n",
"\\mathcal{U} = \\exp\\left[\\sum_{pq, \\sigma} \\log(\\mathbf{U})_{pq}\\, a^\\dagger_{p\\sigma} a_{q\\sigma}\\right].\n",
"$$\n",
"\n",
"`ffsim.optimize_orbitals` devuelve la **matriz** $\\mathbf{U}$, y aplicarla a la\n",
"base orbital (mediante `hamiltonian.rotated`) equivale a aplicar $\\mathcal{U}$ al\n",
"estado. Consulta la\n",
"[explicación sobre la rotación orbital de ffsim](https://qiskit-community.github.io/ffsim/explanations/orbital-rotation.html)\n",
"para obtener más detalles.\n",
"\n",
"Dado que la rotación de los orbitales modifica el hamiltoniano que percibe el subespacio, alternamos\n",
"dos pasos hasta que la energía deja de mejorar:\n",
"\n",
"1. **Diagonalizar** el hamiltoniano en la base actual sobre el conjunto fijo de\n",
" configuraciones.\n",
"2. **Optimiza los orbitales** determinando la rotación que minimice la energía del\n",
" estado resultante y, a continuación, gira las integrales a la nueva base.\n",
"\n",
"Delegamos el paso de rotación orbital a\n",
"[`ffsim.optimize_orbitals`](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.optimize_orbitals),\n",
"que calcula la rotación que minimiza la energía a partir de las matrices de densidad reducidas\n",
"(RDM) de un cuerpo y de dos cuerpos del estado. Véase\n",
"[el art. II A 4](https://arxiv.org/pdf/2405.05068) para más detalles.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "006745b5",
"metadata": {},
"source": [
"\n",
"\n",
"### ¿Por qué la optimización orbital resulta útil en este caso?\n",
"\n",
"La base de orbitales moleculares (OM) SCF es estacionaria con respecto a las rotaciones orbitales para el problema -CI *completo* . Sin embargo, el método SQD opera en un pequeño subespacio truncado (en este caso, unos\n",
"cientos de cadenas de CI de entre unos 19 millones de determinantes de CI completos), para el cual la base MO\n",
"no suele ser óptima, por lo que la rotación de los orbitales reduce la energía que el subespacio\n",
"puede representar.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "46eb2f5a",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-16T01:37:42.597573Z",
"iopub.status.busy": "2026-07-16T01:37:42.597399Z",
"iopub.status.idle": "2026-07-16T01:37:42.654809Z",
"shell.execute_reply": "2026-07-16T01:37:42.654267Z"
}
},
"outputs": [],
"source": [
"import ffsim\n",
"from pyscf import fci\n",
"\n",
"# ffsim's ``MolecularHamiltonian`` uses the same \"chemist\" ordering for the two-body\n",
"# tensor as PySCF's ``eri``, and stores the nuclear repulsion energy as the constant\n",
"# term so that expectation values come out as total energies."
]
},
{
"cell_type": "markdown",
"id": "93179dc4",
"metadata": {},
"source": [
"\n",
"\n",
"### Diagonalización alternativa y optimización orbital\n",
"\n",
"Mantenemos el subespacio de diagonalización **fijado** a las configuraciones descubiertas por SQD\n",
"anteriormente, de modo que cada iteración aísle el efecto de la rotación de los orbitales. En cada\n",
"iteración:\n",
"\n",
"1. **Diagonaliza** el hamiltoniano sobre el subespacio fijo en la base actual, utilizando\n",
" el solucionador de CI seleccionado de PySCF's.\n",
"2. **Crea los RDM**\n",
" del estado resultante, que es todo lo que hace`ffsim.optimize_orbitals`\n",
" falta.\n",
"3. **Optimiza los orbitales** : `ffsim.optimize_orbitals` devuelve la rotación que minimiza la energía,\n",
" la cual aplicamos a las integrales para pasar a la base mejorada.\n",
"\n",
"Registramos la energía *antes de* cada paso de optimización. Dado que la base mejora en cada\n",
"iteración, esta secuencia disminuye de forma monótona hacia la mejor energía alcanzable en\n",
"el subespacio fijo.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "a0783e88",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-16T01:37:42.657394Z",
"iopub.status.busy": "2026-07-16T01:37:42.657213Z",
"iopub.status.idle": "2026-07-16T01:38:48.374706Z",
"shell.execute_reply": "2026-07-16T01:38:48.373697Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iteration 0: energy = -108.98452447\n",
"Iteration 1: energy = -108.99981993\n",
"Iteration 2: energy = -109.00585329\n",
"Iteration 3: energy = -109.00816569\n",
"Iteration 4: energy = -109.00936616\n",
"Iteration 5: energy = -109.01014322\n",
"Iteration 6: energy = -109.01069439\n",
"Iteration 7: energy = -109.01109308\n",
"Iteration 8: energy = -109.01138928\n",
"Iteration 9: energy = -109.01161411\n"
]
}
],
"source": [
"# Fix the diagonalization subspace to the configurations found by SQD.\n",
"ci_strings = (result.sci_state.ci_strs_a, result.sci_state.ci_strs_b)\n",
"nelec = (num_elec_a, num_elec_b)\n",
"\n",
"# Start from the MO basis in which we ran SQD.\n",
"hamiltonian_opt = ffsim.MolecularHamiltonian(\n",
" hcore, eri, constant=nuclear_repulsion_energy\n",
")\n",
"\n",
"num_iters = 10\n",
"for i in range(num_iters):\n",
" # Diagonalize over the fixed subspace in the current basis.\n",
" myci = fci.selected_ci.SelectedCI()\n",
" myci = fci.addons.fix_spin_(myci, ss=spin_sq)\n",
" _, amplitudes = fci.selected_ci.kernel_fixed_space(\n",
" myci,\n",
" hamiltonian_opt.one_body_tensor,\n",
" hamiltonian_opt.two_body_tensor,\n",
" num_orbitals,\n",
" nelec,\n",
" ci_strs=ci_strings,\n",
" )\n",
"\n",
" # Build the RDMs and record the energy before re-optimizing the orbitals.\n",
" dm1, dm2 = myci.make_rdm12(amplitudes, num_orbitals, nelec)\n",
" rdm = ffsim.ReducedDensityMatrix(dm1, dm2)\n",
" energy = rdm.expectation(hamiltonian_opt).real\n",
" print(f\"Iteration {i}: energy = {energy:.8f}\")\n",
"\n",
" # Rotate the Hamiltonian into the energy-minimizing basis for the next iteration.\n",
" # optimize_orbitals returns the unitary matrix U minimizing\n",
" # rdm.rotated(U).expectation(hamiltonian), equivalently\n",
" # rdm.expectation(hamiltonian.rotated(U.conj().T)), so we rotate by U^dagger.\n",
" orbital_rotation = ffsim.optimize_orbitals(rdm, hamiltonian_opt)\n",
" hamiltonian_opt = hamiltonian_opt.rotated(orbital_rotation.T.conj())"
]
},
{
"cell_type": "markdown",
"id": "a61dbda8",
"metadata": {},
"source": [
"\n",
"\n",
"### Compara los resultados\n",
"\n",
"La optimización orbital mejora la estimación del subespacio fijo, reduciendo en gran medida la diferencia con respecto a\n",
"la energía exacta, sin dejar de situarse por encima de ella.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "762b0903",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-16T01:38:48.377620Z",
"iopub.status.busy": "2026-07-16T01:38:48.377404Z",
"iopub.status.idle": "2026-07-16T01:38:53.111168Z",
"shell.execute_reply": "2026-07-16T01:38:53.110312Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact energy: -109.04667178\n",
"SQD energy (MO): -108.98469255\n",
"Energy after OO: -109.01178727\n"
]
}
],
"source": [
"# Diagonalize once more in the final optimized basis to report the improved energy.\n",
"myci = fci.selected_ci.SelectedCI()\n",
"myci = fci.addons.fix_spin_(myci, ss=spin_sq)\n",
"_, amplitudes = fci.selected_ci.kernel_fixed_space(\n",
" myci,\n",
" hamiltonian_opt.one_body_tensor,\n",
" hamiltonian_opt.two_body_tensor,\n",
" num_orbitals,\n",
" nelec,\n",
" ci_strs=ci_strings,\n",
")\n",
"dm1, dm2 = myci.make_rdm12(amplitudes, num_orbitals, nelec)\n",
"energy_after_oo = (\n",
" ffsim.ReducedDensityMatrix(dm1, dm2).expectation(hamiltonian_opt).real\n",
")\n",
"\n",
"print(f\"Exact energy: {exact_energy:.8f}\")\n",
"print(f\"SQD energy (MO): {sqd_energy:.8f}\")\n",
"print(f\"Energy after OO: {energy_after_oo:.8f}\")"
]
},
{
"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
}