{
"cells": [
{
"cell_type": "markdown",
"id": "9e40af77-7f0f-4dd6-ab0a-420cf396050e",
"metadata": {},
"source": [
"---\n",
"title: \"Diagonalização quantum baseada em amostras de um Hamiltoniano químico\"\n",
"description: \"Use o algoritmo de diagonalização quântica baseado em amostras para simular uma molécula de nitrogênio usando hardware quântico ruidoso.\"\n",
"---\n",
"\n",
"{/* cspell:ignore fontdict fontsize milli LUCJ CCSD ccsd hcore pvdz */}\n",
"\n",
"\n",
"\n",
"# Diagonalização quantum baseada em amostras de um Hamiltoniano químico\n",
"\n",
"*Estimativa de uso: menos de um minuto em um processador Heron r2 (OBSERVAÇÃO: esta é apenas uma estimativa. Seu tempo de execução pode variar)*\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9701917b",
"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](/docs/addons/qiskit-addon-sqd) para aproximar a energia do estado fundamental de um sistema molecular utilizando sequências de bits amostradas a partir de uma unidade de processamento quântico (QPU).\n",
"* Como usar [o ffsim](https://github.com/qiskit-community/ffsim) para construir um circuito Jastrow de cluster unitário local (LUCJ) para simulação em química quântica.\n",
"\n",
"\n",
"\n",
"## Pré-requisitos\n",
"\n",
"Recomendamos que os usuários se familiarizem com os seguintes tópicos antes de seguir com este tutorial:\n",
"\n",
"* Química quântica e segunda quantização\n",
"* Utilizando a primitiva Sampler para obter amostras de circuitos quânticos\n",
"\n",
"\n",
"\n",
"## Segundo plano\n",
"\n",
"Neste tutorial, mostramos como realizar o pós-processamento de amostras quânticas com ruído para aproximar o estado fundamental da molécula de nitrogênio $\\text{N}_2$ no comprimento de ligação de equilíbrio, utilizando o [complemento SQD do Qiskit](https://github.com/Qiskit/qiskit-addon-sqd) para implementar [o algoritmo de diagonalização quântica baseada em amostras (SQD)](https://arxiv.org/abs/2405.05068). Mais detalhes sobre o software podem ser encontrados na [documentação](/docs/addons/qiskit-addon-sqd) correspondente, incluindo um [exemplo simples](/docs/addons/qiskit-addon-sqd/guides/quickstart) para começar.\n",
"\n",
"Este tutorial é recomendado para usuários familiarizados com a química quântica: mais especificamente, com o cálculo das energias do estado fundamental de uma molécula. Para obter um guia detalhado sobre o fluxo de trabalho, consulte o [curso](/learning/courses/quantum-diagonalization-algorithms) sobre o algoritmo de diagonalização quântica.\n",
"\n",
"A SQD é uma técnica para determinar os autovalores e autovetores de operadores quânticos, como o hamiltoniano de um sistema quântico, por meio da combinação da computação quântica com a computação clássica distribuída. A computação distribuída clássica é utilizada para processar amostras obtidas de um processador quântico e para projetar e diagonalizar um hamiltoniano alvo num subespaço que essas amostras definem. Um fluxo de trabalho baseado em SQD segue as seguintes etapas:\n",
"\n",
"1. Escolha um ansatz de circuito e aplique-o em um computador quântico a um estado de referência (nesse caso, o estado [Hartree-Fock](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method) ).\n",
"2. Amostra de cadeias de bits do estado quântico resultante.\n",
"3. Execute o procedimento *de recuperação de configuração autoconsistente* nas cadeias de bits para obter a aproximação do estado fundamental.\n",
"\n",
"Sabe-se que a SQD funciona bem quando o estado próprio alvo é esparso: a função de onda é suportada em um conjunto de estados básicos $\\mathcal{S} = \\{|x\\rangle \\}$ cujo tamanho não aumenta exponencialmente com o tamanho do problema.\n",
"\n",
"\n",
"\n",
"### Química quântica\n",
"\n",
"O Hamiltoniano de um sistema molecular pode ser escrito como\n",
"\n",
"$$\n",
"\\hat{H} = \\sum_{ \\substack{pr\\\\\\sigma} } h_{pr} \\, \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{r\\sigma}\n",
"+ \\frac12\n",
"\\sum_{ \\substack{prqs\\\\\\sigma\\tau} }\n",
"h_{prqs} \\,\n",
"\\hat{a}^\\dagger_{p\\sigma}\n",
"\\hat{a}^\\dagger_{q\\tau}\n",
"\\hat{a}_{s\\tau}\n",
"\\hat{a}_{r\\sigma},\n",
"$$\n",
"\n",
"onde $h_{pr}$ e $h_{prqs}$ são números complexos chamados de integrais moleculares que podem ser calculados a partir da especificação da molécula usando um programa de computador. Neste tutorial, calculamos as integrais usando o pacote de software [PySCF](https://pyscf.org/) pacote de software.\n",
"\n",
"Para obter detalhes sobre como o Hamiltoniano molecular é derivado, consulte um livro-texto sobre química quântica (por exemplo, *Modern Quantum Chemistry*, de Szabo e Ostlund). Para obter uma explicação de alto nível sobre como os problemas de química quântica são mapeados em computadores quânticos, confira a palestra [*Mapping Problems to Qubits (Mapeando problemas para Qubits*](https://youtube.com/watch?v=TyFU6r8uEsE\\&t=900) ) da Qiskit Global Summer School 2024.\n",
"\n",
"\n",
"\n",
"### Abordagem do cluster unitário local de Jastrow (LUCJ)\n",
"\n",
"O SQD requer um ansatz de circuito quântico do qual extrair amostras. Neste tutorial, utilizaremos o modelo [do cluster unitário local de Jastrow (LUCJ),](https://pubs.rsc.org/en/content/articlelanding/2023/sc/d3sc02516k) devido à sua combinação de fundamentação física e facilidade de implementação em hardware. Usaremos [o ffsim](https://qiskit-community.github.io/ffsim/) para construir o circuito de ansatz.\n",
"\n",
"A abordagem LUCJ adapta-se a QPUs com conectividade de qubits restrita. Os orbitais de spin são mapeados para qubits de forma que a hipótese não exija o roteamento por meio de portas SWAP. IBM® O hardware possui uma topologia de qubits em rede hexagonal densa; nesse caso, podemos adotar um padrão em “ziguezague”, ilustrado abaixo. Nesse padrão, os orbitais com o mesmo spin são mapeados para qubits com topologia em linha (círculos vermelhos e azuis), e há uma conexão entre orbitais de spin diferente a cada quarto orbital espacial, sendo essa conexão facilitada por um qubit auxiliar (círculos roxos).\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"### Recuperação de configuração auto-consistente\n",
"\n",
"O procedimento de recuperação de configuração autoconsistente foi projetado para extrair o máximo de sinal possível de amostras quânticas com ruído. Como o Hamiltoniano molecular conserva o número de partículas e o spin Z, faz sentido escolher um ansatz de circuito que também conserve essas simetrias. Quando aplicado ao estado Hartree-Fock, o estado resultante tem um número de partículas e um spin Z fixos na configuração sem ruído. Portanto, as metades spin- $\\alpha$ e spin- $\\beta$ de qualquer bitstring amostrada a partir desse estado devem ter o mesmo [peso de Hamming](https://en.wikipedia.org/wiki/Hamming_weight) que no estado Hartree-Fock. Devido à presença de ruído nos processadores quânticos atuais, algumas cadeias de bits medidas violarão essa propriedade. Uma forma simples de pós-seleção descartaria essas cadeias de bits, mas isso é um desperdício porque essas cadeias de bits ainda podem conter algum sinal. O procedimento de recuperação autoconsistente tenta recuperar parte desse sinal no pós-processamento. O procedimento é iterativo e requer, como entrada, uma estimativa das ocupações médias de cada orbital no estado fundamental, que é primeiro computada a partir das amostras brutas. O procedimento é executado em um loop, e cada iteração tem as seguintes etapas:\n",
"\n",
"1. Para cada bitstring que violar as simetrias especificadas, inverta seus bits com um procedimento probabilístico projetado para aproximar o bitstring da estimativa atual das ocupações orbitais médias, para obter um novo bitstring.\n",
"2. Coletar todas as cadeias de bits antigas e novas que satisfaçam as simetrias e subamostrar subconjuntos de tamanho fixo, escolhidos antecipadamente.\n",
"3. Para cada subconjunto de cadeias de bits, projete o Hamiltoniano no subespaço abrangido pelos vetores de base correspondentes (consulte a [seção anterior](#quantum-chemistry) para obter uma descrição desses vetores de base) e calcule uma estimativa do estado fundamental do Hamiltoniano projetado em um computador clássico.\n",
"4. Atualize a estimativa das ocupações orbitais médias com a estimativa do estado fundamental com a energia mais baixa.\n",
"\n",
"\n",
"\n",
"### Diagrama do fluxo de trabalho SQD\n",
"\n",
"O fluxo de trabalho do SQD está representado no diagrama a seguir:\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "88422c4b",
"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.75 ou mais recente (`pip install ffsim`)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "c6e44a31",
"metadata": {},
"source": [
"\n",
"\n",
"## Instalação\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "6e51c3d8",
"metadata": {},
"outputs": [],
"source": [
"import math\n",
"\n",
"import ffsim\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pyscf\n",
"import pyscf.cc\n",
"import pyscf.mcscf\n",
"from qiskit import QuantumCircuit, QuantumRegister\n",
"from qiskit.primitives import StatevectorSampler\n",
"from qiskit.providers.fake_provider import GenericBackendV2\n",
"from qiskit_ibm_runtime import QiskitRuntimeService\n",
"from qiskit_ibm_runtime import SamplerV2 as Sampler"
]
},
{
"cell_type": "markdown",
"id": "4bc6ee26-4371-4cd2-80a7-60752bf8775d",
"metadata": {},
"source": [
"\n",
"\n",
"## Exemplo de simulador em pequena escala\n",
"\n",
"Neste tutorial, vamos encontrar uma aproximação do estado fundamental de uma molécula de nitrogênio próxima à sua distância de ligação de equilíbrio. Primeiramente, utilizamos um pequeno conjunto de bases do tipo “ STO-6G ” para simular o experimento e garantir que tudo esteja funcionando corretamente.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "afeb054c",
"metadata": {},
"source": [
"\n",
"\n",
"### Passo 1: Mapear entradas clássicas para um problema quântico\n",
"\n",
"Primeiro, especificamos a molécula e suas propriedades.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "b821e660",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"converged SCF energy = -108.464957764796\n",
"CASCI E = -108.595987350986 E(CI) = -32.4115475088426 S^2 = 0.0000000\n",
"norb = 8\n",
"nelec = (5, 5)\n"
]
}
],
"source": [
"# Specify molecule properties\n",
"spin_sq = 0\n",
"\n",
"# Build N2 molecule\n",
"mol = pyscf.gto.Mole()\n",
"mol.build(\n",
" atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]],\n",
" basis=\"sto-6g\",\n",
" symmetry=\"Dooh\",\n",
")\n",
"\n",
"# Define active space\n",
"n_frozen = 2\n",
"active_space = range(n_frozen, mol.nao_nr())\n",
"\n",
"# Get molecular integrals\n",
"scf = pyscf.scf.RHF(mol).run()\n",
"norb = len(active_space)\n",
"n_electrons = int(sum(scf.mo_occ[active_space]))\n",
"n_alpha = (n_electrons + mol.spin) // 2\n",
"n_beta = (n_electrons - mol.spin) // 2\n",
"nelec = (n_alpha, n_beta)\n",
"cas = pyscf.mcscf.CASCI(scf, norb, nelec)\n",
"mo = cas.sort_mo(active_space, base=0)\n",
"hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n",
"eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)\n",
"\n",
"# Compute exact energy using FCI\n",
"reference_energy = cas.run().e_tot\n",
"\n",
"print(f\"norb = {norb}\")\n",
"print(f\"nelec = {nelec}\")"
]
},
{
"cell_type": "markdown",
"id": "96bfe018",
"metadata": {},
"source": [
"Antes de construir o circuito ansatz LUCJ, primeiro realizamos um cálculo CCSD na seguinte célula de código. As [amplitudes $t_1$ e $t_2$](https://en.wikipedia.org/wiki/Coupled_cluster#Cluster_operator) desse cálculo serão usadas para inicializar os parâmetros do ansatz.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "efe83d98",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"E(CCSD) = -108.5933309085008 E_corr = -0.1283731437052354\n"
]
}
],
"source": [
"# Get CCSD t2 amplitudes for initializing the ansatz\n",
"ccsd = pyscf.cc.CCSD(\n",
" scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]\n",
").run()\n",
"t1 = ccsd.t1\n",
"t2 = ccsd.t2"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "f4d882fa",
"metadata": {},
"source": [
"Agora, usamos [o ffsim](https://github.com/qiskit-community/ffsim) para criar o circuito de ansatz. Como nossa molécula possui um estado de Hartree-Fock de camada fechada, utilizamos a variante com equilíbrio de spin do ansatz UCJ, [UCJOpSpinBalanced](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced). Definimos `optimize=True` no método `from_t_amplitudes` para habilitar a dupla fatoração \"comprimida\" das amplitudes de $t_2$ (consulte [a abordagem do conjunto unitário local de Jastrow (LUCJ)](https://qiskit-community.github.io/ffsim/explanations/lucj.html#Parameter-initialization-from-CCSD) na documentação do ffsim para obter mais detalhes).\n",
"\n",
"Como o ansatz LUCJ se adapta à conectividade disponível da QPU, precisamos inicializar o backend da QPU antes de criar o ansatz. Por enquanto, vamos criar um backend genérico com um mapa de acoplamento hexagonal forte e um conjunto de portas no qual o ansatz LUCJ se decompõe naturalmente. Em seguida, usaremos `ffsim.qiskit.generate_lucj_pass_manager` para criar um gerenciador de passagens especializado em transpilá-lo o ansatz LUCJ para o backend especificado, de acordo com o layout “zig-zag” descrito na [seção](#local-unitary-cluster-jastrow-lucj-ansatz) de contextualização sobre o ansatz LUCJ. Esta função utiliza uma heurística de pontuação para minimizar os erros associados ao layout selecionado, o que é importante se o seu backend for uma QPU real ou um simulador com um modelo de ruído. Além de retornar o gerenciador de passagem, esta função também retorna os pares de acoplamento alfa-beta que podem ser implementados no hardware. Se nem todos os pares puderem ser implementados, será exibido um aviso.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dd69a86c",
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"\n",
"from qiskit.transpiler import CouplingMap\n",
"\n",
"warnings.formatwarning = lambda msg, *args, **kwargs: f\"Warning: {msg}\\n\"\n",
"\n",
"# Set ansatz properties\n",
"n_reps = 1\n",
"pairs_aa = [(p, p + 1) for p in range(norb - 1)]\n",
"\n",
"# Let generate_lucj_pass_manager determine the alpha-beta interactions\n",
"pairs_ab = None\n",
"\n",
"# Initialize backend\n",
"coupling_map = CouplingMap.from_heavy_hex(3)\n",
"backend = GenericBackendV2(\n",
" coupling_map.size(),\n",
" coupling_map=coupling_map,\n",
" basis_gates=[\"cp\", \"xx_plus_yy\", \"p\", \"x\", \"swap\"],\n",
")\n",
"\n",
"# Create pass manager\n",
"pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(\n",
" backend=backend,\n",
" norb=norb,\n",
" connectivity=\"heavy-hex\",\n",
" interaction_pairs=(pairs_aa, pairs_ab),\n",
" optimization_level=3,\n",
")\n",
"\n",
"# Create the LUCJ ansatz operator\n",
"ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
" t2=t2,\n",
" t1=t1,\n",
" n_reps=n_reps,\n",
" interaction_pairs=(pairs_aa, pairs_ab),\n",
" # Setting optimize=True enables the \"compressed\" factorization\n",
" optimize=True,\n",
" # Limit the number of optimization iterations to prevent the code cell\n",
" # from running too long. Removing this line may improve results.\n",
" options=dict(maxiter=1000),\n",
")\n",
"\n",
"# create an empty quantum circuit\n",
"qubits = QuantumRegister(2 * norb, name=\"q\")\n",
"circuit = QuantumCircuit(qubits)\n",
"\n",
"# prepare Hartree-Fock state as the reference state and append it\n",
"# to the quantum circuit\n",
"circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)\n",
"\n",
"# apply the UCJ operator to the reference state\n",
"circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n",
"circuit.measure_all()"
]
},
{
"cell_type": "markdown",
"id": "db11bf6d",
"metadata": {},
"source": [
"\n",
"\n",
"### Passo 2: Otimizar para execução em hardware quântico\n",
"\n",
"Em seguida, otimizamos o circuito para um hardware específico. Normalmente, essa etapa envolve a inicialização do backend de hardware e de um gerenciador de passagens para esse backend. No entanto, como a abordagem LUCJ está adaptada à conectividade do hardware, já realizamos essas ações na etapa anterior. Tudo o que resta fazer é executar o gerenciador de passagens no circuito para compilá-lo para um circuito ISA que possa ser executado diretamente na QPU.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "7d554aa5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Gate counts: OrderedDict({'xx_plus_yy': 86, 'p': 16, 'measure': 16, 'cp': 15, 'x': 10, 'swap': 2, 'barrier': 1})\n"
]
}
],
"source": [
"isa_circuit = pass_manager.run(circuit)\n",
"print(f\"Gate counts: {isa_circuit.count_ops()}\")"
]
},
{
"cell_type": "markdown",
"id": "0cc1edef",
"metadata": {},
"source": [
"\n",
"\n",
"### Passo 3: Execute usando Qiskit primitives\n",
"\n",
"Depois de otimizar o circuito para execução em hardware, estamos prontos para executá-lo no hardware de destino e coletar amostras para a estimativa de energia do estado fundamental. Como temos apenas um circuito, usaremos [o modo de execução de trabalho](/docs/guides/execution-modes) do site Qiskit Runtime e executaremos nosso circuito.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "93c1cef3-298e-4deb-8512-769fe94cd5a5",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Warning: Trying to add QuantumRegister to a QuantumCircuit having a layout\n"
]
}
],
"source": [
"rng = np.random.default_rng()\n",
"sampler = StatevectorSampler(seed=rng)\n",
"job = sampler.run([isa_circuit], shots=100_000)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "332ecab3-77e6-473f-b0e7-af30f983393a",
"metadata": {},
"outputs": [],
"source": [
"primitive_result = job.result()\n",
"pub_result = primitive_result[0]"
]
},
{
"cell_type": "markdown",
"id": "6df05b6e",
"metadata": {},
"source": [
"\n",
"\n",
"### Etapa 4: Pós-processamento e retorno do resultado no formato clássico desejado\n",
"\n",
"Uma métrica útil para avaliar a qualidade da saída da QPU é o número de configurações válidas retornadas. Uma configuração válida tem o número correto de partículas e spin Z, o que significa que a metade direita da cadeia de bits tem peso de Hamming igual ao número de elétrons com spin para cima, e a metade esquerda tem peso de Hamming igual ao número de elétrons com spin para baixo. A célula a seguir calcula a fração de configurações amostradas que são válidas.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "718f8517",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fraction of sampled configurations that are valid: 1.0\n"
]
}
],
"source": [
"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",
"bit_array = pub_result.data.meas\n",
"num_valid = sum(\n",
" is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n",
")\n",
"valid_fraction = num_valid / bit_array.num_shots\n",
"print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")"
]
},
{
"cell_type": "markdown",
"id": "7b126f3a",
"metadata": {},
"source": [
"Todas as sequências de bits são válidas porque estamos fazendo a amostragem do circuito em um simulador sem ruído. Ao executar em uma QPU ruidosa, a fração será menor que um, mas esperamos que seja maior do que a fração que se esperaria se as sequências de bits fossem amostradas de forma aleatória e uniforme, o que é calculado na célula a seguir.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6b3e4bca",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625\n"
]
}
],
"source": [
"expected_fraction_random = (\n",
" math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n",
")\n",
"print(\n",
" f\"Expected fraction of valid configurations from uniformly random bitstrings: \"\n",
" f\"{expected_fraction_random}\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "eb704101-0fe8-4d12-b572-b1d844e35a90",
"metadata": {},
"source": [
"Agora, estimamos a energia do estado fundamental do Hamiltoniano usando a função `diagonalize_fermionic_hamiltonian` . Essa função executa o procedimento de recuperação de configuração autoconsistente para refinar iterativamente as amostras quânticas com ruído para melhorar a estimativa de energia. Passamos uma função de retorno de chamada para que possamos salvar os resultados intermediários para análise posterior. Consulte [a documentação da API](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) para obter explicações sobre os argumentos para `diagonalize_fermionic_hamiltonian`.\n",
"\n",
"Aqui, usamos o `initial_occupancies` argumento para `diagonalize_fermionic_hamiltonian` especificar a configuração de Hartree-Fock como a estimativa inicial para as ocupações orbitais no estado fundamental. Essa abordagem é sensata para sistemas em que o estado fundamental tem suporte significativo na configuração Hartree-Fock, mas pode não ser apropriada em outras situações, embora métodos computacionais mais avançados possam produzir melhores estimativas iniciais nesses casos. Especificar `initial_occupancies` também permite que a recuperação da configuração seja executada mesmo que nenhuma configuração válida tenha sido amostrada, como pode ser o caso ao amostrar um circuito grande em um QPU ruidoso. Sem esse argumento, a recuperação da configuração falharia e geraria um erro se nenhuma configuração válida fosse fornecida.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "2f32a352",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iteration 1\n",
"\tSubsample 0\n",
"\t\tEnergy: -108.59275573641656\n",
"\t\tSubspace dimension: 900\n",
"\tSubsample 1\n",
"\t\tEnergy: -108.59275573641656\n",
"\t\tSubspace dimension: 900\n",
"\tSubsample 2\n",
"\t\tEnergy: -108.59275573641656\n",
"\t\tSubspace dimension: 900\n",
"Iteration 2\n",
"\tSubsample 0\n",
"\t\tEnergy: -108.59275573641656\n",
"\t\tSubspace dimension: 900\n",
"\tSubsample 1\n",
"\t\tEnergy: -108.59275573641656\n",
"\t\tSubspace dimension: 900\n",
"\tSubsample 2\n",
"\t\tEnergy: -108.59275573641656\n",
"\t\tSubspace dimension: 900\n",
"Final energy: -108.59275573641656\n",
"Final energy error: 0.0032316145694579745\n"
]
}
],
"source": [
"from functools import partial\n",
"\n",
"from qiskit_addon_sqd.fermion import (\n",
" SCIResult,\n",
" diagonalize_fermionic_hamiltonian,\n",
" solve_sci_batch,\n",
")\n",
"\n",
"# SQD options\n",
"energy_tol = 1e-3\n",
"occupancies_tol = 1e-3\n",
"max_iterations = 5\n",
"\n",
"# Eigenstate solver options\n",
"num_batches = 3\n",
"samples_per_batch = 300\n",
"symmetrize_spin = True\n",
"carryover_threshold = 1e-4\n",
"max_cycle = 200\n",
"\n",
"# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies\n",
"initial_occupancies = (\n",
" np.array([1] * n_alpha + [0] * (norb - n_alpha)),\n",
" np.array([1] * n_beta + [0] * (norb - n_beta)),\n",
")\n",
"\n",
"# Pass options to the built-in eigensolver. If you just want to use the defaults,\n",
"# you can omit this step, in which case you would not specify the sci_solver argument\n",
"# in the call to diagonalize_fermionic_hamiltonian below.\n",
"sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)\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 + nuclear_repulsion_energy}\")\n",
" print(\n",
" f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
" )\n",
"\n",
"\n",
"result = diagonalize_fermionic_hamiltonian(\n",
" hcore,\n",
" eri,\n",
" bit_array,\n",
" samples_per_batch=samples_per_batch,\n",
" norb=norb,\n",
" nelec=nelec,\n",
" num_batches=num_batches,\n",
" energy_tol=energy_tol,\n",
" occupancies_tol=occupancies_tol,\n",
" max_iterations=max_iterations,\n",
" sci_solver=sci_solver,\n",
" symmetrize_spin=symmetrize_spin,\n",
" initial_occupancies=initial_occupancies,\n",
" carryover_threshold=carryover_threshold,\n",
" callback=callback,\n",
" seed=rng,\n",
")\n",
"\n",
"final_energy = result.energy + nuclear_repulsion_energy\n",
"energy_error = final_energy - reference_energy\n",
"print(f\"Final energy: {final_energy}\")\n",
"print(f\"Final energy error: {energy_error}\")"
]
},
{
"cell_type": "markdown",
"id": "9d78906b-4759-4506-9c69-85d4e67766b3",
"metadata": {},
"source": [
"\n",
"\n",
"#### Visualize os resultados\n",
"\n",
"O primeiro gráfico mostra que, nesta simulação, já estamos próximos `1 mH` da resposta exata após a primeira iteração (normalmente, aceita-se que a precisão química seja de `1 kcal/mol`$\\approx$`1.6 mH`). Trata-se, porém, de um sistema pequeno e, como as amostras não apresentam ruído, não é necessário recuperar a configuração. Em um sistema de maior porte executado em uma QPU ruidosa, podem ser necessárias várias iterações de recuperação da configuração, e a precisão final pode ser inferior. Geralmente, é possível melhorar a energia permitindo mais iterações de recuperação da configuração ou aumentando o número de amostras por lote.\n",
"\n",
"O segundo gráfico mostra a ocupação média de cada orbital espacial após a iteração final. Podemos ver que tanto os elétrons de spin para cima quanto os de spin para baixo ocupam os primeiros cinco orbitais com alta probabilidade em nossas soluções.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "caffd888-e89c-4aa9-8bae-4d1bb723b35e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Data for energies plot\n",
"x1 = range(len(result_history))\n",
"min_e = [\n",
" min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy\n",
" for result in result_history\n",
"]\n",
"e_diff = [abs(e - reference_energy) for e in min_e]\n",
"yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]\n",
"\n",
"# Chemical accuracy (+/- 1 milli-Hartree)\n",
"chem_accuracy = 0.001\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, e_diff, label=\"energy error\", marker=\"o\")\n",
"axs[0].set_xticks(x1)\n",
"axs[0].set_xticklabels(x1)\n",
"axs[0].set_yticks(yt1)\n",
"axs[0].set_yticklabels(yt1)\n",
"axs[0].set_yscale(\"log\")\n",
"axs[0].set_ylim(1e-4)\n",
"axs[0].axhline(\n",
" y=chem_accuracy,\n",
" color=\"#BF5700\",\n",
" linestyle=\"--\",\n",
" label=\"chemical accuracy\",\n",
")\n",
"axs[0].set_title(\"Approximated Ground State Energy Error vs SQD Iterations\")\n",
"axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n",
"axs[0].set_ylabel(\"Energy Error (Ha)\", 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",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "ce0eecb3-8a23-4118-aa1e-a28afcec6334",
"metadata": {},
"source": [
"\n",
"\n",
"## Exemplo de hardware em grande escala\n",
"\n",
"Agora, vamos executar um exemplo maior em hardware quântico real. Aqui, derivaremos um espaço ativo para a molécula de nitrogênio a partir do conjunto de bases cc-pVDZ.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "24ca3090-3f3b-4efb-a482-70b2e1b5d062",
"metadata": {},
"source": [
"\n",
"\n",
"### Etapas 1 a 4\n",
"\n",
"Aqui, reunimos todas as etapas em um único fluxo de trabalho em maior escala, que é então executado em hardware quântico real.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3858949c-a55d-4ff8-a0fc-54fb53e131b5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"converged SCF energy = -108.929838385609\n",
"norb = 26\n",
"nelec = (5, 5)\n",
"E(CCSD) = -109.2177884185544 E_corr = -0.2879500329450045\n",
"Using backend ibm_boston\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].\n",
"Removing interaction (24, 24) from the end.\n",
"Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].\n",
"Removing interaction (20, 20) from the end.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Gate counts: OrderedDict({'sx': 7039, 'rz': 6990, 'cz': 1858, 'x': 61, 'measure': 52, 'barrier': 1})\n",
"Fraction of sampled configurations that are valid: 0.02124\n",
"Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07\n",
"Iteration 1\n",
"\tSubsample 0\n",
"\t\tEnergy: -109.13889134249762\n",
"\t\tSubspace dimension: 120409\n",
"\tSubsample 1\n",
"\t\tEnergy: -109.11785470455858\n",
"\t\tSubspace dimension: 110889\n",
"\tSubsample 2\n",
"\t\tEnergy: -109.13234360554011\n",
"\t\tSubspace dimension: 130321\n",
"Iteration 2\n",
"\tSubsample 0\n",
"\t\tEnergy: -109.16392179579177\n",
"\t\tSubspace dimension: 223729\n",
"\tSubsample 1\n",
"\t\tEnergy: -109.16281938332986\n",
"\t\tSubspace dimension: 223729\n",
"\tSubsample 2\n",
"\t\tEnergy: -109.16955816711932\n",
"\t\tSubspace dimension: 233289\n",
"Iteration 3\n",
"\tSubsample 0\n",
"\t\tEnergy: -109.17905772999075\n",
"\t\tSubspace dimension: 324900\n",
"\tSubsample 1\n",
"\t\tEnergy: -109.17532445048462\n",
"\t\tSubspace dimension: 357604\n",
"\tSubsample 2\n",
"\t\tEnergy: -109.1733168689756\n",
"\t\tSubspace dimension: 348100\n",
"Iteration 4\n",
"\tSubsample 0\n",
"\t\tEnergy: -109.18437778820451\n",
"\t\tSubspace dimension: 474721\n",
"\tSubsample 1\n",
"\t\tEnergy: -109.18450164209159\n",
"\t\tSubspace dimension: 476100\n",
"\tSubsample 2\n",
"\t\tEnergy: -109.18493571190754\n",
"\t\tSubspace dimension: 487204\n",
"Iteration 5\n",
"\tSubsample 0\n",
"\t\tEnergy: -109.18616522497996\n",
"\t\tSubspace dimension: 622521\n",
"\tSubsample 1\n",
"\t\tEnergy: -109.18652868888333\n",
"\t\tSubspace dimension: 644809\n",
"\tSubsample 2\n",
"\t\tEnergy: -109.18753326484406\n",
"\t\tSubspace dimension: 585225\n",
"Final energy: -109.18753326484406\n",
"Final energy error: 0.040495951813099396\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# ------------------------------ Step 1 ------------------------------\n",
"# Build N2 molecule\n",
"mol = pyscf.gto.Mole()\n",
"mol.build(\n",
" atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]],\n",
" basis=\"cc-pvdz\",\n",
" symmetry=\"Dooh\",\n",
")\n",
"\n",
"# Define active space\n",
"n_frozen = 2\n",
"active_space = range(n_frozen, mol.nao_nr())\n",
"\n",
"# Get molecular integrals\n",
"scf = pyscf.scf.RHF(mol).run()\n",
"norb = len(active_space)\n",
"n_electrons = int(sum(scf.mo_occ[active_space]))\n",
"n_alpha = (n_electrons + mol.spin) // 2\n",
"n_beta = (n_electrons - mol.spin) // 2\n",
"nelec = (n_alpha, n_beta)\n",
"cas = pyscf.mcscf.CASCI(scf, norb, nelec)\n",
"mo = cas.sort_mo(active_space, base=0)\n",
"hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n",
"eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)\n",
"\n",
"# Store reference energy from SCI calculation performed separately\n",
"reference_energy = -109.22802921665716\n",
"\n",
"print(f\"norb = {norb}\")\n",
"print(f\"nelec = {nelec}\")\n",
"\n",
"# Get CCSD t2 amplitudes for initializing the ansatz\n",
"ccsd = pyscf.cc.CCSD(\n",
" scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]\n",
").run()\n",
"t1 = ccsd.t1\n",
"t2 = ccsd.t2\n",
"\n",
"# Set ansatz properties\n",
"n_reps = 1\n",
"pairs_aa = [(p, p + 1) for p in range(norb - 1)]\n",
"\n",
"# Let generate_lucj_pass_manager determine the alpha-beta interactions\n",
"pairs_ab = None\n",
"\n",
"# Initialize backend\n",
"service = QiskitRuntimeService()\n",
"backend = service.least_busy(\n",
" operational=True, simulator=False, min_num_qubits=133\n",
")\n",
"print(f\"Using backend {backend.name}\")\n",
"\n",
"# Create pass manager\n",
"pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(\n",
" backend=backend,\n",
" norb=norb,\n",
" connectivity=\"heavy-hex\",\n",
" interaction_pairs=(pairs_aa, pairs_ab),\n",
" optimization_level=3,\n",
")\n",
"\n",
"# Create the LUCJ ansatz operator\n",
"ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
" t2=t2,\n",
" t1=t1,\n",
" n_reps=n_reps,\n",
" interaction_pairs=(pairs_aa, pairs_ab),\n",
" # Setting optimize=True enables the \"compressed\" factorization\n",
" optimize=True,\n",
" # Limit the number of optimization iterations to prevent the code cell\n",
" # from running too long. Removing this line may improve results.\n",
" options=dict(maxiter=1000),\n",
")\n",
"\n",
"# create an empty quantum circuit\n",
"qubits = QuantumRegister(2 * norb, name=\"q\")\n",
"circuit = QuantumCircuit(qubits)\n",
"\n",
"# prepare Hartree-Fock state as the reference state and append it\n",
"# to the quantum circuit\n",
"circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)\n",
"\n",
"# apply the UCJ operator to the reference state\n",
"circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n",
"circuit.measure_all()\n",
"\n",
"\n",
"# ------------------------------ Step 2 ------------------------------\n",
"\n",
"isa_circuit = pass_manager.run(circuit)\n",
"print(f\"Gate counts: {isa_circuit.count_ops()}\")\n",
"\n",
"\n",
"# ------------------------------ Step 3 ------------------------------\n",
"sampler = Sampler(mode=backend)\n",
"sampler.options.environment.job_tags = [\"TUT_SQD\"]\n",
"job = sampler.run([isa_circuit], shots=100_000)\n",
"primitive_result = job.result()\n",
"pub_result = primitive_result[0]\n",
"\n",
"\n",
"# ------------------------------ Step 4 ------------------------------\n",
"\n",
"bit_array = pub_result.data.meas\n",
"num_valid = sum(\n",
" is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n",
")\n",
"valid_fraction = num_valid / bit_array.num_shots\n",
"print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")\n",
"expected_fraction_random = (\n",
" math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n",
")\n",
"print(\n",
" f\"Expected fraction of valid configurations from uniformly random bitstrings: \"\n",
" f\"{expected_fraction_random}\"\n",
")\n",
"# SQD options\n",
"energy_tol = 1e-3\n",
"occupancies_tol = 1e-3\n",
"max_iterations = 5\n",
"\n",
"# Eigenstate solver options\n",
"num_batches = 3\n",
"samples_per_batch = 300\n",
"symmetrize_spin = True\n",
"carryover_threshold = 1e-4\n",
"max_cycle = 200\n",
"\n",
"# Use the Hartree-Fock configuration as an initial guess for the\n",
"# orbital occupancies\n",
"initial_occupancies = (\n",
" np.array([1] * n_alpha + [0] * (norb - n_alpha)),\n",
" np.array([1] * n_beta + [0] * (norb - n_beta)),\n",
")\n",
"\n",
"# Pass options to the built-in eigensolver. If you just want to use the defaults,\n",
"# you can omit this step, in which case you would not specify the\n",
"# sci_solver argument in the call to diagonalize_fermionic_hamiltonian below.\n",
"sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=max_cycle)\n",
"\n",
"# List to capture intermediate results\n",
"result_history = []\n",
"\n",
"\n",
"result = diagonalize_fermionic_hamiltonian(\n",
" hcore,\n",
" eri,\n",
" bit_array,\n",
" samples_per_batch=samples_per_batch,\n",
" norb=norb,\n",
" nelec=nelec,\n",
" num_batches=num_batches,\n",
" energy_tol=energy_tol,\n",
" occupancies_tol=occupancies_tol,\n",
" max_iterations=max_iterations,\n",
" sci_solver=sci_solver,\n",
" symmetrize_spin=symmetrize_spin,\n",
" initial_occupancies=initial_occupancies,\n",
" carryover_threshold=carryover_threshold,\n",
" callback=callback,\n",
" seed=rng,\n",
")\n",
"\n",
"final_energy = result.energy + nuclear_repulsion_energy\n",
"energy_error = final_energy - reference_energy\n",
"print(f\"Final energy: {final_energy}\")\n",
"print(f\"Final energy error: {energy_error}\")\n",
"\n",
"# Data for energies plot\n",
"x1 = range(len(result_history))\n",
"min_e = [\n",
" min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy\n",
" for result in result_history\n",
"]\n",
"e_diff = [abs(e - reference_energy) for e in min_e]\n",
"yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]\n",
"\n",
"# Chemical accuracy (+/- 1 milli-Hartree)\n",
"chem_accuracy = 0.001\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, e_diff, label=\"energy error\", marker=\"o\")\n",
"axs[0].set_xticks(x1)\n",
"axs[0].set_xticklabels(x1)\n",
"axs[0].set_yticks(yt1)\n",
"axs[0].set_yticklabels(yt1)\n",
"axs[0].set_yscale(\"log\")\n",
"axs[0].set_ylim(1e-4)\n",
"axs[0].axhline(\n",
" y=chem_accuracy,\n",
" color=\"#BF5700\",\n",
" linestyle=\"--\",\n",
" label=\"chemical accuracy\",\n",
")\n",
"axs[0].set_title(\"Approximated Ground State Energy Error vs SQD Iterations\")\n",
"axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n",
"axs[0].set_ylabel(\"Energy Error (Ha)\", 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",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "405e89ea-57da-4021-bb18-91e8d583d310",
"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 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",
" * [Amplie os fluxos de trabalho de química SQD com o solucionador Dice](/docs/addons/qiskit-addon-sqd/guides/integrate-dice-solver) — uma página que mostra como usar o software Dice, mais eficiente, para a diagonalização\n",
" * [Documentação da API do complemento SQD](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) – referência para a `diagonalize_fermionic_hamiltonian` função\n",
" * [*Química além da escala da diagonalização exata em um supercomputador quantístico*](https://www.science.org/doi/10.1126/sciadv.adu9991) — 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": 60
},
"nbformat": 4,
"nbformat_minor": 5
}