{ "cells": [ { "cell_type": "markdown", "id": "ee388de5-2507-48ff-97a6-6777698d6256", "metadata": {}, "source": [ "---\n", "title: \"Diagonalizzazione quantistica di Krylov basata su campioni di un modello a reticolo fermionico\"\n", "description: \"Utilizza l'algoritmo di diagonalizzazione quantistica basato su campioni per simulare il modello di Anderson a singola impurità utilizzando hardware quantistico rumoroso.\"\n", "---\n", "\n", "{/* cspell:ignore fontdict fontsize nocc SQKD DMRG textrm varepsilon vecs pqrs ijkl */}\n", "\n", "\n", "\n", "# Diagonalizzazione quantistica di Krylov basata su campioni di un modello a reticolo fermionico\n", "\n", "*Stima di utilizzo: Nove secondi su un processore Heron r2 (NOTA: questa è solo una stima. Il tempo di esecuzione potrebbe variare)*\n", "\n" ] }, { "cell_type": "markdown", "id": "a1b2c3d4-e5f6-7890-abcd-ef1234567890", "metadata": {}, "source": [ "\n", "\n", "## Risultati di apprendimento\n", "\n", "Dopo aver seguito questo tutorial, gli utenti dovrebbero aver compreso:\n", "\n", "* Come utilizzare [l](https://github.com/Qiskit/qiskit-addon-sqd) 'add-on SQD Qiskit per approssimare l'energia dello stato fondamentale di un modello reticolare utilizzando stringhe di bit campionate da un'unità di elaborazione quantistica (QPU).\n", "* Come utilizzare [ffsim](https://github.com/qiskit-community/ffsim) per costruire circuiti di evoluzione temporale per la simulazione fermionica.\n", "* Come combinare campioni provenienti da più circuiti per la post-elaborazione con l'algoritmo di diagonalizzazione di Krylov basato sui campioni (SKQD).\n", "\n", "\n", "\n", "## Prerequisiti\n", "\n", "Consigliamo agli utenti di acquisire familiarità con i seguenti argomenti prima di seguire questo tutorial:\n", "\n", "* [Diagonalizzazione quantistica di un Hamiltoniano chimico basata su campioni](/docs/tutorials/sample-based-quantum-diagonalization)\n", "* [Diagonalizzazione quantistica di Krylov degli hamiltoniani su reticolo](/docs/tutorials/krylov-quantum-diagonalization)\n", "* [Qiskit primitives](/docs/guides/primitives)\n", "\n" ] }, { "cell_type": "markdown", "id": "dc5cc74e-06bf-45ac-a69b-81778138e08f", "metadata": {}, "source": [ "\n", "\n", "## Sfondo\n", "\n", "Questo tutorial mostra come utilizzare la diagonalizzazione quantistica basata su campioni (SQD) per stimare l'energia di stato fondamentale di un modello reticolare fermionico. In particolare, studiamo il modello di Anderson monodimensionale a singola impurità (SIAM), utilizzato per descrivere le impurità magnetiche incorporate nei metalli.\n", "\n", "Questa esercitazione segue un flusso di lavoro simile a quello dell'esercitazione correlata [Diagonalizzazione quantistica a campione di un'hamiltoniana chimica](/docs/tutorials/sample-based-quantum-diagonalization). Tuttavia, una differenza fondamentale sta nel modo in cui vengono costruiti i circuiti quantistici. L'altro tutorial utilizza un ansatz variazionale euristico, interessante per gli hamiltoniani della chimica con potenzialmente milioni di termini di interazione. D'altra parte, questo tutorial utilizza circuiti che approssimano l'evoluzione del tempo tramite l'hamiltoniana. Tali circuiti possono essere profondi, il che rende questo approccio migliore per le applicazioni ai modelli reticolari. I vettori di stato preparati da questi circuiti formano la base di un [sottospazio di Krylov](https://en.wikipedia.org/wiki/Krylov_subspace) e, di conseguenza, l'algoritmo converge in modo dimostrabile ed efficiente allo stato fondamentale, sotto opportune ipotesi.\n", "\n", "L'approccio utilizzato in questa esercitazione può essere visto come una combinazione delle tecniche utilizzate in SQD e nella [diagonalizzazione quantistica di Krylov (KQD)](https://arxiv.org/abs/2407.14431). L'approccio combinato viene talvolta definito diagonalizzazione quantistica di Krylov basata su campioni (SQKD). Per un tutorial sul metodo KQD, vedere [la diagonalizzazione quantistica di Krylov degli hamiltoniani reticolari](/docs/tutorials/krylov-quantum-diagonalization).\n", "\n", "Questa esercitazione si basa sul lavoro [\"Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization\"](https://arxiv.org/abs/2501.09702), a cui si rimanda per maggiori dettagli.\n", "\n", "\n", "\n", "### Modello di Anderson a singola impurità (SIAM)\n", "\n", "L'hamiltoniana SIAM monodimensionale è una somma di tre termini:\n", "\n", "$$\n", "H = H_{\\textrm{imp}}+ H_\\textrm{bath} + H_\\textrm{hyb},\n", "$$\n", "\n", "Dove\n", "\n", "$$\n", "\\begin{align*}\n", " H_\\textrm{imp} &= \\varepsilon \\left( \\hat{n}_{d\\uparrow} + \\hat{n}_{d\\downarrow} \\right) + U \\hat{n}_{d\\uparrow}\\hat{n}_{d\\downarrow}, \\\\\n", " H_\\textrm{bath} &= -t \\sum_{\\substack{\\mathbf{j} = 0\\\\ \\sigma\\in \\{\\uparrow, \\downarrow\\}}}^{L-1} \\left(\\hat{c}^\\dagger_{\\mathbf{j}, \\sigma}\\hat{c}_{\\mathbf{j}+1, \\sigma} + \\hat{c}^\\dagger_{\\mathbf{j}+1, \\sigma}\\hat{c}_{\\mathbf{j}, \\sigma} \\right), \\\\\n", " H_\\textrm{hyb} &= V\\sum_{\\sigma \\in \\{\\uparrow, \\downarrow \\}} \\left(\\hat{d}^\\dagger_\\sigma \\hat{c}_{0, \\sigma} + \\hat{c}^\\dagger_{0, \\sigma} \\hat{d}_{\\sigma} \\right).\n", "\\end{align*}\n", "$$\n", "\n", "Qui, $c^\\dagger_{\\mathbf{j},\\sigma}/c_{\\mathbf{j},\\sigma}$ sono gli operatori fermionici di creazione/annientamento per il sito del bagno $\\mathbf{j}^{\\textrm{th}}$ con spin $\\sigma$, $\\hat{d}^\\dagger_{\\sigma}/\\hat{d}_{\\sigma}$ sono gli operatori di creazione/annientamento per il modo dell'impurità, e $\\hat{n}_{d\\sigma} = \\hat{d}^\\dagger_{\\sigma} \\hat{d}_{\\sigma}$. $t$, $U$, e $V$ sono numeri reali che descrivono le interazioni di hopping, on-site, ibridazione, e $\\varepsilon$ è un numero reale che specifica il potenziale chimico.\n", "\n", "Si noti che l'hamiltoniana è un'istanza specifica dell'hamiltoniana generica interazione-elettrone,\n", "\n", "$$\n", "\\begin{align*}\n", " H &= \\sum_{\\substack{p, q \\\\ \\sigma}} h_{pq} \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{q\\sigma} + \\sum_{\\substack{p, q, r, s \\\\ \\sigma \\tau}} \\frac{h_{pqrs}}{2} \\hat{a}^\\dagger_{p\\sigma} \\hat{a}^\\dagger_{q\\tau} \\hat{a}_{s\\tau} \\hat{a}_{r\\sigma} \\\\\n", " &= H_1 + H_2,\n", "\\end{align*}\n", "$$\n", "\n", "dove $H_1$ consiste di termini a un corpo, che sono quadratici negli operatori di creazione e annichilazione fermionici, e $H_2$ consiste di termini a due corpi, che sono quartici. Per il SIAM,\n", "\n", "$$\n", "H_2 = U \\hat{n}_{d\\uparrow}\\hat{n}_{d\\downarrow}\n", "$$\n", "\n", "e $H_1$ contiene il resto dei termini dell'hamiltoniana. Per rappresentare programmaticamente l'hamiltoniana, memorizziamo la matrice $h_{pq}$ e il tensore $h_{pqrs}$.\n", "\n", "\n", "\n", "### Basi di posizione e quantità di moto\n", "\n", "A causa della simmetria traslazionale approssimativa in $H_\\textrm{bath}$, non ci aspettiamo che lo stato fondamentale sia rado nella base di posizione (la base orbitale in cui l'Hamiltoniana è specificata sopra). Le prestazioni di SQD sono garantite solo se lo stato fondamentale è rado, cioè ha un peso significativo solo su un piccolo numero di stati base computazionali. Per migliorare la sparsità dello stato fondamentale, eseguiamo la simulazione nella base orbitale in cui $H_\\textrm{bath}$ è diagonale. Chiamiamo questa base la *base del momento*. Poiché $H_\\textrm{bath}$ è un'hamiltoniana fermionica quadratica, può essere efficientemente diagonalizzata da una rotazione orbitale.\n", "\n", "\n", "\n", "### Evoluzione temporale approssimativa mediante l'Hamiltoniano\n", "\n", "Per approssimare l'evoluzione temporale dell'hamiltoniana, utilizziamo una decomposizione di Trotter-Suzuki del secondo ordine,\n", "\n", "$$\n", " e^{-i \\Delta t H} \\approx e^{-i\\frac{\\Delta t}{2} H_2} e^{-i\\Delta t H_1} e^{-i\\frac{\\Delta t}{2} H_2}.\n", "$$\n", "\n", "Sotto la [trasformazione di Jordan-Wigner](https://en.wikipedia.org/wiki/Jordan%E2%80%93Wigner_transformation), l'evoluzione temporale di $H_2$ equivale a un singolo gate [CPhase](/docs/api/qiskit/qiskit.circuit.library.CPhaseGate) tra gli orbitali di spin-up e spin-down nel sito dell'impurità. Poiché $H_1$ è un'hamiltoniana fermionica quadratica, l'evoluzione temporale di $H_1$ equivale a una rotazione orbitale.\n", "\n", "Gli stati base di Krylov $\\{ |\\psi_k\\rangle \\}_{k=0}^{D-1}$, dove $D$ è la dimensione del sottospazio di Krylov, sono formati dall'applicazione ripetuta di un singolo passo di Trotter, quindi\n", "\n", "$$\n", " |\\psi_k\\rangle \\approx \\left[e^{-i\\frac{\\Delta t}{2} H_2} e^{-i\\Delta t H_1} e^{-i\\frac{\\Delta t}{2} H_2} \\right]^k\\ket{\\psi_0}.\n", "$$\n", "\n", "Nel seguente flusso di lavoro basato su SQD, campioneremo da questo insieme di circuiti e post-processeremo l'insieme combinato di bitstring con SQD. Questo approccio contrasta con quello utilizzato nel tutorial correlato [Diagonalizzazione quantistica basata su campioni di un'hamiltoniana chimica](/docs/tutorials/sample-based-quantum-diagonalization), dove i campioni sono stati estratti da un singolo circuito variazionale euristico.\n", "\n" ] }, { "cell_type": "markdown", "id": "df66d697-102c-4a2c-80f3-f67fdda05573", "metadata": {}, "source": [ "\n", "\n", "## Requisiti\n", "\n", "Prima di iniziare questa esercitazione, assicuratevi di aver installato quanto segue:\n", "\n", "* Qiskit SDK v1.0 o versioni successive, con supporto [alla visualizzazione](/docs/api/qiskit/visualization)\n", "* Qiskit Runtime v0.22 o successivamente (`pip install qiskit-ibm-runtime`)\n", "* Componente aggiuntivo SQD Qiskit v0.11 o versioni successive (`pip install qiskit-addon-sqd`)\n", "* ffsim v0.0.72 o versioni successive (`pip install ffsim`)\n", "\n" ] }, { "cell_type": "markdown", "id": "c3d4e5f6-a7b8-9012-cdef-123456789012", "metadata": {}, "source": [ "\n", "\n", "## Esempio di simulatore su piccola scala\n", "\n" ] }, { "cell_type": "markdown", "id": "8540487a-8033-49c2-9f30-022336105f64", "metadata": {}, "source": [ "\n", "\n", "### Passaggio 1: mappare il problema su un circuito quantistico\n", "\n" ] }, { "cell_type": "markdown", "id": "4e6652bd-c97d-4f2f-98a7-d44857805cf1", "metadata": {}, "source": [ "In primo luogo, generiamo l'hamiltoniana SIAM nella base di posizione. L'hamiltoniana è rappresentata dalla matrice $h_{pq}$ e dal tensore $h_{pqrs}$. Poi, la ruotiamo nella base del momento. Nella base di posizione, collochiamo l'impurità nel primo sito. Tuttavia, quando si passa alla base di quantità di moto, si sposta l'impurità in un sito centrale per facilitare le interazioni con gli altri orbitali.\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "b5cb9c28-4721-4141-8665-96885038e210", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pyscf.fci\n", "\n", "\n", "def siam_hamiltonian(\n", " norb: int,\n", " hopping: float,\n", " onsite: float,\n", " hybridization: float,\n", " chemical_potential: float,\n", ") -> tuple[np.ndarray, np.ndarray]:\n", " \"\"\"Hamiltonian for the single-impurity Anderson model.\"\"\"\n", " # Place the impurity on the first site\n", " impurity_orb = 0\n", "\n", " # One body matrix elements in the \"position\" basis\n", " h1e = np.zeros((norb, norb))\n", " np.fill_diagonal(h1e[:, 1:], -hopping)\n", " np.fill_diagonal(h1e[1:, :], -hopping)\n", " h1e[impurity_orb, impurity_orb + 1] = -hybridization\n", " h1e[impurity_orb + 1, impurity_orb] = -hybridization\n", " h1e[impurity_orb, impurity_orb] = chemical_potential\n", "\n", " # Two body matrix elements in the \"position\" basis\n", " h2e = np.zeros((norb, norb, norb, norb))\n", " h2e[impurity_orb, impurity_orb, impurity_orb, impurity_orb] = onsite\n", "\n", " return h1e, h2e\n", "\n", "\n", "def momentum_basis(norb: int) -> np.ndarray:\n", " \"\"\"Get the orbital rotation to change from the position to the momentum basis.\"\"\"\n", " n_bath = norb - 1\n", "\n", " # Orbital rotation that diagonalizes the bath (non-interacting system)\n", " hopping_matrix = np.zeros((n_bath, n_bath))\n", " np.fill_diagonal(hopping_matrix[:, 1:], -1)\n", " np.fill_diagonal(hopping_matrix[1:, :], -1)\n", " _, vecs = np.linalg.eigh(hopping_matrix)\n", "\n", " # Expand to include impurity\n", " orbital_rotation = np.zeros((norb, norb))\n", " # Impurity is on the first site\n", " orbital_rotation[0, 0] = 1\n", " orbital_rotation[1:, 1:] = vecs\n", "\n", " # Move the impurity to the center\n", " new_index = n_bath // 2\n", " perm = np.r_[1 : (new_index + 1), 0, (new_index + 1) : norb]\n", " orbital_rotation = orbital_rotation[:, perm]\n", "\n", " return orbital_rotation\n", "\n", "\n", "def rotated(\n", " h1e: np.ndarray, h2e: np.ndarray, orbital_rotation: np.ndarray\n", ") -> tuple[np.ndarray, np.ndarray]:\n", " \"\"\"Rotate the orbital basis of a Hamiltonian.\"\"\"\n", " h1e_rotated = np.einsum(\n", " \"ab,Aa,Bb->AB\",\n", " h1e,\n", " orbital_rotation,\n", " orbital_rotation.conj(),\n", " optimize=\"greedy\",\n", " )\n", " h2e_rotated = np.einsum(\n", " \"abcd,Aa,Bb,Cc,Dd->ABCD\",\n", " h2e,\n", " orbital_rotation,\n", " orbital_rotation.conj(),\n", " orbital_rotation,\n", " orbital_rotation.conj(),\n", " optimize=\"greedy\",\n", " )\n", " return h1e_rotated, h2e_rotated\n", "\n", "\n", "# Total number of spatial orbitals, including the bath sites and the impurity\n", "# This should be an even number\n", "norb = 8\n", "\n", "# System is half-filled\n", "nelec = (norb // 2, norb // 2)\n", "# One orbital is the impurity, the rest are bath sites\n", "n_bath = norb - 1\n", "\n", "# Hamiltonian parameters\n", "hybridization = 1.0\n", "hopping = 1.0\n", "onsite = 10.0\n", "chemical_potential = -0.5 * onsite\n", "\n", "# Generate Hamiltonian in position basis\n", "h1e, h2e = siam_hamiltonian(\n", " norb=norb,\n", " hopping=hopping,\n", " onsite=onsite,\n", " hybridization=hybridization,\n", " chemical_potential=chemical_potential,\n", ")\n", "\n", "# Rotate to momentum basis\n", "orbital_rotation = momentum_basis(norb)\n", "h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())\n", "# In the momentum basis, the impurity is placed in the center\n", "impurity_index = n_bath // 2\n", "\n", "# Use PySCF to compute the exact ground state energy\n", "reference_energy, _ = pyscf.fci.direct_spin1.kernel(h1e, h2e, norb, nelec)" ] }, { "cell_type": "markdown", "id": "e888edf2-7865-41a8-be57-6ddb72dd0cc7", "metadata": {}, "source": [ "Quindi, generiamo i circuiti per produrre gli stati della base di Krylov.\n", "Per ogni specie di spin, lo stato iniziale $\\ket{\\psi_0}$ è dato dalla sovrapposizione di tutte le possibili eccitazioni dei tre elettroni più vicini al livello di Fermi nei 4 modi vuoti più vicini a partire dallo stato $|00\\cdots 0011 \\cdots 11\\rangle$, e realizzato mediante l'applicazione di sette [XXPlusYYGates](/docs/api/qiskit/qiskit.circuit.library.XXPlusYYGate).\n", "Gli stati evoluti nel tempo sono prodotti da applicazioni successive di un passo di Trotter del secondo ordine.\n", "\n", "Per una descrizione più dettagliata di questo modello e di come sono stati progettati i circuiti, si rimanda a [\"Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization\"](https://arxiv.org/abs/2501.09702).\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "0f729f86-1814-4d5a-ae65-3f9614e103b3", "metadata": {}, "outputs": [], "source": [ "from typing import Sequence\n", "\n", "import ffsim\n", "import scipy\n", "from qiskit import QuantumCircuit, QuantumRegister\n", "from qiskit.circuit import CircuitInstruction, Qubit\n", "from qiskit.circuit.library import CPhaseGate, XGate, XXPlusYYGate\n", "\n", "\n", "def prepare_initial_state(qubits: Sequence[Qubit], norb: int, nocc: int):\n", " \"\"\"Prepare initial state.\"\"\"\n", " assert norb >= 8\n", " x_gate = XGate()\n", " rot = XXPlusYYGate(0.5 * np.pi, -0.5 * np.pi)\n", " for i in range(nocc):\n", " yield CircuitInstruction(x_gate, [qubits[i]])\n", " yield CircuitInstruction(x_gate, [qubits[norb + i]])\n", " for i in range(3):\n", " for j in range(nocc - i - 1, nocc + i, 2):\n", " yield CircuitInstruction(rot, [qubits[j], qubits[j + 1]])\n", " yield CircuitInstruction(\n", " rot, [qubits[norb + j], qubits[norb + j + 1]]\n", " )\n", " yield CircuitInstruction(rot, [qubits[j + 1], qubits[j + 2]])\n", " yield CircuitInstruction(\n", " rot, [qubits[norb + j + 1], qubits[norb + j + 2]]\n", " )\n", "\n", "\n", "def trotter_step(\n", " qubits: Sequence[Qubit],\n", " time_step: float,\n", " one_body_evolution: np.ndarray,\n", " h2e: np.ndarray,\n", " impurity_index: int,\n", " norb: int,\n", "):\n", " \"\"\"A Trotter step.\"\"\"\n", " # Assume the two-body interaction is just the on-site interaction of the impurity\n", " onsite = h2e[\n", " impurity_index, impurity_index, impurity_index, impurity_index\n", " ]\n", " # Two-body evolution for half the time\n", " yield CircuitInstruction(\n", " CPhaseGate(-0.5 * time_step * onsite),\n", " [qubits[impurity_index], qubits[norb + impurity_index]],\n", " )\n", " # One-body evolution for the full time\n", " yield CircuitInstruction(\n", " ffsim.qiskit.OrbitalRotationJW(norb, one_body_evolution), qubits\n", " )\n", " # Two-body evolution for half the time\n", " yield CircuitInstruction(\n", " CPhaseGate(-0.5 * time_step * onsite),\n", " [qubits[impurity_index], qubits[norb + impurity_index]],\n", " )\n", "\n", "\n", "# Time step\n", "time_step = 0.2\n", "# Number of Krylov basis states\n", "krylov_dim = 8\n", "\n", "# Initialize circuit\n", "qubits = QuantumRegister(2 * norb, name=\"q\")\n", "circuit = QuantumCircuit(qubits)\n", "\n", "# Generate initial state\n", "for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):\n", " circuit.append(instruction)\n", "circuit.measure_all()\n", "\n", "# Create list of circuits, starting with the initial state circuit\n", "circuits = [circuit.copy()]\n", "\n", "# Add time evolution circuits to the list\n", "one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)\n", "for i in range(krylov_dim - 1):\n", " # Remove measurements\n", " circuit.remove_final_measurements()\n", " # Append another Trotter step\n", " for instruction in trotter_step(\n", " qubits,\n", " time_step,\n", " one_body_evolution,\n", " h2e_momentum,\n", " impurity_index,\n", " norb,\n", " ):\n", " circuit.append(instruction)\n", " # Measure qubits\n", " circuit.measure_all()\n", " # Add a copy of the circuit to the list\n", " circuits.append(circuit.copy())" ] }, { "cell_type": "code", "execution_count": 3, "id": "9f2cc4d4-ecac-457a-bcae-558319668e1f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "circuits[0].draw(\"mpl\", scale=0.4, fold=-1)" ] }, { "cell_type": "code", "execution_count": 4, "id": "827976ec-4815-4707-80b1-e13fb2fef309", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "circuits[-1].draw(\"mpl\", scale=0.4, fold=-1)" ] }, { "cell_type": "markdown", "id": "b8ca6be4-61d9-47be-8099-8712c7ecc774", "metadata": {}, "source": [ "\n", "\n", "### Fase 2: Ottimizzazione del problema per l'esecuzione quantistica\n", "\n" ] }, { "cell_type": "markdown", "id": "a3304e1b-9c7c-4212-8744-d1c62292eced", "metadata": {}, "source": [ "Successivamente, ottimizziamo il circuito per un hardware specifico. Per ora, creeremo un backend generico con un numero specificato di qubit e un insieme di porte in cui i circuiti di evoluzione temporale si scompongono naturalmente.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "2d2fdbff-1e22-45af-a2eb-c334e4328c59", "metadata": {}, "outputs": [], "source": [ "from qiskit.providers.fake_provider import GenericBackendV2\n", "\n", "backend = GenericBackendV2(\n", " 2 * norb, basis_gates=[\"cp\", \"xx_plus_yy\", \"p\", \"x\"]\n", ")" ] }, { "cell_type": "markdown", "id": "ce3e6a52-7b99-49b6-8294-b005efa59cfc", "metadata": {}, "source": [ "A questo punto, usiamo Qiskit per transpilare i circuiti nel backend di destinazione.\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "c8643533-9fec-40bf-a307-da8839b1e444", "metadata": {}, "outputs": [], "source": [ "from qiskit.transpiler import generate_preset_pass_manager\n", "\n", "pass_manager = generate_preset_pass_manager(\n", " optimization_level=3, backend=backend\n", ")\n", "isa_circuits = pass_manager.run(circuits)" ] }, { "cell_type": "markdown", "id": "6cfd3eea-e2d9-40a5-a449-1d3d790a5f2d", "metadata": {}, "source": [ "\n", "\n", "### Passaggio 3: Eseguire utilizzando Qiskit primitives\n", "\n" ] }, { "cell_type": "markdown", "id": "ad48e3e6-1013-45a1-942c-63766fed5819", "metadata": {}, "source": [ "Dopo aver ottimizzato i circuiti per l'esecuzione hardware, siamo pronti a eseguirli sull'hardware di destinazione e a raccogliere campioni per la stima dell'energia dello stato di massa. Dopo aver utilizzato la primitiva Sampler per campionare le stringhe di bit di ciascun circuito, combiniamo tutti i risultati in un unico dizionario di conteggi e tracciamo le 20 stringhe più comunemente campionate.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "80eee553-60d6-4258-88ab-d8d120418c36", "metadata": {}, "outputs": [], "source": [ "from qiskit.visualization import plot_histogram\n", "from qiskit.primitives import StatevectorSampler\n", "\n", "# Sample from the circuits\n", "sampler = StatevectorSampler()\n", "job = sampler.run(isa_circuits, shots=500)" ] }, { "cell_type": "code", "execution_count": 8, "id": "10af4663-7375-4b50-bae6-9f3d5106457b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from qiskit.primitives import BitArray\n", "\n", "# Combine the shots from the individual Trotter circuits\n", "bit_array = BitArray.concatenate_shots(\n", " [result.data.meas for result in job.result()]\n", ")\n", "\n", "plot_histogram(bit_array.get_counts(), number_to_keep=20)" ] }, { "cell_type": "markdown", "id": "2aa74455-d16b-4ac3-a354-54a79d5c5759", "metadata": {}, "source": [ "\n", "\n", "### Fase 4: Post-elaborazione e restituzione del risultato nel formato classico desiderato\n", "\n" ] }, { "cell_type": "markdown", "id": "d3f7713c-a4e6-407f-94de-0e5cfeb1134c", "metadata": {}, "source": [ "Ora eseguiamo l'algoritmo SQD utilizzando la funzione `diagonalize_fermionic_hamiltonian` . Per spiegazioni sugli argomenti di questa funzione, consultare la [documentazione API](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian).\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "7609d1e1-e8ef-48e1-a965-97927f403163", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iteration 1\n", "\tSubsample 0\n", "\t\tEnergy: -13.4222953188441\n", "\t\tSubspace dimension: 529\n", "\tSubsample 1\n", "\t\tEnergy: -13.42237556285828\n", "\t\tSubspace dimension: 784\n", "\tSubsample 2\n", "\t\tEnergy: -13.422045397387413\n", "\t\tSubspace dimension: 529\n", "Iteration 2\n", "\tSubsample 0\n", "\t\tEnergy: -13.422379583305478\n", "\t\tSubspace dimension: 900\n", "\tSubsample 1\n", "\t\tEnergy: -13.422376197704326\n", "\t\tSubspace dimension: 841\n", "\tSubsample 2\n", "\t\tEnergy: -13.422421162849295\n", "\t\tSubspace dimension: 1089\n", "Iteration 3\n", "\tSubsample 0\n", "\t\tEnergy: -13.422421164670345\n", "\t\tSubspace dimension: 1156\n", "\tSubsample 1\n", "\t\tEnergy: -13.422421492737689\n", "\t\tSubspace dimension: 1156\n", "\tSubsample 2\n", "\t\tEnergy: -13.422421205869572\n", "\t\tSubspace dimension: 1156\n", "Iteration 4\n", "\tSubsample 0\n", "\t\tEnergy: -13.422421494558726\n", "\t\tSubspace dimension: 1225\n", "\tSubsample 1\n", "\t\tEnergy: -13.422421492737689\n", "\t\tSubspace dimension: 1156\n", "\tSubsample 2\n", "\t\tEnergy: -13.422421492737689\n", "\t\tSubspace dimension: 1156\n" ] } ], "source": [ "from qiskit_addon_sqd.fermion import (\n", " SCIResult,\n", " diagonalize_fermionic_hamiltonian,\n", ")\n", "\n", "# List to capture intermediate results\n", "result_history = []\n", "\n", "\n", "def callback(results: list[SCIResult]):\n", " result_history.append(results)\n", " iteration = len(result_history)\n", " print(f\"Iteration {iteration}\")\n", " for i, result in enumerate(results):\n", " print(f\"\\tSubsample {i}\")\n", " print(f\"\\t\\tEnergy: {result.energy}\")\n", " print(\n", " f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n", " )\n", "\n", "\n", "rng = np.random.default_rng(24)\n", "result = diagonalize_fermionic_hamiltonian(\n", " h1e_momentum,\n", " h2e_momentum,\n", " bit_array,\n", " samples_per_batch=100,\n", " norb=norb,\n", " nelec=nelec,\n", " num_batches=3,\n", " max_iterations=5,\n", " symmetrize_spin=True,\n", " callback=callback,\n", " seed=rng,\n", ")" ] }, { "cell_type": "markdown", "id": "3dee9c61-fc42-48e9-8888-af8fc831cd5c", "metadata": {}, "source": [ "La seguente cella di codice visualizza i risultati. Il primo grafico mostra l'energia calcolata in funzione del numero di iterazioni di recupero della configurazione, mentre il secondo grafico mostra l'occupazione media di ciascun orbitale spaziale dopo l'iterazione finale. Trattandosi di un problema così semplice, già la prima iterazione ci avvicina molto all'energia esatta (si noti la scala dell'asse y).\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "b6879566-8bf5-4c28-bfb6-b2686692e3d3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Reference energy: -13.42249\n", "SQD energy: -13.42242\n", "Absolute error: 0.00007\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "min_es = [\n", " min(result, key=lambda res: res.energy).energy\n", " for result in result_history\n", "]\n", "min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])\n", "\n", "# Data for energies plot\n", "x1 = range(len(result_history))\n", "\n", "# Data for avg spatial orbital occupancy\n", "y2 = np.sum(result.orbital_occupancies, axis=0)\n", "x2 = range(len(y2))\n", "\n", "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n", "\n", "# Plot energies\n", "axs[0].plot(x1, min_es, label=\"energy\", marker=\"o\")\n", "axs[0].set_xticks(x1)\n", "axs[0].set_xticklabels(x1)\n", "axs[0].axhline(\n", " y=reference_energy,\n", " color=\"#BF5700\",\n", " linestyle=\"--\",\n", " label=\"reference energy\",\n", ")\n", "axs[0].set_title(\"Approximated Ground State Energy vs SQD Iterations\")\n", "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n", "axs[0].set_ylabel(\"Energy\", fontdict={\"fontsize\": 12})\n", "axs[0].legend()\n", "\n", "# Plot orbital occupancy\n", "axs[1].bar(x2, y2, width=0.8)\n", "axs[1].set_xticks(x2)\n", "axs[1].set_xticklabels(x2)\n", "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n", "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n", "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n", "\n", "print(f\"Reference energy: {reference_energy:.5f}\")\n", "print(f\"SQD energy: {min_e:.5f}\")\n", "print(f\"Absolute error: {abs(min_e - reference_energy):.5f}\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "0c9f4976-d770-426a-822e-e6756c2cfbe7", "metadata": {}, "source": [ "\n", "\n", "### Verificare l'energia\n", "\n", "Si garantisce che l'energia restituita da SQD costituisca un limite superiore dell'energia effettiva dello stato fondamentale. È possibile verificare il valore dell'energia poiché SQD restituisce anche i coefficienti del vettore di stato che approssima lo stato fondamentale. È possibile calcolare l'energia a partire dal vettore di stato utilizzando le matrici di densità ridotte a una e a due particelle, come illustrato nella seguente cella di codice.\n", "\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "e2b9de72-61cf-49d3-a1b5-f043e4b16956", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Recomputed energy: -13.42242\n" ] } ], "source": [ "rdm1 = result.sci_state.rdm(rank=1, spin_summed=True)\n", "rdm2 = result.sci_state.rdm(rank=2, spin_summed=True)\n", "\n", "energy = np.sum(h1e_momentum * rdm1) + 0.5 * np.sum(h2e_momentum * rdm2)\n", "\n", "print(f\"Recomputed energy: {energy:.5f}\")" ] }, { "cell_type": "markdown", "id": "5221928a-79ff-4f54-90b4-1fe4d7739aae", "metadata": {}, "source": [ "\n", "\n", "## Esempio di hardware su larga scala\n", "\n", "Ora eseguiamo un esempio più complesso su una QPU reale.\n", "Per l'energia di riferimento, utilizziamo i risultati di un calcolo [DMRG](https://en.wikipedia.org/wiki/Density_matrix_renormalization_group) effettuato separatamente.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "933037d8-847e-4986-80da-5ac8d677b2ff", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using backend ibm_boston\n", "Iteration 1\n", "\tSubsample 0\n", "\t\tEnergy: -28.63965951544449\n", "\t\tSubspace dimension: 9801\n", "\tSubsample 1\n", "\t\tEnergy: -28.625588929202006\n", "\t\tSubspace dimension: 9409\n", "\tSubsample 2\n", "\t\tEnergy: -28.647371834135498\n", "\t\tSubspace dimension: 8281\n", "Iteration 2\n", "\tSubsample 0\n", "\t\tEnergy: -28.67213260849567\n", "\t\tSubspace dimension: 29584\n", "\tSubsample 1\n", "\t\tEnergy: -28.670340686158816\n", "\t\tSubspace dimension: 27225\n", "\tSubsample 2\n", "\t\tEnergy: -28.669976379525988\n", "\t\tSubspace dimension: 31329\n", "Iteration 3\n", "\tSubsample 0\n", "\t\tEnergy: -28.68622875601382\n", "\t\tSubspace dimension: 36100\n", "\tSubsample 1\n", "\t\tEnergy: -28.698569623143126\n", "\t\tSubspace dimension: 34225\n", "\tSubsample 2\n", "\t\tEnergy: -28.694848533971882\n", "\t\tSubspace dimension: 33856\n", "Iteration 4\n", "\tSubsample 0\n", "\t\tEnergy: -28.69883392844593\n", "\t\tSubspace dimension: 42025\n", "\tSubsample 1\n", "\t\tEnergy: -28.701289495200996\n", "\t\tSubspace dimension: 38025\n", "\tSubsample 2\n", "\t\tEnergy: -28.699319594978245\n", "\t\tSubspace dimension: 45369\n", "Iteration 5\n", "\tSubsample 0\n", "\t\tEnergy: -28.701936886834154\n", "\t\tSubspace dimension: 51076\n", "\tSubsample 1\n", "\t\tEnergy: -28.702468711812013\n", "\t\tSubspace dimension: 53824\n", "\tSubsample 2\n", "\t\tEnergy: -28.702298147575938\n", "\t\tSubspace dimension: 52900\n", "Reference energy: -28.70660\n", "SQD energy: -28.70247\n", "Absolute error: 0.00413\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from qiskit_ibm_runtime import SamplerV2 as Sampler\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "\n", "# Model parameters\n", "norb = 20\n", "nelec = (norb // 2, norb // 2)\n", "n_bath = norb - 1\n", "hybridization = 1.0\n", "hopping = 1.0\n", "onsite = 10.0\n", "chemical_potential = -0.5 * onsite\n", "\n", "# Generate Hamiltonian and orbital rotation\n", "h1e, h2e = siam_hamiltonian(\n", " norb=norb,\n", " hopping=hopping,\n", " onsite=onsite,\n", " hybridization=hybridization,\n", " chemical_potential=chemical_potential,\n", ")\n", "orbital_rotation = momentum_basis(norb)\n", "h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())\n", "impurity_index = n_bath // 2\n", "\n", "# Set reference energy to DMRG value computed separately\n", "reference_energy = -28.70659686\n", "\n", "# Algorithm parameters\n", "time_step = 0.2\n", "krylov_dim = 8\n", "\n", "# Construct circuits\n", "qubits = QuantumRegister(2 * norb, name=\"q\")\n", "circuit = QuantumCircuit(qubits)\n", "for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):\n", " circuit.append(instruction)\n", "circuit.measure_all()\n", "circuits = [circuit.copy()]\n", "one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)\n", "for i in range(krylov_dim - 1):\n", " circuit.remove_final_measurements()\n", " for instruction in trotter_step(\n", " qubits,\n", " time_step,\n", " one_body_evolution,\n", " h2e_momentum,\n", " impurity_index,\n", " norb,\n", " ):\n", " circuit.append(instruction)\n", " circuit.measure_all()\n", " circuits.append(circuit.copy())\n", "\n", "# Initialize hardware backend\n", "service = QiskitRuntimeService()\n", "backend = service.least_busy(\n", " operational=True, simulator=False, min_num_qubits=127\n", ")\n", "print(f\"Using backend {backend.name}\")\n", "\n", "# Transpile to backend\n", "pass_manager = generate_preset_pass_manager(\n", " optimization_level=3, backend=backend\n", ")\n", "isa_circuits = pass_manager.run(circuits)\n", "\n", "# Sample from the circuits\n", "sampler = Sampler(backend)\n", "sampler.options.environment.job_tags = [\"TUT_SKQD\"]\n", "job = sampler.run(isa_circuits, shots=500)\n", "\n", "# Combine the shots from the individual Trotter circuits\n", "bit_array = BitArray.concatenate_shots(\n", " [result.data.meas for result in job.result()]\n", ")\n", "\n", "# Run configuration recovery and diagonalization\n", "result_history = []\n", "\n", "\n", "def callback(results: list[SCIResult]):\n", " result_history.append(results)\n", " iteration = len(result_history)\n", " print(f\"Iteration {iteration}\")\n", " for i, result in enumerate(results):\n", " print(f\"\\tSubsample {i}\")\n", " print(f\"\\t\\tEnergy: {result.energy}\")\n", " print(\n", " f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n", " )\n", "\n", "\n", "rng = np.random.default_rng(24)\n", "result = diagonalize_fermionic_hamiltonian(\n", " h1e_momentum,\n", " h2e_momentum,\n", " bit_array,\n", " samples_per_batch=100,\n", " norb=norb,\n", " nelec=nelec,\n", " num_batches=3,\n", " max_iterations=5,\n", " symmetrize_spin=True,\n", " callback=callback,\n", " seed=rng,\n", ")\n", "\n", "\n", "# Plot results\n", "min_es = [\n", " min(result, key=lambda res: res.energy).energy\n", " for result in result_history\n", "]\n", "min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])\n", "x1 = range(len(result_history))\n", "y2 = np.sum(result.orbital_occupancies, axis=0)\n", "x2 = range(len(y2))\n", "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n", "axs[0].plot(x1, min_es, label=\"energy\", marker=\"o\")\n", "axs[0].set_xticks(x1)\n", "axs[0].set_xticklabels(x1)\n", "axs[0].axhline(\n", " y=reference_energy,\n", " color=\"#BF5700\",\n", " linestyle=\"--\",\n", " label=\"reference energy\",\n", ")\n", "axs[0].set_title(\"Approximated Ground State Energy vs SQD Iterations\")\n", "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n", "axs[0].set_ylabel(\"Energy\", fontdict={\"fontsize\": 12})\n", "axs[0].legend()\n", "axs[1].bar(x2, y2, width=0.8)\n", "axs[1].set_xticks(x2)\n", "axs[1].set_xticklabels(x2)\n", "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n", "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n", "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n", "print(f\"Reference energy: {reference_energy:.5f}\")\n", "print(f\"SQD energy: {min_e:.5f}\")\n", "print(f\"Absolute error: {abs(min_e - reference_energy):.5f}\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "482ebea3-84b8-471b-bddc-b282e23162ad", "metadata": {}, "source": [ "\n", "\n", "## Passi successivi\n", "\n", "\n", " Se questo lavoro ti è sembrato interessante, potrebbero interessarti anche i seguenti contenuti:\n", "\n", " * [Diagonalizzazione quantistica basata su campioni di un hamiltoniano chimico](/docs/tutorials/sample-based-quantum-diagonalization) : un tutorial correlato che utilizza un approccio variazionale euristico al posto dei circuiti di Trotter\n", " * [Diagonalizzazione quantistica di Krylov degli hamiltoniani su reticolo](/docs/tutorials/krylov-quantum-diagonalization) : un tutorial sul metodo KQD\n", " * [Documentazione dell'API dell'add-on SQD](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian) - Riferimento per la `diagonalize_fermionic_hamiltonian` funzione\n", " * [*Algoritmo quantistico-centrico per la diagonalizzazione di Krylov basata su campioni*](https://arxiv.org/abs/2501.09702) - l'articolo su cui si basa questo tutorial\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "id": "a1b8767d", "source": "© IBM Corp., 2017-2026" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" }, "hours": 1.5, "qpuSeconds": 9 }, "nbformat": 4, "nbformat_minor": 4 }