{ "cells": [ { "cell_type": "markdown", "id": "ee388de5-2507-48ff-97a6-6777698d6256", "metadata": {}, "source": [ "---\n", "title: \"Diagonalização quantum de Krylov baseada em amostras de um modelo de rede fermiônica\"\n", "description: \"Use o algoritmo de diagonalização quântica baseado em amostras para simular o modelo Anderson de impureza única usando hardware quântico ruidoso.\"\n", "---\n", "\n", "{/* cspell:ignore fontdict fontsize nocc SQKD DMRG textrm varepsilon vecs pqrs ijkl */}\n", "\n", "\n", "\n", "# Diagonalização quantum de Krylov baseada em amostras de um modelo de rede fermiônica\n", "\n", "*Estimativa de uso: Nove segundos em um processador Heron r2 (OBSERVAÇÃO: esta é apenas uma estimativa. Seu tempo de execução pode variar)*\n", "\n" ] }, { "cell_type": "markdown", "id": "a1b2c3d4-e5f6-7890-abcd-ef1234567890", "metadata": {}, "source": [ "\n", "\n", "## Resultados do aprendizado\n", "\n", "Após concluir este tutorial, os usuários deverão compreender:\n", "\n", "* Como usar o [complemento SQD Qiskit](https://github.com/Qiskit/qiskit-addon-sqd) para estimar a energia do estado fundamental de um modelo de rede cristalina utilizando sequências de bits obtidas de uma unidade de processamento quântico (QPU).\n", "* Como usar [o ffsim](https://github.com/qiskit-community/ffsim) para construir circuitos de evolução temporal para simulação de fermiónicos.\n", "* Como combinar amostras de vários circuitos para pós-processamento com o algoritmo de diagonalização de Krylov baseado em amostras (SKQD).\n", "\n", "\n", "\n", "## Pré-requisitos\n", "\n", "Recomendamos que os usuários estejam familiarizados com os seguintes tópicos antes de seguir com este tutorial:\n", "\n", "* [Diagonalização quantum baseada em amostras de um Hamiltoniano químico](/docs/tutorials/sample-based-quantum-diagonalization)\n", "* [Diagonalização quântica de Krylov de hamiltonianos de rede](/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", "## Segundo plano\n", "\n", "Este tutorial mostra como usar a diagonalização quântica baseada em amostra (SQD) para estimar a energia do estado fundamental de um modelo de rede fermiônica. Especificamente, estudamos o modelo unidimensional de impureza única de Anderson (SIAM), que é usado para descrever impurezas magnéticas incorporadas em metais.\n", "\n", "Este tutorial segue um fluxo de trabalho semelhante ao tutorial relacionado [Diagonalização quântica baseada em amostras de um Hamiltoniano químico](/docs/tutorials/sample-based-quantum-diagonalization). No entanto, uma diferença fundamental está na forma como os circuitos quânticos são construídos. O outro tutorial usa um ansatz variacional heurístico, que é atraente para Hamiltonianos químicos com potencialmente milhões de termos de interação. Por outro lado, este tutorial usa circuitos que aproximam a evolução do tempo pelo Hamiltoniano. Esses circuitos podem ser profundos, o que torna essa abordagem melhor para aplicações em modelos de rede. Os vetores de estado preparados por esses circuitos formam a base de um [subespaço de Krylov](https://en.wikipedia.org/wiki/Krylov_subspace) e, como resultado, o algoritmo converge de forma comprovada e eficiente para o estado fundamental, sob as suposições adequadas.\n", "\n", "A abordagem usada neste tutorial pode ser vista como uma combinação das técnicas usadas na SQD e na [diagonalização quântica de Krylov (KQD)](https://arxiv.org/abs/2407.14431). A abordagem combinada às vezes é chamada de diagonalização quântica de Krylov baseada em amostra (SQKD). Consulte [Diagonalização quântica de Krylov de Hamiltonianos de rede](/docs/tutorials/krylov-quantum-diagonalization) para obter um tutorial sobre o método KQD.\n", "\n", "Este tutorial é baseado no trabalho [\"Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization\"](https://arxiv.org/abs/2501.09702), que pode ser consultado para obter mais detalhes.\n", "\n", "\n", "\n", "### Modelo de Anderson de impureza única (SIAM)\n", "\n", "O Hamiltoniano SIAM unidimensional é uma soma de três termos:\n", "\n", "$$\n", "H = H_{\\textrm{imp}}+ H_\\textrm{bath} + H_\\textrm{hyb},\n", "$$\n", "\n", "em que\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", "Aqui, $c^\\dagger_{\\mathbf{j},\\sigma}/c_{\\mathbf{j},\\sigma}$ são os operadores fermiônicos de criação/aniquilação para o local de banho $\\mathbf{j}^{\\textrm{th}}$ com spin $\\sigma$, $\\hat{d}^\\dagger_{\\sigma}/\\hat{d}_{\\sigma}$ são operadores de criação/aniquilação para o modo de impureza e $\\hat{n}_{d\\sigma} = \\hat{d}^\\dagger_{\\sigma} \\hat{d}_{\\sigma}$. $t$, $U$ e $V$ são números reais que descrevem as interações de hopping, no local e de hibridização, e $\\varepsilon$ é um número real que especifica o potencial químico.\n", "\n", "Observe que o Hamiltoniano é uma instância específica do Hamiltoniano genérico de interação-elétron,\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", "em que $H_1$ consiste em termos de um corpo, que são quadráticos nos operadores de criação e aniquilação fermiônica, e $H_2$ consiste em termos de dois corpos, que são quárticos. Para a SIAM,\n", "\n", "$$\n", "H_2 = U \\hat{n}_{d\\uparrow}\\hat{n}_{d\\downarrow}\n", "$$\n", "\n", "e $H_1$ contém o restante dos termos no Hamiltoniano. Para representar o Hamiltoniano de forma programática, armazenamos a matriz $h_{pq}$ e o tensor $h_{pqrs}$.\n", "\n", "\n", "\n", "### Bases de posição e impulso\n", "\n", "Devido à simetria translacional aproximada em $H_\\textrm{bath}$, não esperamos que o estado fundamental seja esparso na base de posição (a base orbital na qual o Hamiltoniano é especificado acima). O desempenho do SQD é garantido somente se o estado básico for esparso, ou seja, ele tem peso significativo em apenas um pequeno número de estados da base computacional. Para melhorar a esparsidade do estado fundamental, realizamos a simulação na base orbital na qual $H_\\textrm{bath}$ é diagonal. Chamamos essa base de *base momentânea*. Como o $H_\\textrm{bath}$ é um Hamiltoniano fermiônico quadrático, ele pode ser diagonalizado com eficiência por uma rotação orbital.\n", "\n", "\n", "\n", "### Evolução temporal aproximada pelo hamiltoniano\n", "\n", "Para aproximar a evolução do tempo pelo Hamiltoniano, usamos uma decomposição de Trotter-Suzuki de segunda ordem,\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", "Sob a [transformação de Jordan-Wigner](https://en.wikipedia.org/wiki/Jordan%E2%80%93Wigner_transformation), a evolução do tempo por $H_2$ equivale a uma única porta [CPhase](/docs/api/qiskit/qiskit.circuit.library.CPhaseGate) entre os orbitais de spin-up e spin-down no local da impureza. Como o $H_1$ é um Hamiltoniano fermiônico quadrático, a evolução do tempo pelo $H_1$ equivale a uma rotação orbital.\n", "\n", "Os estados da base de Krylov $\\{ |\\psi_k\\rangle \\}_{k=0}^{D-1}$, em que $D$ é a dimensão do subespaço de Krylov, são formados pela aplicação repetida de uma única etapa de Trotter, de modo que\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", "No fluxo de trabalho baseado em SQD a seguir, faremos uma amostragem desse conjunto de circuitos e pós-processaremos o conjunto combinado de bitstrings com SQD. Essa abordagem contrasta com a usada no tutorial relacionado [Diagonalização quântica baseada em amostras de um Hamiltoniano químico](/docs/tutorials/sample-based-quantum-diagonalization), em que as amostras foram extraídas de um único circuito variacional heurístico.\n", "\n" ] }, { "cell_type": "markdown", "id": "df66d697-102c-4a2c-80f3-f67fdda05573", "metadata": {}, "source": [ "\n", "\n", "## Requisitos\n", "\n", "Antes de iniciar este tutorial, verifique se você tem os seguintes itens instalados:\n", "\n", "* Qiskit SDK v1.0 ou posterior, com suporte [para visualização](/docs/api/qiskit/visualization)\n", "* Qiskit Runtime v0.22 ou mais tarde (`pip install qiskit-ibm-runtime`)\n", "* Complemento SQD Qiskit v0.11 ou posterior (`pip install qiskit-addon-sqd`)\n", "* ffsim v0.0.72 ou posterior (`pip install ffsim`)\n", "\n" ] }, { "cell_type": "markdown", "id": "c3d4e5f6-a7b8-9012-cdef-123456789012", "metadata": {}, "source": [ "\n", "\n", "## Exemplo de simulador em pequena escala\n", "\n" ] }, { "cell_type": "markdown", "id": "8540487a-8033-49c2-9f30-022336105f64", "metadata": {}, "source": [ "\n", "\n", "### Passo 1: Mapear o problema para um circuito quântico\n", "\n" ] }, { "cell_type": "markdown", "id": "4e6652bd-c97d-4f2f-98a7-d44857805cf1", "metadata": {}, "source": [ "Primeiro, geramos o Hamiltoniano SIAM na base de posição. O Hamiltoniano é representado pela matriz $h_{pq}$ e o tensor $h_{pqrs}$. Em seguida, nós o rotacionamos para a base do momento. Na base de posição, colocamos a impureza no primeiro local. No entanto, quando giramos para a base de momento, movemos a impureza para um local central para facilitar as interações com outros orbitais.\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": [ "Em seguida, geramos os circuitos para produzir os estados da base de Krylov.\n", "Para cada espécie de spin, o estado inicial $\\ket{\\psi_0}$ é dado pela superposição de todas as possíveis excitações dos três elétrons mais próximos do nível de Fermi nos 4 modos vazios mais próximos, começando pelo estado $|00\\cdots 0011 \\cdots 11\\rangle$, e realizado pela aplicação de sete [XXPlusYYGates](/docs/api/qiskit/qiskit.circuit.library.XXPlusYYGate).\n", "Os estados evoluídos no tempo são produzidos por aplicações sucessivas de uma etapa de Trotter de segunda ordem.\n", "\n", "Para obter uma descrição mais detalhada desse modelo e de como os circuitos são projetados, consulte [\"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", "### Etapa 2: Otimizar o problema para execução quântica\n", "\n" ] }, { "cell_type": "markdown", "id": "a3304e1b-9c7c-4212-8744-d1c62292eced", "metadata": {}, "source": [ "Em seguida, otimizamos o circuito para um hardware específico. Por enquanto, vamos criar um backend genérico com um número específico de qubits e um conjunto de portas no qual os circuitos de evolução temporal se decompõem 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": [ "Agora, usamos o Qiskit para transpilar os circuitos para o backend de destino.\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", "### Passo 3: Execute usando Qiskit primitives\n", "\n" ] }, { "cell_type": "markdown", "id": "ad48e3e6-1013-45a1-942c-63766fed5819", "metadata": {}, "source": [ "Depois de otimizar os circuitos para execução em hardware, estamos prontos para executá-los no hardware de destino e coletar amostras para a estimativa de energia do estado fundamental. Depois de usar a primitiva Sampler para obter amostras de bitstrings de cada circuito, combinamos todos os resultados em um único dicionário de contagens e plotamos as 20 principais bitstrings mais comumente amostradas.\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", "### Etapa 4: Pós-processamento e retorno do resultado ao formato clássico desejado\n", "\n" ] }, { "cell_type": "markdown", "id": "d3f7713c-a4e6-407f-94de-0e5cfeb1134c", "metadata": {}, "source": [ "Agora, executamos o algoritmo SQD usando a função `diagonalize_fermionic_hamiltonian` . Consulte a [documentação da API](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian) para obter explicações sobre os argumentos dessa função.\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": [ "A célula de código a seguir representa graficamente os resultados. O primeiro gráfico mostra a energia calculada em função do número de iterações de recuperação da configuração, e o segundo gráfico mostra a ocupação média de cada orbital espacial após a iteração final. Como se trata de um problema tão pequeno, a primeira iteração já nos aproxima bastante da energia exata (observe a escala do eixo 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", "### Verifique a energia\n", "\n", "A energia retornada pelo SQD é garantidamente um limite superior da energia real do estado fundamental. O valor da energia pode ser verificado, pois o SQD também retorna os coeficientes do vetor de estado que aproxima o estado fundamental. É possível calcular a energia a partir do vetor de estado utilizando suas matrizes de densidade reduzida de uma e duas partículas, conforme demonstrado na célula de código a seguir.\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", "## Exemplo de hardware em grande escala\n", "\n", "Agora, vamos executar um exemplo maior em uma QPU real.\n", "Para a energia de referência, utilizamos os resultados de um cálculo [DMRG](https://en.wikipedia.org/wiki/Density_matrix_renormalization_group) realizado separadamente.\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", "## Próximas etapas\n", "\n", "\n", " Se você achou este trabalho interessante, talvez se interesse pelo seguinte material:\n", "\n", " * [Diagonalização quântica baseada em amostras de um hamiltoniano químico](/docs/tutorials/sample-based-quantum-diagonalization) — um tutorial relacionado que utiliza um enfoque variacional heurístico em vez de circuitos de Trotter\n", " * [Diagonalização quântica de Krylov de hamiltonianos de rede](/docs/tutorials/krylov-quantum-diagonalization) - um tutorial sobre o método KQD\n", " * [Documentação da API do complemento SQD](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian) - referência para a `diagonalize_fermionic_hamiltonian` função\n", " * [*Algoritmo centrado em quantos para a diagonalização de Krylov baseada em amostras*](https://arxiv.org/abs/2501.09702) — o artigo no qual este tutorial se baseia\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 }