{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "048b37e0",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"SqDRIFT algoritmo para estimativa do estado fundamental\"\n",
        "description: \"SqDRIFT combina dois algoritmos de destaque, o “ qDRIFT ” e o SKQD, para o problema da estimativa do estado fundamental, ao mesmo tempo em que reduz a profundidade do circuito.\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore cisolver ECORE combinatorially multiset */}\n",
        "\n",
        "<span id=\"sqdrift-algorithm-for-ground-state-estimation\" />\n",
        "\n",
        "# SqDRIFT algoritmo para estimativa do estado fundamental\n",
        "\n",
        "Estimativa de tempo de *execução: 180 segundos em um processador Heron r3 (OBSERVAÇÃO: trata-se apenas de uma estimativa. (O tempo de execução pode variar.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "sqdrift-cpp-note",
      "metadata": {},
      "source": [
        "<Admonition type=\"note\" title=\"Está procurando a versão em C++?\">\n",
        "  Este tutorial utiliza o site Python. Para a implementação em C++, incluindo o código-fonte e as instruções de compilação, consulte o [tutorial “ SqDRIFT ” em C++](https://github.com/Qiskit/documentation/tree/main/docs/tutorials/assets/sqdrift/cpp).\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5edbd059",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## Resultados do aprendizado\n",
        "\n",
        "* Aprenda a criar circuitos com menor profundidade em comparação com a técnica de Trotterização\n",
        "* Conheça um fluxo de trabalho completo para estimativa do estado fundamental utilizando o “ qDRIFT ” e o SQD\n",
        "* Aprenda a usar `qiskit-fermions` em conjunto com outros complementos do Qiskit para implementar esse fluxo de trabalho\n",
        "\n",
        "Este tutorial é apresentado como um caderno do Python para fins didáticos.\n",
        "\n",
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## Pré-requisitos\n",
        "\n",
        "* Leia a visão geral sobre [a diagonalização quântica baseada em amostras (SQD)](/docs/addons/qiskit-addon-sqd)\n",
        "* Leia a aula sobre [a Diagonalização Quântica de Krylov Baseada em Amostras (SKQD)](/learning/courses/quantum-diagonalization-algorithms/skqd)\n",
        "\n",
        "<span id=\"background\" />\n",
        "\n",
        "## Segundo plano\n",
        "\n",
        "[SqDRIFT](https://arxiv.org/abs/2508.02578) é uma variante do SKQD que substitui a necessidade de escolher um ansatz a partir do qual amostrar cadeias de bits por um conjunto de circuitos de evolução temporal construídos diretamente a partir do hamiltoniano alvo. Isso é alcançado por meio da subamostragem de operadores de evolução temporal menores a partir do hamiltoniano, com base em seus coeficientes, o que é conhecido como método de trotterização de qDRIFT.\n",
        "\n",
        "Este tutorial utiliza o [Qiskit Fermions](/docs/addons/qiskit-fermions) para criar circuitos fermiónicos mais naturais para o algoritmo “ [qDRIFT](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.070503) ”, seguido pela aplicação de etapas de layout e síntese fermiónicas antes de integrar os circuitos ao pipeline tradicional do Qiskit para execução em hardware.\n",
        "\n",
        "Seja o hamiltoniano da forma:\n",
        "\n",
        "$$\n",
        "H = \\sum_{i=1}^{N} c_i h_i\n",
        "$$\n",
        "\n",
        "onde, sem perda de generalidade, exigimos que $c_i > 0$ e que o maior valor próprio de $h_i$ seja igual, em valor absoluto, a $1$. Qualquer prefator com sinal ou complexo é absorvido por $h_i$, de modo que os coeficientes $c_i$ são pesos estritamente positivos, enquanto os $h_i$ determinam a direção de cada termo. Aqui, $N$ é o número de termos (ou, após o agrupamento, o número de grupos) no hamiltoniano; trata-se de uma propriedade do hamiltoniano e é diferente do número de operadores incluídos em um único circuito, denotado por $n$ a seguir.\n",
        "\n",
        "O algoritmo “ qDRIFT ” realiza, então, para o tempo-alvo $t$, algum operador $V_k$, em que $k$ vai de $1 \\cdots K$ e representa o circuito $k_{th}$ SqDRIFT, definido como:\n",
        "\n",
        "$$\n",
        "V_k = \\prod_{j=1}^{n} e^{-i h_{k_j} \\lambda t / n }\n",
        "$$\n",
        "\n",
        "Aqui, $n$ é o número de operadores amostrados por circuito e $K$ é o número de circuitos no conjunto. O produto é calculado sobre os sorteios $n$, e não sobre todos os termos hamiltonianos $N$; e, como os termos são sorteados com reposição, o mesmo $h_i$ pode aparecer mais de uma vez em um único $V_k$.\n",
        "\n",
        "A quantidade:\n",
        "\n",
        "$$\n",
        "\\lambda = \\sum_{i=1}^{N} c_i\n",
        "$$\n",
        "\n",
        "é a norma de $L_1$ dos coeficientes; assim, cada uma das etapas $n$ evolui durante o mesmo intervalo de tempo $\\lambda t / n$, independentemente do termo que tenha sido sorteado. A uniformidade do ângulo de passo é a característica marcante de qDRIFT: : um coeficiente influencia o resultado pela *frequência com que* seu termo é extraído, e não pelo grau de rotação desse termo. Os índices são extraídos da distribuição:\n",
        "\n",
        "$$\n",
        "P[k_i] = \\frac{c_i}{\\lambda}\n",
        "$$\n",
        "\n",
        "portanto, a série $(k_1, \\ldots, k_n)$ é uma sequência aleatória de índices de termos extraídos dessa distribuição. Como os $c_i$ são positivos e somam $\\lambda$, trata-se de uma distribuição de probabilidade normalizada, e a esperança do canal resultante ao longo dos sorteios aleatórios se aproxima da evolução sob $H$, com um erro que diminui à medida que $n$ cresce. Observe que o erro de aproximação depende de $\\lambda$ e não do número de termos $N$.\n",
        "\n",
        "(O artigo “ SqDRIFT ” refere-se ao número de termos como “ $\\mathcal{N}$ ” e ao comprimento da sequência como “ $N$ ”; usamos aqui “ $N$ ” e “ $n$ ” para manter as duas noções claramente distintas.)\n",
        "\n",
        "Este tutorial mostra como gerar um conjunto desses circuitos aleatórios. Depois de criarmos esses circuitos, da mesma forma que criamos um subespaço de Krylov para diferentes operadores, amostramos sequências de bits a partir de vários desses operadores com diferentes parâmetros de tempo. Isso garante uma maior sobreposição entre os vetores do estado fundamental e as sequências de bits amostradas.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "cda4ef2a",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## Requisitos\n",
        "\n",
        "Antes de iniciar este tutorial, certifique-se de ter instalado\n",
        "\n",
        "* Um ambiente virtual do Python (>= 3.10 )\n",
        "* pip>=25.1\n",
        "* qiskit \\~= 2.5\n",
        "* qiskit-fermions==0.1.0 (Observe que o nome está no plural)\n",
        "* numpy\n",
        "* pyscf\n",
        "* qiskit-aer\n",
        "* qiskit-ibm-runtime\n",
        "* qiskit-addon-sqd\n",
        "\n",
        "Você pode instalar todos os pacotes necessários com o seguinte comando:\n",
        "\n",
        "```\n",
        "pip install \"qiskit~=2.5\" \"qiskit-fermions==0.1.0\" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy\n",
        "```\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e24facd1",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## Instalação\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "1e5d4294",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Third-party scientific computing\n",
        "import numpy as np\n",
        "\n",
        "# PySCF\n",
        "from pyscf import tools, ao2mo, fci\n",
        "\n",
        "# Qiskit core\n",
        "from qiskit import transpile\n",
        "from qiskit.primitives import BitArray\n",
        "\n",
        "# Qiskit Aer\n",
        "from qiskit_aer import AerSimulator\n",
        "\n",
        "# IBM Quantum Compute Service\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler\n",
        "\n",
        "# Qiskit Fermions\n",
        "from qiskit_fermions.operators.library import FCIDump\n",
        "from qiskit_fermions.operators import FermionOperator\n",
        "from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms\n",
        "from qiskit_fermions.operators.terms.grouping import (\n",
        "    group_terms_by_electronic_structure,\n",
        ")\n",
        "from qiskit_fermions.operators.terms.ordering import canonical_order\n",
        "from qiskit_fermions.circuit import FermionicCircuit\n",
        "from qiskit_fermions.circuit.library import Evolution\n",
        "from qiskit_fermions.transpiler import FermionicPassManager\n",
        "from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager\n",
        "from qiskit_fermions.transpiler.passes import QDriftTrotterization\n",
        "from qiskit_fermions.circuit.library import InitializeModes\n",
        "\n",
        "# Qiskit addon SQD\n",
        "from qiskit_addon_sqd.fermion import (\n",
        "    diagonalize_fermionic_hamiltonian,\n",
        "    SCIResult,\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "65478ff1",
      "metadata": {},
      "source": [
        "<span id=\"simulator-example\" />\n",
        "\n",
        "## Exemplo de simulador\n",
        "\n",
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### Etapa 1: Mapeie entradas clássicas para um problema quântico\n",
        "\n",
        "**Como ler e preparar o FCIDump**\n",
        "\n",
        "Para este tutorial, vamos carregar o hamiltoniano da estrutura eletrônica do nitrogênio ( N2 ). Existem também outras maneiras de criar operadores fermiónicos.  Consulte a documentação em [`qiskit_fermions.operators.library`](https://qiskit.github.io/qiskit-fermions/stable/0.1/pydoc/qiskit_fermions.operators.library.html#module-qiskit_fermions.operators.library).\n",
        "\n",
        "**Sobre este FCIDump.** O arquivo descreve `N2_sto_3g` uma molécula de nitrogênio ( $N_2$ ) na base mínima STO-3G, com uma separação interatômica de 1.09 $\\AA$, o comprimento de ligação de equilíbrio experimental. Seu cabeçalho declara `NORB=10`, `NELEC=14`, e `MS2=0`: 10 orbitais espaciais (portanto, 20 orbitais de spin e 20 qubits na representação de Jordan-Wigner), 14 elétrons em um singuleto de spin, ou seja, sete elétrons em $\\alpha$ e sete em $\\beta$. Todos os orbitais recebem o índice de simetria 1, ou seja, não é utilizada nenhuma simetria de grupo pontual. Por se tratar de um dump de STO-3G em espaço completo, nenhum orbital está congelado e o espaço de correlação é pequeno o suficiente para que uma energia de referência FCI exata possa ser calculada classicamente para fins de comparação, conforme mostrado na próxima célula.\n",
        "\n",
        "É possível regenerar um arquivo equivalente com o comando `PySCF:`\n",
        "\n",
        "```python\n",
        "from pyscf import gto, scf, tools\n",
        "\n",
        "mol = gto.M(atom=\"N 0 0 0; N 0 0 1.09\", basis=\"sto-3g\", symmetry=False)\n",
        "mf = scf.RHF(mol).run()\n",
        "tools.fcidump.from_scf(mf, \"N2_sto_3g\")\n",
        "```\n",
        "\n",
        "Como as integrais dependem dos orbitais SCF convergentes, um arquivo regenerado pode diferir do arquivo original quanto à fase ou à ordem dos orbitais; as energias totais não são afetadas.\n",
        "\n",
        "**Obtenção do arquivo.** Encontre o FCIDump neste [repositório: GitHub](https://github.com/Qiskit/documentation/tree/main/docs/tutorials/assets/sqdrift/fcidump_files). Você pode executar a célula abaixo para importá-la para o local esperado pelo restante do tutorial.\n",
        "\n",
        "Primeiro, usamos a função fornecida `cisolver` pelo pyscf para obter a energia de referência. Essa é a verdadeira energia do estado fundamental da molécula com a qual estamos trabalhando. Para isso, vamos primeiro declarar e `norb` `nelec`, que correspondem ao número de orbitais e ao número de elétrons, respectivamente. Em seguida, declaramos e `h1e` `h2e`, que são, respectivamente, os integrais de um e de dois elétrons. Todos esses dados também serão utilizados posteriormente para o SQD.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "f15f2d83",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g\n"
          ]
        }
      ],
      "source": [
        "import os\n",
        "from urllib.request import urlopen\n",
        "\n",
        "# The FCIDump is stored with this tutorial in the Qiskit documentation repository.\n",
        "FCIDUMP_URL = \"https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g\"\n",
        "FCIDUMP_PATH = \"assets/sqdrift/fcidump_files/N2_sto_3g\"\n",
        "\n",
        "if not os.path.exists(FCIDUMP_PATH):\n",
        "    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)\n",
        "    with urlopen(FCIDUMP_URL) as response:\n",
        "        contents = response.read()\n",
        "    with open(FCIDUMP_PATH, \"wb\") as f:\n",
        "        f.write(contents)\n",
        "    print(f\"Downloaded FCIDump to {FCIDUMP_PATH}\")\n",
        "else:\n",
        "    print(f\"Using existing FCIDump at {FCIDUMP_PATH}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "41b178cf",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Parsing assets/sqdrift/fcidump_files/N2_sto_3g\n",
            "Reference FCI Energy  = -107.6481842917 Ha\n",
            "Nuclear Repulsion Energy = 23.7887003074 Ha\n"
          ]
        }
      ],
      "source": [
        "name = \"assets/sqdrift/fcidump_files/N2_sto_3g\"\n",
        "\n",
        "fcidump = tools.fcidump.read(name)\n",
        "\n",
        "# Extract metadata from the FCIDump header\n",
        "norb = fcidump[\"NORB\"]  # number of spatial orbitals\n",
        "nelec = fcidump[\"NELEC\"]  # total number of electrons\n",
        "e_nuc = fcidump[\"ECORE\"]  # nuclear repulsion / core energy\n",
        "ms2 = fcidump[\"MS2\"]  # 2S (spin)\n",
        "\n",
        "num_elec_a = (nelec + ms2) // 2  # alpha electrons\n",
        "num_elec_b = (nelec - ms2) // 2  # beta  electrons\n",
        "\n",
        "# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)\n",
        "h1e = fcidump[\"H1\"]  # shape (norb, norb)\n",
        "h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)\n",
        "    1, fcidump[\"H2\"], norb\n",
        ")\n",
        "\n",
        "cisolver = fci.direct_spin1.FCI()\n",
        "cisolver.max_cycle = 200\n",
        "cisolver.conv_tol = 1e-12\n",
        "\n",
        "e_fci, _ = cisolver.kernel(\n",
        "    h1e,\n",
        "    h2e,\n",
        "    norb,\n",
        "    (num_elec_a, num_elec_b),\n",
        "    ecore=e_nuc,  # adds nuclear repulsion to the final energy\n",
        ")\n",
        "\n",
        "reference_energy = e_fci\n",
        "\n",
        "print(f\"Reference FCI Energy  = {reference_energy:.10f} Ha\")\n",
        "\n",
        "nuclear_repulsion_energy = fcidump[\"ECORE\"]\n",
        "print(f\"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7818e18a",
      "metadata": {},
      "source": [
        "**Carregando o hamiltoniano**\n",
        "\n",
        "Com os dados necessários preparados, lemos o hamiltoniano do arquivo FCI em um formato compatível com `qiskit-fermions`\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "565b92fc",
      "metadata": {},
      "outputs": [],
      "source": [
        "fcidump = FCIDump.from_file(name)\n",
        "hamiltonian = FermionOperator.from_fcidump(fcidump)\n",
        "num_modes = 2 * fcidump.norb"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d88bb917",
      "metadata": {},
      "source": [
        "**Fluxos de trabalho de férmions com `qiskit-fermions`**\n",
        "\n",
        "Primeiramente, mapearemos o hamiltoniano para um modelo de circuito fermiónico utilizando `qiskit-fermions`, que oferece passagens de transpiler e portas específicas para circuitos fermiónicos. Esses elementos serão utilizados posteriormente, antes das etapas tradicionais de transpilagem do Qiskit para este fluxo de trabalho.\n",
        "\n",
        "**Agrupamento de termos**\n",
        "\n",
        "Para garantir a reprodutibilidade dos resultados, primeiro utilizamos `canonical_order` para ordenar os termos com base apenas em sua estrutura. A ordem dos operadores na lista `canon` é, portanto, fixa. Isso garante a reprodutibilidade dos operadores criados, pois a função `pass `QDriftTrotterization` ` que utilizaremos nas próximas amostras seleciona índices aleatórios para criar os operadores `qDRIFT`.\n",
        "\n",
        "Nesta etapa, aproveitamos as diversas simetrias presentes no hamiltoniano da estrutura eletrônica, agrupando termos relacionados com coeficientes idênticos. Embora isso altere a distribuição dos coeficientes dos operadores a partir da qual o protocolo de qDRIFT e realiza suas amostragens, isso não afeta suas garantias de convergência. Fundamentalmente, o agrupamento de termos relacionados por simetria resulta em uma compensação favorável dos termos de Pauli e em uma profundidade de circuito globalmente menor ao se realizar a evolução temporal de um estado sob a ação desses termos.\n",
        "\n",
        "`qiskit-fermions` fornece a função `group_terms_by_electronic_structure` que realiza esse agrupamento para nós.\n",
        "\n",
        "Observe que a fórmula `group_terms_by_electronic_structure` pressupõe termos [com ordem normal](https://qiskit.github.io/qiskit-fermions/stable/0.1/stubs/qiskit_fermions.operators.FermionOperator.html#qiskit_fermions.operators.FermionOperator.normal_ordered).\n",
        "\n",
        "**Filtragem de termos diagonais**\n",
        "\n",
        "Removemos os termos diagonais do hamiltoniano usado para gerar os circuitos, de modo que $n$ qDRIFT os intervalos de amostragem sejam dedicados aos termos que movimentam a população entre as configurações. É melhor filtrar esses termos do hamiltoniano neste momento, antes que o `Evolution` portão seja construído na próxima etapa.\n",
        "\n",
        "Os termos em questão são aqueles que estão na diagonal na base de ocupação-número, ou seja, os produtos dos operadores numéricos $a^\\dagger_i a_i$. Três tipos de termos se enquadram nessa descrição:\n",
        "\n",
        "* o **deslocamento de energia constante**, um produto de operadores de número zero, cuja evolução temporal contribui apenas com uma fase global;\n",
        "* os **operadores de número individual** $n_i$, cuja evolução temporal se reduz a rotações de um único qubit $Z$;\n",
        "* os **produtos de ordem superior**, como $n_i n_j$.\n",
        "\n",
        "Por si só, nenhuma dessas ações transfere população entre configurações de ocupação e número; elas atuam apenas nas fases das configurações já existentes. No entanto, elas não são inertes: essas fases relativas contribuem para a interferência gerada pelos termos de excitação mais adiante no circuito; portanto, filtrá-las altera a evolução que é efetivamente gerada e pode alterar a distribuição amostral. Trata-se de uma aproximação deliberada na etapa de geração do circuito, feita para concentrar a amostragem nos termos de excitação, e não de uma etapa que deixe a distribuição amostrada inalterada. Ao contrário do agrupamento de simetria acima, que mantém intactas as garantias de convergência do método “ qDRIFT ”, esse filtro altera o operador que está sendo evoluído. Portanto, os circuitos não se aproximam mais da evolução sob o hamiltoniano completo, e os limites de erro de “ qDRIFT ” se aplicam ao operador filtrado, e não ao original. Isso é aceitável neste contexto porque os circuitos são apenas uma heurística de amostragem usada para propor configurações: nenhum termo é perdido na própria estimativa de energia, uma vez que o filtro se aplica apenas ao hamiltoniano utilizado para construir os circuitos, enquanto a diagonalização clássica, realizada posteriormente, utiliza o hamiltoniano completo, incluindo os termos diagonais. A precisão do SQD depende desse passo clássico, que permanece variacional no subespaço amostrado, independentemente de como as configurações foram propostas.\n",
        "\n",
        "A função `filter_diagonal_terms()` remove esses termos de um operador diretamente. Ela os identifica a partir de sua estrutura ordenada normalmente — o multiconjunto de modos de criação correspondendo ao multiconjunto de modos de aniquilação —, portanto, só é válida para um operador que já esteja ordenado normalmente. Essa suposição não é verificada durante a execução.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "464a4470",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "5060\n"
          ]
        }
      ],
      "source": [
        "# Apply automatic grouping\n",
        "canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))\n",
        "exit_code = group_terms_by_electronic_structure(\n",
        "    canon, num_modes, two_body_physicist_order=False\n",
        ")\n",
        "filter_diagonal_terms(canon)\n",
        "\n",
        "print(len(canon.groups))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "fe33ada4",
      "metadata": {},
      "source": [
        "Agora que agrupamos os termos no hamiltoniano, vamos definir os seguintes parâmetros para gerar o conjunto de circuitos:\n",
        "\n",
        "* O número de circuitos a serem gerados: `num_circuits`\n",
        "* O comprimento de cada circuito em termos de grupos de excitação: `num_exc`\n",
        "* O motivo para os diferentes tempos de evolução: `times`\n",
        "\n",
        "**Criação de circuitos fermiónicos**\n",
        "\n",
        "Agora, vamos criar circuitos fermiónicos para cada um dos intervalos de tempo. Cada circuito será composto por um único portão de evolução, com o tempo de evolução que definimos anteriormente. O operador de evolução é o hamiltoniano. Posteriormente, executamos etapas de transpilagem nesses circuitos para criar circuitos d qDRIFT.\n",
        "\n",
        "**Preparação do Ansatz**\n",
        "\n",
        "Preparamos o estado de Hartree-Fock utilizando a classe `InitializeModes` . Para o nitrogênio, o processo consiste simplesmente em aplicar portas X aos primeiros qubits `num_elec_a` e, em seguida, aos `num_elec_b` qubits, sendo que, no caso do nitrogênio, ambos os valores são iguais a sete. Esse estado representa os sete elétrons $\\alpha$ e os sete elétrons $\\beta$ do nitrogênio.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "140abcc6",
      "metadata": {},
      "outputs": [],
      "source": [
        "# SqDRIFT parameters\n",
        "times = [1.0, 10.0]  # Total evolution times used for the subspace creation\n",
        "num_exc = 10  # Number of excitation groups per circuit\n",
        "num_circuits = 200  # Number of circuits to generate\n",
        "\n",
        "\n",
        "init_circuits = []\n",
        "\n",
        "hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))\n",
        "\n",
        "for time in times:\n",
        "    evo_gate = Evolution(num_modes, canon, time)\n",
        "    circ = FermionicCircuit(num_modes)\n",
        "    circ.append(hf_gate, circ.modes)\n",
        "    circ.append(evo_gate, circ.modes)\n",
        "    init_circuits.append(circ)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c18ed9a5",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-problem-for-quantum-hardware-execution\" />\n",
        "\n",
        "### Etapa 2: Otimizar o problema para execução em hardware quântico\n",
        "\n",
        "Agora que temos nossos circuitos, vamos primeiro usar as etapas disponíveis em `qiskit-fermions` para realizar otimizações no nível fermionico e, em seguida, fazer a transpilagem do nosso circuito para o backend de nossa escolha. Como se trata de um experimento em simulador, vamos fazer isso primeiro para o AerSimulator.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7de434b9",
      "metadata": {},
      "source": [
        "**Cálculo do peso para cada grupo**\n",
        "\n",
        "Nesta etapa, realizamos a amostragem “ qDRIFT ” dos termos de forma estocástica, com probabilidades proporcionais aos seus coeficientes no hamiltoniano. A etapa do transpiler `qDRIFT` faz isso por nós. Agora podemos criar circuitos mais simples que podem ser executados no hardware de forma mais eficiente, apesar da conectividade limitada dos qubits, mesmo quando o hamiltoniano contém acoplamentos de longo alcance e termos superiores ao quadrático.\n",
        "Após o agrupamento de termos, ele seleciona os operadores com base em seus pesos. Para cada operador $h_i$, o peso $W_{h_i}$ é definido da seguinte forma:\n",
        "\n",
        "$$\n",
        "W_{h_i} = |c_i| / \\lambda\n",
        "$$\n",
        "\n",
        "**Otimizações fermiónicas e nativas do hardware**\n",
        "\n",
        "A função retorna `generate_preset_jw_pass_manager()` um objeto `MultiStagePassManager` que recebe um argumento e `FermionicCircuit` gera um circuito final otimizado, que podemos transpilá-lo para ser executado em nosso hardware. Substituímos sua etapa de otimização padrão por uma que `FermionicPassManager` contenha nossa passagem `QDriftTrotterization` :\n",
        "\n",
        "* A etapa `QDriftTrotterization` utiliza internamente o cálculo de pesos e a amostragem para gerar os circuitos que usaremos para a amostragem\n",
        "* A etapa `RelabelModes` é mais uma etapa de otimização que pode ser usada para permutar os modos fermiónicos, a fim de otimizar a conectividade entre os qubits e reduzir a profundidade das portas; leia mais na [referência](https://qiskit.github.io/qiskit-fermions/stable/0.1/stubs/qiskit_fermions.transpiler.passes.RelabelModes.html#qiskit_fermions.transpiler.passes.RelabelModes) da API\n",
        "\n",
        "As etapas restantes são executadas `MultiStagePassManager` automaticamente e cuidam de todo o mapeamento de férmions para qubits:\n",
        "\n",
        "* [F2QLayout](https://qiskit.github.io/qiskit-fermions/stable/0.1/stubs/qiskit_fermions.transpiler.passes.TrivialF2QLayout.html) : O gerenciador de passagens predefinido aplica a `TrivialF2QLayout` passagem, que mapeia de forma trivial os bits fermiônicos $n$ para os qubits $n$.\n",
        "* [F2QSynth](https://qiskit.github.io/qiskit-fermions/stable/0.1/stubs/qiskit_fermions.transpiler.passes.F2QSynthesis.html) : Uma etapa de transpilagem para mapear instruções de circuitos baseados em férmions para instruções baseadas em qubits.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "34dbca48",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "400\n"
          ]
        }
      ],
      "source": [
        "qdrift = QDriftTrotterization(num_exc, rng=19)\n",
        "\n",
        "pm = generate_preset_jw_pass_manager()\n",
        "pm.optimization = FermionicPassManager([qdrift])\n",
        "\n",
        "sqdrift_circuits = []\n",
        "for circ in init_circuits:\n",
        "    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))\n",
        "\n",
        "for circ in sqdrift_circuits:\n",
        "    circ.measure_all()\n",
        "\n",
        "print(len(sqdrift_circuits))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6f71b1cb",
      "metadata": {},
      "source": [
        "Agora que concluímos as otimizações no nível dos férmions, podemos transpilá-los para execução no simulador.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "de130514",
      "metadata": {},
      "outputs": [],
      "source": [
        "simulator = AerSimulator()\n",
        "shots = 100\n",
        "\n",
        "transpiled_circuits = transpile(sqdrift_circuits, simulator)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "98db4d80",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-qiskit-primitives\" />\n",
        "\n",
        "### Etapa 3: Executar usando o comando `Qiskit primitives`\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "74496e28",
      "metadata": {},
      "source": [
        "Agora que já temos nossos circuitos, podemos executá-los usando o Qiskit primitives no site AerSimulator. Vamos somar todas as contagens dos diferentes circuitos. Nós os convertemos em vetores booleanos antes de, por fim, realizar o pós-processamento com o SQD.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "24b9df39",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Executing 400 circuits with 100 shots each...\n",
            "400 length before post processing\n"
          ]
        }
      ],
      "source": [
        "print(\n",
        "    f\"Executing {len(transpiled_circuits)} circuits with {shots} shots each...\"\n",
        ")\n",
        "\n",
        "job = simulator.run(transpiled_circuits, shots=shots)\n",
        "result = job.result()\n",
        "\n",
        "all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]\n",
        "\n",
        "print(len(all_counts), \"length before post processing\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b797c220",
      "metadata": {},
      "source": [
        "<span id=\"step-4-post-process-and-return-result-in-desired-classical-format\" />\n",
        "\n",
        "### Etapa 4: Realizar o pós-processamento e apresentar o resultado no formato clássico desejado\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2129d0ac",
      "metadata": {},
      "source": [
        "**Utilização de cadeias de bits para SQD**\n",
        "\n",
        "Agora podemos aplicar o esquema de diagonalização às sequências de bits selecionadas para encontrar o menor valor próprio que corresponderá à energia do estado fundamental da molécula. Criamos uma função de retorno de chamada, declaramos as ocupações iniciais e definimos os parâmetros antes de, finalmente, executar o esquema de diagonalização. A função de retorno é usada para exibir a iteração atual e a estimativa atual do valor próprio a cada iteração.\n",
        "\n",
        "Por fim, para obter a estimativa do estado fundamental, somamos o `nuclear_repulsion_energy` à energia resultante.\n",
        "\n",
        "**Observação** : a dimensão do subespaço não é fixa entre as iterações, mesmo no simulador sem ruído — cada subamostra gera um conjunto diferente de configurações, e a etapa de recuperação reestrutura o conjunto entre as iterações; portanto, a dimensão relatada varia de uma subamostra para outra. A amostragem silenciosa, por si só, não determina a dimensão do subespaço selecionado. A execução no hardware, no entanto, tende a gerar subespaços sistematicamente maiores, pois as imagens com ruído quebram a simetria do número de partículas e a recuperação da configuração as transforma em vetores de base adicionais. Por isso, também apresentaremos outra etapa para a poda de cadeias de bits na seção de hardware.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "f7c5b2ed",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "40000\n",
            "  Alpha electrons: 7\n",
            "  Beta electrons: 7\n",
            "  Number of orbitals: 10\n",
            "  Number of spin orbitals (qubits): 20\n",
            "Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)\n",
            "\n",
            "Running SQD with configuration recovery...\n",
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -107.64767025226178\n",
            "\t\tSubspace dimension: 5538\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -107.64772799119115\n",
            "\t\tSubspace dimension: 5670\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -107.64765512281548\n",
            "\t\tSubspace dimension: 5767\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -107.64795948524682\n",
            "\t\tSubspace dimension: 6080\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -107.64806617355072\n",
            "\t\tSubspace dimension: 6300\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -107.64802260640258\n",
            "\t\tSubspace dimension: 6308\n",
            "FINAL SQD RESULTS\n",
            "Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712\n",
            " 0.99301927 0.0373282  0.0373266  0.00944508]\n",
            "Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737\n",
            " 0.99302145 0.03733536 0.03733399 0.0094585 ]\n",
            "Reference Energy: -107.6481842917 Ha\n",
            "Computed Energy:  -107.6480661736 Ha\n",
            "Error:            1.1811817564e-04 Ha\n"
          ]
        }
      ],
      "source": [
        "combined_counts = {}\n",
        "for counts in all_counts:\n",
        "    for bitstring, count in counts.items():\n",
        "        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count\n",
        "\n",
        "bit_array = BitArray.from_counts(combined_counts)\n",
        "print(bit_array.num_shots)\n",
        "\n",
        "print(f\"  Alpha electrons: {num_elec_a}\")\n",
        "print(f\"  Beta electrons: {num_elec_b}\")\n",
        "print(f\"  Number of orbitals: {norb}\")\n",
        "print(f\"  Number of spin orbitals (qubits): {2*norb}\")\n",
        "print(f\"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}\")\n",
        "\n",
        "# SQD parameters\n",
        "samples_per_batch = 300\n",
        "num_batches = 3\n",
        "max_iterations = 5\n",
        "\n",
        "initial_occupancies = (\n",
        "    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha\n",
        "    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta\n",
        ")\n",
        "\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 + nuclear_repulsion_energy}\")\n",
        "        print(\n",
        "            f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
        "        )\n",
        "\n",
        "\n",
        "# Run SQD with configuration recovery\n",
        "print(\"\\nRunning SQD with configuration recovery...\")\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    h1e,\n",
        "    h2e,\n",
        "    bit_array,\n",
        "    samples_per_batch=samples_per_batch,\n",
        "    norb=norb,\n",
        "    nelec=(num_elec_a, num_elec_b),\n",
        "    num_batches=num_batches,\n",
        "    energy_tol=1e-3,\n",
        "    occupancies_tol=1e-3,\n",
        "    max_iterations=max_iterations,\n",
        "    initial_occupancies=initial_occupancies,\n",
        "    seed=42,\n",
        "    callback=callback,\n",
        ")\n",
        "\n",
        "computed_energy = result.energy + nuclear_repulsion_energy\n",
        "\n",
        "print(\"FINAL SQD RESULTS\")\n",
        "print(f\"Orbital occupancies (alpha): {result.orbital_occupancies[0]}\")\n",
        "print(f\"Orbital occupancies (beta): {result.orbital_occupancies[1]}\")\n",
        "\n",
        "\n",
        "energy_error = abs(computed_energy - reference_energy)\n",
        "print(f\"Reference Energy: {reference_energy:.10f} Ha\")\n",
        "print(f\"Computed Energy:  {computed_energy:.10f} Ha\")\n",
        "print(f\"Error:            {energy_error:.10e} Ha\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "84705007",
      "metadata": {},
      "source": [
        "<span id=\"hardware-example\" />\n",
        "\n",
        "## Exemplo de hardware\n",
        "\n",
        "Este exemplo utiliza 20 qubits (10 orbitais espaciais). Essa escolha é uma medida de conveniência para um tutorial que deve ser executado rapidamente, e não um limite rígido para o método.\n",
        "\n",
        "O custo da etapa clássica não é determinado diretamente pelo número de qubits. O SQD diagonaliza o hamiltoniano projetado no subespaço gerado pelas configurações *amostradas*; portanto, o que determina o custo clássico é a dimensão desse subespaço selecionado — determinada, neste caso, por `samples_per_batch`, `num_batches`, e pelo número de configurações distintas que os circuitos realmente produzem — juntamente com a álgebra linear esparsa necessária para aplicar o hamiltoniano projetado. O espaço completo de CI cresce combinatoriamente com os orbitais e os elétrons, mas o subespaço selecionado é uma pequena parcela ajustável desse espaço, e controlamos seu tamanho diretamente. Consequentemente, o número de qubits e a dificuldade clássica podem ser variados de forma relativamente independente: um espaço orbital mais amplo, amostrado em um subespaço modesto, pode ser mais econômico do que um sistema menor diagonalizado em um espaço orbital muito grande.\n",
        "\n",
        "Na prática, portanto, o tamanho viável do sistema depende da dimensão do subespaço necessária para a precisão desejada e da memória e dos núcleos disponíveis para o solucionador de valores próprios. Espaços orbitais maiores geralmente exigem um subespaço maior para se alcançar precisão química, e é isso que, em última instância, motiva o uso de recursos distribuídos — consulte [o qiskit-addon-sqd-hpc](https://qiskit.github.io/qiskit-addon-sqd-hpc/) para saber como ampliar essa etapa. Em vez de definir um limite fixo, a abordagem prática consiste em acompanhar a dimensão do subespaço relatada e a convergência da energia ao longo das iterações, aumentando o tamanho do subespaço até que a energia pare de melhorar ou até esgotar a memória disponível.\n",
        "\n",
        "*Observação:* Devido ao erro de amostragem causado pelo ruído no hardware, o subespaço criado para a diagonalização na execução no hardware será maior do que o obtido ao usar o simulador. Embora aumente a dimensão do subespaço que desejamos diagonalizar, o fluxo de trabalho ainda nos fornece uma resposta precisa devido à robustez do SQD em relação ao ruído.\n",
        "\n",
        "**Eliminação de strings espúrias**\n",
        "\n",
        "Aqui, podemos optar por realizar uma etapa adicional. Quando tivermos todas as sequências de bits resultantes das execuções do circuito, podemos filtrar as sequências inválidas antes de executar o SQD ou prosseguir sem fazer a poda. Geralmente, é preferível não realizar a poda em execuções de hardware, pois isso mantém as tentativas com simetria quebrada disponíveis para a recuperação da configuração, que pode repará-las, transformando-as em configurações válidas e, assim, ampliar o subespaço, em vez de descartar essas tentativas imediatamente.\n",
        "\n",
        "Como o nitrogênio só pode ter sete $\\alpha$ e sete $\\beta$ elétrons, quaisquer sequências de bits que tenham mais ou menos do que sete 1s na primeira e na segunda metade da saída podem ser descartadas. Definimos uma função que verifica se as sequências de bits são válidas e, caso não sejam, as descarta. Depois de filtrarmos as sequências de bits espúrias, o restante é enviado para o esquema de diagonalização. Use o sinalizador `PRUNE` abaixo para alternar entre os dois comportamentos.\n",
        "\n",
        "Lembre-se de que a poda é apenas uma das várias opções que definem o subespaço final, juntamente com o número de circuitos, o conjunto de tempos de evolução e a filtragem de termos diagonais. Comparar uma execução podada com uma não podada só é informativo se todos os outros fatores forem mantidos constantes; o [guia complementar do C++](https://github.com/Qiskit/documentation/tree/main/docs/tutorials/assets/sqdrift/cpp) aborda esse assunto com mais detalhes, já que ele utiliza pós-seleção em vez de recuperação e também difere nesses outros parâmetros.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8276be0a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Parsing assets/sqdrift/fcidump_files/N2_sto_3g\n",
            "Reference FCI Energy  = -107.6481842917 Ha\n",
            "Nuclear Repulsion Energy = 23.7887003074 Ha\n",
            "5060\n",
            "400\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": []
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Selected backend: ibm_aachen (156 qubits)\n",
            "40000\n",
            "Electron configuration:\n",
            "  Total electrons: 14\n",
            "  Alpha electrons: 7\n",
            "  Beta electrons: 7\n",
            "  Number of orbitals: 10\n",
            "  Number of spin orbitals (qubits): 20\n",
            "Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)\n",
            "\n",
            "Running SQD with configuration recovery...\n",
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -107.64593072647523\n",
            "\t\tSubspace dimension: 7221\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -107.6458270048177\n",
            "\t\tSubspace dimension: 7209\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -107.64007673117075\n",
            "\t\tSubspace dimension: 7138\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -107.64757372124944\n",
            "\t\tSubspace dimension: 9009\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -107.64674060104392\n",
            "\t\tSubspace dimension: 8245\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -107.64731360491942\n",
            "\t\tSubspace dimension: 8178\n",
            "Iteration 3\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -107.64765518770588\n",
            "\t\tSubspace dimension: 8835\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -107.64767975712016\n",
            "\t\tSubspace dimension: 8649\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -107.64761634415606\n",
            "\t\tSubspace dimension: 8648\n",
            "FINAL SQD RESULTS\n",
            "Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295\n",
            " 0.99298797 0.03728154 0.03728186 0.00938359]\n",
            "Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174\n",
            " 0.99298424 0.03728121 0.03728169 0.00939159]\n",
            "Reference Energy: -107.6481842917 Ha\n",
            "Computed Energy:  -107.6476797571 Ha\n",
            "Error:            5.0453460619e-04 Ha\n"
          ]
        }
      ],
      "source": [
        "name = \"assets/sqdrift/fcidump_files/N2_sto_3g\"\n",
        "\n",
        "fcidump = tools.fcidump.read(name)\n",
        "\n",
        "# Extract metadata from the FCIDump header\n",
        "norb = fcidump[\"NORB\"]  # number of spatial orbitals\n",
        "nelec = fcidump[\"NELEC\"]  # total number of electrons\n",
        "e_nuc = fcidump[\"ECORE\"]  # nuclear repulsion / core energy\n",
        "ms2 = fcidump[\"MS2\"]  # 2S (spin)\n",
        "\n",
        "num_elec_a = (nelec + ms2) // 2  # alpha electrons\n",
        "num_elec_b = (nelec - ms2) // 2  # beta  electrons\n",
        "\n",
        "# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)\n",
        "h1e = fcidump[\"H1\"]  # shape (norb, norb)\n",
        "h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)\n",
        "    1, fcidump[\"H2\"], norb\n",
        ")\n",
        "\n",
        "cisolver = fci.direct_spin1.FCI()\n",
        "cisolver.max_cycle = 200\n",
        "cisolver.conv_tol = 1e-12\n",
        "\n",
        "e_fci, _ = cisolver.kernel(\n",
        "    h1e,\n",
        "    h2e,\n",
        "    norb,\n",
        "    (num_elec_a, num_elec_b),\n",
        "    ecore=e_nuc,  # adds nuclear repulsion to the final energy\n",
        ")\n",
        "\n",
        "reference_energy = e_fci\n",
        "\n",
        "print(f\"Reference FCI Energy  = {reference_energy:.10f} Ha\")\n",
        "\n",
        "nuclear_repulsion_energy = fcidump[\"ECORE\"]\n",
        "print(f\"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha\")\n",
        "\n",
        "fcidump = FCIDump.from_file(name)\n",
        "hamiltonian = FermionOperator.from_fcidump(fcidump)\n",
        "num_modes = 2 * fcidump.norb\n",
        "\n",
        "# Apply automatic grouping\n",
        "canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))\n",
        "exit_code = group_terms_by_electronic_structure(\n",
        "    canon, num_modes, two_body_physicist_order=False\n",
        ")\n",
        "filter_diagonal_terms(canon)\n",
        "\n",
        "print(len(canon.groups))\n",
        "\n",
        "# SqDRIFT parameters\n",
        "times = [1.0, 10.0]  # Total evolution times used for the subspace creation\n",
        "num_exc = 10  # Number of excitation groups per circuit\n",
        "num_circuits = 200  # Number of circuits to generate\n",
        "\n",
        "init_circuits = []\n",
        "hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))\n",
        "\n",
        "for time in times:\n",
        "    evo_gate = Evolution(num_modes, canon, time)\n",
        "    circ = FermionicCircuit(num_modes)\n",
        "    circ.append(hf_gate, circ.modes)\n",
        "    circ.append(evo_gate, circ.modes)\n",
        "    init_circuits.append(circ)\n",
        "\n",
        "# Calculate weights for sampling (one per group)\n",
        "qdrift = QDriftTrotterization(num_exc, rng=19)\n",
        "\n",
        "pm = generate_preset_jw_pass_manager()\n",
        "pm.optimization = FermionicPassManager([qdrift])\n",
        "\n",
        "sqdrift_circuits = []\n",
        "for circ in init_circuits:\n",
        "    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))\n",
        "\n",
        "for circ in sqdrift_circuits:\n",
        "    circ.measure_all()\n",
        "\n",
        "print(len(sqdrift_circuits))\n",
        "\n",
        "# This example assumes you have saved your IBM Quantum Platform account locally.\n",
        "service = QiskitRuntimeService(channel=\"ibm_quantum_platform\")\n",
        "\n",
        "# Select backend (choose based on qubit requirements)\n",
        "backend = service.least_busy(\n",
        "    operational=True,\n",
        "    simulator=False,\n",
        "    min_num_qubits=2 * norb,\n",
        ")\n",
        "\n",
        "print(f\"Selected backend: {backend.name} ({backend.num_qubits} qubits)\")\n",
        "\n",
        "# Transpile for hardware\n",
        "transpiled_circuits = transpile(\n",
        "    sqdrift_circuits,\n",
        "    backend=backend,\n",
        "    optimization_level=3,\n",
        "    seed_transpiler=42,\n",
        ")\n",
        "\n",
        "shots = 100\n",
        "\n",
        "sampler = Sampler(mode=backend)\n",
        "\n",
        "sampler.options.environment.job_tags = [\"TUT-SqDRIFT\"]\n",
        "\n",
        "job = sampler.run(transpiled_circuits, shots=shots)\n",
        "result = job.result()\n",
        "\n",
        "# Extract counts from SamplerV2 results\n",
        "all_counts = [pub_result.data.meas.get_counts() for pub_result in result]\n",
        "\n",
        "# Set to True to filter out bitstrings that violate electron-number conservation\n",
        "PRUNE = False\n",
        "\n",
        "\n",
        "def is_valid_bitstring(\n",
        "    bitstring: str, norb: int, nelec: tuple[int, int]\n",
        ") -> bool:\n",
        "    n_alpha, n_beta = nelec\n",
        "    return (\n",
        "        len(bitstring) == 2 * norb\n",
        "        and bitstring[norb:].count(\"1\") == n_alpha\n",
        "        and bitstring[:norb].count(\"1\") == n_beta\n",
        "    )\n",
        "\n",
        "\n",
        "if PRUNE:\n",
        "    all_counts_filtered = []\n",
        "    for counts in all_counts:\n",
        "        filtered_count = {}\n",
        "        for key in counts:\n",
        "            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):\n",
        "                continue\n",
        "            elif key not in filtered_count.keys():\n",
        "                filtered_count[key] = counts[key]\n",
        "            else:\n",
        "                filtered_count[key] += counts[key]\n",
        "        all_counts_filtered.append(filtered_count)\n",
        "    all_counts = all_counts_filtered\n",
        "\n",
        "combined_counts = {}\n",
        "for counts in all_counts:\n",
        "    for bitstring, count in counts.items():\n",
        "        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count\n",
        "\n",
        "bit_array = BitArray.from_counts(combined_counts)\n",
        "print(bit_array.num_shots)\n",
        "\n",
        "print(\"Electron configuration:\")\n",
        "print(f\"  Total electrons: {nelec}\")\n",
        "print(f\"  Alpha electrons: {num_elec_a}\")\n",
        "print(f\"  Beta electrons: {num_elec_b}\")\n",
        "print(f\"  Number of orbitals: {norb}\")\n",
        "print(f\"  Number of spin orbitals (qubits): {2*norb}\")\n",
        "print(f\"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}\")\n",
        "\n",
        "# SQD parameters\n",
        "samples_per_batch = 300\n",
        "num_batches = 3\n",
        "max_iterations = 5\n",
        "\n",
        "initial_occupancies = (\n",
        "    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha\n",
        "    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta\n",
        ")\n",
        "\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 + nuclear_repulsion_energy}\")\n",
        "        print(\n",
        "            f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
        "        )\n",
        "\n",
        "\n",
        "# Run SQD with configuration recovery\n",
        "print(\"\\nRunning SQD with configuration recovery...\")\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    h1e,\n",
        "    h2e,\n",
        "    bit_array,\n",
        "    samples_per_batch=samples_per_batch,\n",
        "    norb=norb,\n",
        "    nelec=(num_elec_a, num_elec_b),\n",
        "    num_batches=num_batches,\n",
        "    energy_tol=1e-3,\n",
        "    occupancies_tol=1e-3,\n",
        "    max_iterations=max_iterations,\n",
        "    initial_occupancies=initial_occupancies,\n",
        "    seed=42,\n",
        "    callback=callback,\n",
        ")\n",
        "\n",
        "computed_energy = result.energy + nuclear_repulsion_energy\n",
        "\n",
        "print(\"FINAL SQD RESULTS\")\n",
        "print(f\"Orbital occupancies (alpha): {result.orbital_occupancies[0]}\")\n",
        "print(f\"Orbital occupancies (beta): {result.orbital_occupancies[1]}\")\n",
        "\n",
        "\n",
        "energy_error = abs(computed_energy - reference_energy)\n",
        "print(f\"Reference Energy: {reference_energy:.10f} Ha\")\n",
        "print(f\"Computed Energy:  {computed_energy:.10f} Ha\")\n",
        "print(f\"Error:            {energy_error:.10e} Ha\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e18de7d4",
      "metadata": {},
      "source": [
        "<span id=\"next-steps\" />\n",
        "\n",
        "## Próximas etapas\n",
        "\n",
        "<Admonition type=\"note\" title=\"Recomendações\">\n",
        "  Se você achou este trabalho interessante, talvez se interesse pelo material a seguir:\n",
        "\n",
        "  * [Diagonalização quântica de Krylov baseada em amostras de um modelo de rede fermiónica](/docs/tutorials/sample-based-krylov-quantum-diagonalization) — um tutorial relacionado que utiliza circuitos de evolução temporal em vez de um ansatz variacional.\n",
        "  * [Diagonalização quântica baseada em amostras de um hamiltoniano químico](/docs/tutorials/sample-based-quantum-diagonalization) — um tutorial sobre como construir um circuito Jastrow de cluster unitário local (LUCJ) para simulação em química quântica.\n",
        "  * O artigo “ [SqDRIFT](https://arxiv.org/abs/2508.02578) ” — a literatura na qual este tutorial se baseia. (Observe que algumas das otimizações discutidas neste artigo ainda estão em desenvolvimento, e este tutorial está sujeito a alterações no futuro, de acordo com a evolução das bibliotecas utilizadas.)\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}