{
"cells": [
{
"cell_type": "markdown",
"id": "frontmatter",
"metadata": {},
"source": [
"---\n",
"title: \"Verbesserung einer SQD-Schätzung durch Orbitaloptimierung\"\n",
"description: \"Verbesserung einer SQD-Schätzung durch Orbitaloptimierung für die neueste Version der samplebasierten Quantendiagonalisierung (SQD)\"\n",
"---\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "bb5a576d",
"metadata": {},
"source": [
"\n",
"\n",
"# Verbesserung einer SQD-Schätzung durch Orbitaloptimierung\n",
"\n",
"Bei der samplebasierten Quantendiagonalisierung (SQD) wird die Grundzustandsenergie durch\n",
"Diagonalisierung des Hamilton-Operators in einem festen Unterraum elektronischer Konfigurationen approximiert. Diese\n",
"Schätzung hängt von der Orbitalbasis ab, in der der Hamiltonoperator ausgedrückt wird, und\n",
"*die Orbitaloptimierung* (OO) nutzt diese Freiheit, um die Energie zu senken, ohne\n",
"den Unterraum zu vergrößern.\n",
"\n",
"In dieser Anleitung wird SQD auf einem „ $N_2$ “-Molekül angewendet und das Ergebnis anschließend durch eine Orbitaloptimierung verbessert,\n",
"wobei verwendet wird, um [`ffsim`](https://qiskit-community.github.io/ffsim/)\n",
"den Hamilton-Operator darzustellen und die energieminimierende Orbitalrotation zu ermitteln.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "77ab953b",
"metadata": {},
"source": [
"\n",
"\n",
"## SQD ausführen\n",
"\n",
"Wir erstellen die Molekülintegrale für $N_2$ in der Molekülorbital (MO)-Basis, erzeugen\n",
"gleichverteilte Zufallsstichproben und führen eine SQD-Berechnung durch, um eine Näherung für den Grundzustand zu erhalten.\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",
"## Die Orbitalen optimieren\n",
"\n",
"Bei der Orbitaloptimierung wird nach einer Orbitalrotation gesucht, die die Variationsenergie verringert\n",
"\n",
"$$\n",
"E = \\langle \\psi | \\mathcal{U}^\\dagger\\, H\\, \\mathcal{U} | \\psi \\rangle\n",
"$$\n",
"\n",
"der SQD-Grundzustandsnäherung $|\\psi\\rangle$. Eine Orbitalrotation wird durch\n",
"eine $N \\times N$ e unitäre Matrix $\\mathbf{U}$ ( $N$ ist die Anzahl der räumlichen Orbitale) festgelegt, die\n",
"über den Operator auf den Vielteilchenzustand einwirkt\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` gibt die **Matrix** $\\mathbf{U}$ zurück, und deren Anwendung auf die\n",
"Orbitalbasis (über `hamiltonian.rotated`) entspricht der Anwendung von $\\mathcal{U}$ auf den\n",
"Zustand. Weitere Informationen finden Sie in der\n",
"[Erklärung zur Orbitalrotation bei ffsim](https://qiskit-community.github.io/ffsim/explanations/orbital-rotation.html)\n",
".\n",
"\n",
"Da sich durch die Rotation der Orbitalen der für den Unterraum sichtbare Hamiltonoperator ändert, wechseln wir\n",
"zwischen zwei Schritten hin und her, bis sich die Energie nicht mehr verbessert:\n",
"\n",
"1. **Diagonalisieren Sie** den Hamiltonoperator in der aktuellen Basis über die feste Menge der\n",
" Konfigurationen.\n",
"2. **Optimieren Sie die Orbitale**, indem Sie die Rotation ermitteln, die die Energie des\n",
" resultierenden Zustands minimiert, und transformieren Sie anschließend die Integrale in die neue Basis.\n",
"\n",
"Den Schritt der Orbitalrotation delegieren wir an\n",
"[`ffsim.optimize_orbitals`](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.optimize_orbitals),\n",
"das die energieminimierende Rotation anhand der reduzierten Dichtematrizen (RDMs) für ein und zwei Teilchen\n",
"des Zustands ermittelt. Siehe\n",
"[Abschnitt II A 4](https://arxiv.org/pdf/2405.05068) für weitere Einzelheiten.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "006745b5",
"metadata": {},
"source": [
"\n",
"\n",
"### Warum die Orbitaloptimierung hier hilft\n",
"\n",
"Die SCF-Molekülorbitalbasis (MO) ist bezüglich Orbitalrotationen für das *vollständige* -CI-Problem stationär. Allerdings arbeitet SQD in einem kleinen, abgeschnittenen Unterraum (hier einige\n",
"hundert CI-Strings von insgesamt etwa 19 Millionen vollständigen CI-Determinanten), für den die MO-Basis\n",
"im Allgemeinen nicht optimal ist, sodass eine Rotation der Orbitale die Energie senkt, die der Unterraum\n",
"darstellen kann.\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",
"### Alternierende Diagonalisierung und Orbitaloptimierung\n",
"\n",
"Wir halten den Diagonalisierungsunterraum auf die oben durch SQD ermittelten Konfigurationen **festgelegt**,\n",
"sodass bei jeder Iteration der Effekt der Rotation der Orbitale isoliert wird. Jede\n",
"Iteration:\n",
"\n",
"1. **Diagonalisiert** den Hamilton-Operator über den festen Unterraum in der aktuellen Basis unter Verwendung\n",
" des PySCF's -Selected-CI-Lösers.\n",
"2. **Erstellt die RDMs**\n",
" des resultierenden Zustands – mehr ist nicht`ffsim.optimize_orbitals`\n",
" erforderlich.\n",
"3. **Optimiert die Orbitalen** : `ffsim.optimize_orbitals` Gibt die energieminimierende\n",
" Drehung zurück, die wir auf die Integrale anwenden, um in die verbesserte Basis überzugehen.\n",
"\n",
"*Vor* jedem Optimierungsschritt erfassen wir den Energiezustand. Da sich die Basis mit jeder\n",
"Iteration verbessert, nimmt diese Folge monoton in Richtung der besten erreichbaren Energie im\n",
"festgelegten Unterraum ab.\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",
"### Vergleiche die Ergebnisse\n",
"\n",
"Die Orbitaloptimierung verbessert die Schätzung im festen Unterraum und schließt einen Großteil der Lücke zur\n",
"exakten Energie, bleibt dabei jedoch über dieser.\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
}