Algoritmo SqDRIFT para estimativa do estado fundamental

Estimativa de uso: 180 segundos em um processador Heron r3 (OBSERVAÇÃO: esta é apenas uma estimativa. Seu tempo de execução pode variar.)

Resultados de aprendizado

• Aprenda a criar circuitos de menor profundidade em comparação com a trotterização

• Acompanhe um fluxo de trabalho completo de estimativa do estado fundamental usando qDRIFT e SQD

• Aprenda a usar qiskit-fermions em conjunto com outros addons do Qiskit para implementar esse fluxo de trabalho

Este tutorial é apresentado como um notebook Python para fins didáticos.

Pré-requisitos

• Leia a visão geral de Diagonalização quântica baseada em amostras (SQD)

• Leia a lição Diagonalização quântica de Krylov baseada em amostras (SKQD)

Contexto

O SqDRIFT é uma variante do SKQD que substitui a necessidade de escolher um ansatz do qual amostrar bitstrings por um conjunto de circuitos de evolução temporal construídos diretamente a partir do Hamiltoniano alvo. Isso é feito subamostrando operadores de evolução temporal menores do Hamiltoniano com base em seus coeficientes, o que é conhecido como método de trotterização qDRIFT.

Este tutorial usa o Qiskit Fermions para criar os circuitos fermiônicos mais naturais para o algoritmo qDRIFT, seguido do uso de passes de layout e síntese fermiônicos antes de inserir os circuitos no pipeline tradicional do Qiskit para execução em hardware.

Seja o Hamiltoniano da forma:

$$H = \sum_{i=1}^{N} c_i h_i$$

onde, sem perda de generalidade, exigimos $c_i > 0$ e que o maior autovalor de $h_i$ seja igual, em valor absoluto, a $1$. Qualquer prefator com sinal ou complexo é absorvido em $h_i$, de modo que os coeficientes $c_i$ são pesos estritamente positivos, enquanto os $h_i$ carregam 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; é uma propriedade do Hamiltoniano e é distinto do número de operadores amostrados em um único circuito, escrito como $n$ abaixo.

O algoritmo qDRIFT então realiza, para o tempo alvo $t$, algum operador $V_k$, onde $k$ vai de $1 \cdots K$ e designa o $k_{th}$ circuito SqDRIFT, definido como:

$$V_k = \prod_{j=1}^{n} e^{-i h_{k_j} \lambda t / n }$$

Aqui $n$ é o número de operadores amostrados por circuito e $K$ é o número de circuitos no conjunto. O produto percorre os $n$ sorteios, e não todos os $N$ termos do Hamiltoniano, 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$.

A quantidade:

$$\lambda = \sum_{i=1}^{N} c_i$$

é a norma $L_1$ dos coeficientes, de modo que cada um dos $n$ passos evolui pela mesma duração $\lambda t / n$, independentemente de qual termo foi sorteado. A uniformidade do ângulo do passo é a característica marcante do qDRIFT: um coeficiente influencia o resultado por com que frequência seu termo é sorteado, e não por quanto esse termo é girado. Os índices são amostrados da distribuição:

$$P[k_i] = \frac{c_i}{\lambda}$$

de modo que a série $(k_1, \ldots, k_n)$ é uma sequência aleatória de índices de termos sorteados dessa distribuição. Como os $c_i$ são positivos e somam $\lambda$, essa é uma distribuição de probabilidade normalizada, e o valor esperado do canal resultante sobre os sorteios aleatórios aproxima a evolução sob $H$, com um erro que diminui à medida que $n$ cresce. Note que o erro de aproximação depende de $\lambda$ e não do número de termos $N$.

(O artigo do SqDRIFT escreve o número de termos como $\mathcal{N}$ e o comprimento da sequência como $N$; usamos $N$ e $n$ aqui para manter os dois claramente distintos.)

Este tutorial mostra como gerar um conjunto desses circuitos aleatorizados. Depois de criarmos esses circuitos, de forma semelhante a como criamos um subespaço de Krylov para diferentes operadores, amostramos bitstrings 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 bitstrings amostradas.

Requisitos

Antes de começar este tutorial, certifique-se de ter instalado

• Um ambiente virtual Python (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (observe que o nome está no plural)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Você pode instalar todos os pacotes necessários com:

pip install "qiskit~=2.5" "qiskit-fermions==0.1.0" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy

Configuração

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "qiskit_aer": "qiskit-aer", "qiskit_fermions": "qiskit-fermions", "qiskit_ibm_runtime": "qiskit-ibm-runtime"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

# Qiskit core
from qiskit import transpile
from qiskit.primitives import BitArray

# Qiskit Aer
from qiskit_aer import AerSimulator

# IBM Quantum Compute Service
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler

# Qiskit Fermions
from qiskit_fermions.operators.library import FCIDump
from qiskit_fermions.operators import FermionOperator
from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms
from qiskit_fermions.operators.terms.grouping import (
    group_terms_by_electronic_structure,
)
from qiskit_fermions.operators.terms.ordering import canonical_order
from qiskit_fermions.circuit import FermionicCircuit
from qiskit_fermions.circuit.library import Evolution
from qiskit_fermions.transpiler import FermionicPassManager
from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
from qiskit_fermions.transpiler.passes import QDriftTrotterization
from qiskit_fermions.circuit.library import InitializeModes

# Qiskit addon SQD
from qiskit_addon_sqd.fermion import (
    diagonalize_fermionic_hamiltonian,
    SCIResult,
)

Exemplo com simulador

Etapa 1: Mapear entradas clássicas para um problema quântico

Leitura e preparação do FCIDump

Neste tutorial, vamos carregar o Hamiltoniano de estrutura eletrônica do nitrogênio (N2). Há outras formas de criar operadores fermiônicos também. Consulte a documentação em qiskit_fermions.operators.library.

Sobre este FCIDump. O arquivo N2_sto_3g descreve uma molécula de nitrogênio ($N_2$) na base mínima STO-3G, com 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 spin-orbitais e 20 qubits sob Jordan-Wigner), 14 elétrons em um singleto de spin, ou seja, sete elétrons $\alpha$ e sete $\beta$. Todos os orbitais recebem o rótulo de simetria 1, isto é, nenhuma simetria de grupo pontual é explorada. Por ser um dump STO-3G de espaço completo, nenhum orbital é congelado e o espaço de correlação é pequeno o bastante para que uma energia de referência FCI exata possa ser calculada classicamente para comparação, como mostrado na próxima célula.

Um arquivo equivalente pode ser regenerado com o PySCF:

from pyscf import gto, scf, tools

mol = gto.M(atom="N 0 0 0; N 0 0 1.09", basis="sto-3g", symmetry=False)
mf = scf.RHF(mol).run()
tools.fcidump.from_scf(mf, "N2_sto_3g")

Como as integrais dependem dos orbitais SCF convergidos, um arquivo regenerado pode diferir do fornecido na fase ou na ordem dos orbitais; as energias totais não são afetadas.

Obtendo o arquivo. Encontre o FCIDump neste repositório do GitHub. Você pode executar a célula abaixo para baixá-lo para o local que o restante do tutorial espera.

Primeiro usamos o cisolver fornecido 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, primeiro declararemos norb e nelec, que são o número de orbitais e o número de elétrons, respectivamente. Depois declaramos h1e e h2e, que são as integrais de um e de dois elétrons, respectivamente. Todos esses valores também serão usados depois no SQD.

import os
from urllib.request import urlopen

# The FCIDump is stored with this tutorial in the Qiskit documentation repository.
FCIDUMP_URL = "https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g"
FCIDUMP_PATH = "fcidump_files/N2_sto_3g"

if not os.path.exists(FCIDUMP_PATH):
    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)
    with urlopen(FCIDUMP_URL) as response:
        contents = response.read()
    with open(FCIDUMP_PATH, "wb") as f:
        f.write(contents)
    print(f"Downloaded FCIDump to {FCIDUMP_PATH}")
else:
    print(f"Using existing FCIDump at {FCIDUMP_PATH}")

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

Carregando o Hamiltoniano

Com os dados necessários prontos, lemos o Hamiltoniano do arquivo FCI em um formato compatível com qiskit-fermions

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

Fluxos de trabalho fermiônicos com qiskit-fermions

Primeiro mapearemos o Hamiltoniano em um modelo de circuito fermiônico usando qiskit-fermions, que fornece passes de transpilador e portas específicos para circuitos fermiônicos. Eles serão usados depois, antes dos passes de transpilador tradicionais do Qiskit, neste fluxo de trabalho.

Agrupamento de termos

Para garantir a reprodutibilidade dos resultados, primeiro usamos canonical_order para ordenar os termos com base apenas em sua estrutura. A ordem dos operadores na lista canon fica, portanto, fixa. Isso garante a reprodutibilidade dos operadores criados, porque o passe QDriftTrotterization que usaremos adiante sorteia índices aleatórios para criar os operadores qDRIFT.

Nesta etapa, exploramos as muitas simetrias presentes no Hamiltoniano de estrutura eletrônica agrupando termos relacionados com coeficientes idênticos. Embora isso altere a distribuição de coeficientes dos operadores da qual o protocolo qDRIFT amostra, isso não afeta suas garantias de convergência. Crucialmente, agrupar termos relacionados por simetria resulta em um cancelamento favorável de termos de Pauli e em uma profundidade de circuito global menor ao evoluir no tempo um estado sob a ação deles.

O qiskit-fermions fornece a função group_terms_by_electronic_structure, que faz esse agrupamento para nós.

Observe que group_terms_by_electronic_structure assume termos com ordenação normal.

Filtrando termos diagonais

Removemos os termos diagonais do Hamiltoniano usado para gerar os circuitos, de modo que os $n$ espaços de amostragem do qDRIFT sejam gastos em termos que movem população entre configurações. É melhor filtrar esses termos do Hamiltoniano neste ponto, antes de a porta Evolution ser construída na próxima etapa.

Os termos em questão são os que são diagonais na base de números de ocupação, isto é, os produtos de operadores de número $a^\dagger_i a_i$. Três tipos de termo se enquadram nessa descrição:

• o deslocamento constante de energia, um produto de zero operadores de número, cuja evolução temporal contribui apenas com uma fase global;

• os operadores de número individuais $n_i$, cuja evolução temporal se reduz a rotações $Z$ de um único qubit;

• os produtos de ordem superior, como $n_i n_j$.

Por si sós, nenhum deles move população entre configurações de números de ocupação; eles atuam apenas sobre as fases das configurações já presentes. Eles não são inertes, porém: essas fases relativas alimentam a interferência gerada pelos termos de excitação mais adiante no circuito, então filtrá-los altera a evolução que de fato é gerada e pode alterar a distribuição de amostragem. Trata-se de uma aproximação deliberada na etapa de geração de circuitos, feita para concentrar a amostragem nos termos de excitação, e não de uma etapa que deixa a distribuição amostrada intacta. Ao contrário do agrupamento por simetria acima, que deixa intactas as garantias de convergência do qDRIFT, esse filtro altera o operador que está sendo evoluído. Os circuitos, portanto, deixam de aproximar a evolução sob o Hamiltoniano completo, e os limites de erro do qDRIFT se aplicam ao operador filtrado, e não ao original. Isso é aceitável aqui porque os circuitos são apenas uma heurística de amostragem usada para propor configurações: nenhum termo é perdido da própria estimativa de energia, pois o filtro se aplica apenas ao Hamiltoniano usado para construir os circuitos, enquanto a diagonalização clássica posterior usa o Hamiltoniano completo, incluindo os termos diagonais. A precisão do SQD depende dessa etapa clássica, que continua variacional no subespaço amostrado, independentemente de como as configurações foram propostas.

A função filter_diagonal_terms() remove esses termos de um operador no local. Ela os identifica a partir de sua estrutura em ordenação normal — o multiconjunto de modos de criação coincidindo com o multiconjunto de modos de aniquilação — e por isso só é válida em um operador que já esteja em ordenação normal. Essa suposição não é verificada em tempo de execução.

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

5060

Agora que agrupamos os termos do Hamiltoniano, decidiremos os seguintes parâmetros para gerar o conjunto de circuitos:

• O número de circuitos a gerar: num_circuits
• O comprimento de cada circuito em termos de grupos de excitação: num_exc
• O fator para os diferentes tempos de evolução: times

Criando circuitos fermiônicos

Agora criaremos circuitos fermiônicos para cada um dos passos de tempo. Cada circuito consistirá em uma única porta de evolução, com o tempo de evolução que declaramos antes. O operador de evolução é o Hamiltoniano. Depois executamos passes de transpilador nesses circuitos para criar os circuitos qDRIFT.

Preparação do ansatz

Preparamos o estado de Hartree-Fock usando a classe InitializeModes. Para o nitrogênio, o processo consiste simplesmente em aplicar portas X aos primeiros num_elec_a qubits e depois aos num_elec_b qubits, ambos iguais a sete para o nitrogênio. Esse estado representa os sete elétrons $\alpha$ e sete $\beta$ do nitrogênio.

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []

hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

Etapa 2: Otimizar o problema para execução em hardware quântico

Agora que temos nossos circuitos, primeiro usaremos os passes disponíveis em qiskit-fermions para realizar otimizações no nível fermiônico, seguidas da transpilação do nosso circuito para o backend escolhido. Como este é um experimento com simulador, faremos isso primeiro para o AerSimulator. Cálculo de pesos para cada grupo

Nesta etapa, realizamos a amostragem qDRIFT de termos de forma estocástica, com probabilidades proporcionais aos seus coeficientes no Hamiltoniano. O pass do transpilador qDRIFT faz isso por nós. Agora podemos criar circuitos mais rasos, que podem ser executados no hardware com mais eficiência apesar da conectividade limitada dos qubits, mesmo quando o Hamiltoniano contém acoplamentos de longo alcance e termos de ordem superior à quadrática. Após o agrupamento de termos, ele amostra os operadores com base em seus pesos. Para cada operador $h_i$, o peso $W_{h_i}$ é definido da seguinte forma:

$$W_{h_i} = |c_i| / \lambda$$

Como os termos foram agrupados na Etapa 1, cada $h_i$ aqui é um grupo inteiro: $c_i$ é o coeficiente absoluto médio dos termos do grupo $i$, e cada termo do grupo é evoluído com seu coeficiente reduzido ao seu sinal.

Otimizações fermiônicas e nativas do hardware

A função generate_preset_jw_pass_manager() retorna um MultiStagePassManager que recebe um FermionicCircuit e produz um circuito final otimizado que podemos transpilar para executar em nosso hardware. Substituímos seu estágio de otimização padrão por um FermionicPassManager contendo nosso pass QDriftTrotterization:

• O pass QDriftTrotterization usa internamente o cálculo de pesos e a amostragem para gerar os circuitos que usaremos na amostragem

• O pass RelabelModes é outro pass de otimização que pode ser usado 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 da API

Os demais estágios do MultiStagePassManager são executados automaticamente e cuidam de todo o mapeamento de férmions para qubits:

• F2QLayout: O pass manager predefinido aplica o pass TrivialF2QLayout, que mapeia de forma trivial $n$ bits fermiônicos para $n$ qubits.

• F2QSynth: Um pass de transpilação que mapeia instruções de circuito baseadas em férmions para instruções baseadas em qubits.

qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

400

Agora que terminamos as otimizações no nível fermiônico, podemos transpilar os circuitos para execução no simulador.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Etapa 3: Executar usando primitivas do Qiskit

Agora que temos nossos circuitos, podemos executá-los usando primitivas do Qiskit no AerSimulator. Vamos combinar todas as contagens dos diferentes circuitos. Convertemos essas contagens em vetores booleanos antes de, por fim, fazer o pós-processamento com SQD.

print(
    f"Executing {len(transpiled_circuits)} circuits with {shots} shots each..."
)

job = simulator.run(transpiled_circuits, shots=shots)
result = job.result()

all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]

print(len(all_counts), "length before post processing")

Executing 400 circuits with 100 shots each...
400 length before post processing

Etapa 4: Pós-processar e retornar o resultado no formato clássico desejado

Usando bitstrings para o SQD

Agora podemos executar o esquema de diagonalização nas bitstrings selecionadas para encontrar o menor autovalor, que corresponderá à energia do estado fundamental da molécula. Criamos uma função de callback, declaramos as ocupações iniciais e definimos os parâmetros antes de, por fim, executar o esquema de diagonalização. A função de callback é usada para imprimir a iteração atual e a estimativa atual do autovalor a cada iteração.

Por fim, para obter a estimativa do estado fundamental, somamos a nuclear_repulsion_energy à energia resultante.

Nota: A dimensão do subespaço não é fixa entre as iterações, mesmo no simulador sem ruído — cada subamostra sorteia um conjunto diferente de configurações, e a etapa de recuperação remodela o conjunto entre as iterações, de modo que a dimensão reportada varia de uma subamostra para outra. A amostragem sem ruído, por si só, não fixa a dimensão do subespaço selecionado. A execução no hardware, porém, tende a produzir subespaços sistematicamente maiores, porque os shots com ruído quebram a simetria do número de partículas e a recuperação de configurações os transforma em vetores de base adicionais. Por isso, também vamos introduzir outra etapa para podar bitstrings na seção de hardware.

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

40000
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64767025226178
		Subspace dimension: 5538
	Subsample 1
		Energy: -107.64772799119115
		Subspace dimension: 5670
	Subsample 2
		Energy: -107.64765512281548
		Subspace dimension: 5767
Iteration 2
	Subsample 0
		Energy: -107.64795948524682
		Subspace dimension: 6080
	Subsample 1
		Energy: -107.64806617355072
		Subspace dimension: 6300
	Subsample 2
		Energy: -107.64802260640258
		Subspace dimension: 6308
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712
 0.99301927 0.0373282  0.0373266  0.00944508]
Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737
 0.99302145 0.03733536 0.03733399 0.0094585 ]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6480661736 Ha
Error:            1.1811817564e-04 Ha

Exemplo em hardware

Este exemplo usa 20 qubits (10 orbitais espaciais). Essa escolha é uma conveniência para um tutorial que deve ser executado rapidamente, e não um limite rígido do método.

O custo da etapa clássica não é definido diretamente pelo número de qubits. O SQD diagonaliza o Hamiltoniano projetado no subespaço gerado pelas configurações amostradas, de modo que o que determina o custo clássico é a dimensão desse subespaço selecionado — governada aqui por samples_per_batch, num_batches e pelo número de configurações distintas que os circuitos realmente produzem — junto com a álgebra linear esparsa necessária para aplicar o Hamiltoniano projetado. O espaço CI completo cresce de forma combinatória com os orbitais e os elétrons, mas o subespaço selecionado é uma fatia pequena e ajustável dele, e controlamos seu tamanho diretamente. Consequentemente, o número de qubits e a dificuldade clássica podem variar de forma relativamente independente: um espaço orbital mais amplo amostrado em um subespaço modesto pode ser mais barato do que um sistema menor diagonalizado em um subespaço muito grande.

Na prática, portanto, o tamanho de sistema viável depende da dimensão do subespaço de que você precisa para a precisão desejada e da memória e dos núcleos disponíveis para o solucionador de autovalores. Espaços orbitais maiores normalmente exigem um subespaço maior para atingir a precisão química, e é isso que acaba motivando o uso de recursos distribuídos — veja qiskit-addon-sqd-hpc para escalar essa etapa. Em vez de assumir um limite fixo, a abordagem prática é observar a dimensão do subespaço reportada e a convergência da energia ao longo das iterações, e aumentar o tamanho do subespaço até que a energia pare de melhorar ou você esgote a memória disponível.

Nota: Devido ao erro de amostragem causado pelo ruído do hardware, o subespaço criado para a diagonalização na execução em hardware será maior do que o obtido com o simulador. Embora isso aumente a dimensão do subespaço que queremos diagonalizar, o fluxo de trabalho ainda nos dá uma resposta precisa graças à robustez do SQD frente ao ruído.

Poda de strings espúrias

Aqui podemos optar por realizar uma etapa adicional. Quando temos todas as bitstrings das execuções dos circuitos, podemos filtrar as bitstrings inválidas antes de executar o SQD ou seguir em frente sem podar. Pular a poda geralmente é preferível em execuções em hardware, porque deixa os shots com simetria quebrada disponíveis para a recuperação de configurações, que pode reparar esses shots transformando-os em configurações válidas e, assim, ampliar o subespaço em vez de descartá-los de vez.

Como o nitrogênio só pode ter sete elétrons $\alpha$ e sete elétrons $\beta$, qualquer bitstring que tenha mais ou menos de sete 1s na primeira e na segunda metade da saída pode ser descartada. Definimos uma função que verifica se as bitstrings são válidas e, se não forem, as descarta. Depois de filtrar as bitstrings espúrias, as restantes são enviadas ao esquema de diagonalização. Use a flag PRUNE abaixo para alternar entre os dois comportamentos.

Lembre-se de que a poda é apenas uma de várias escolhas que moldam o subespaço final, ao lado do número de circuitos, do conjunto de tempos de evolução e da filtragem de termos diagonais. Comparar uma execução com poda a uma sem poda só é informativo se todo o resto for mantido fixo; a versão em C++ deste tutorial discute isso com mais detalhes, já que ela faz pós-seleção em vez de recuperação e também difere nesses outros parâmetros.

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []
hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

# Calculate weights for sampling (one per group)
qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService(channel="ibm_quantum_platform")

# Select backend (choose based on qubit requirements)
backend = service.least_busy(
    operational=True,
    simulator=False,
    min_num_qubits=2 * norb,
)

print(f"Selected backend: {backend.name} ({backend.num_qubits} qubits)")

# Transpile for hardware
transpiled_circuits = transpile(
    sqdrift_circuits,
    backend=backend,
    optimization_level=3,
    seed_transpiler=42,
)

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-SqDRIFT"]

job = sampler.run(transpiled_circuits, shots=shots)
result = job.result()

# Extract counts from SamplerV2 results
all_counts = [pub_result.data.meas.get_counts() for pub_result in result]

# Set to True to filter out bitstrings that violate electron-number conservation
PRUNE = False

def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )

if PRUNE:
    all_counts_filtered = []
    for counts in all_counts:
        filtered_count = {}
        for key in counts:
            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):
                continue
            elif key not in filtered_count.keys():
                filtered_count[key] = counts[key]
            else:
                filtered_count[key] += counts[key]
        all_counts_filtered.append(filtered_count)
    all_counts = all_counts_filtered

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print("Electron configuration:")
print(f"  Total electrons: {nelec}")
print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha
5060
400

Selected backend: ibm_aachen (156 qubits)
40000
Electron configuration:
  Total electrons: 14
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64593072647523
		Subspace dimension: 7221
	Subsample 1
		Energy: -107.6458270048177
		Subspace dimension: 7209
	Subsample 2
		Energy: -107.64007673117075
		Subspace dimension: 7138
Iteration 2
	Subsample 0
		Energy: -107.64757372124944
		Subspace dimension: 9009
	Subsample 1
		Energy: -107.64674060104392
		Subspace dimension: 8245
	Subsample 2
		Energy: -107.64731360491942
		Subspace dimension: 8178
Iteration 3
	Subsample 0
		Energy: -107.64765518770588
		Subspace dimension: 8835
	Subsample 1
		Energy: -107.64767975712016
		Subspace dimension: 8649
	Subsample 2
		Energy: -107.64761634415606
		Subspace dimension: 8648
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295
 0.99298797 0.03728154 0.03728186 0.00938359]
Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174
 0.99298424 0.03728121 0.03728169 0.00939159]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6476797571 Ha
Error:            5.0453460619e-04 Ha

Próximos passos

Recomendações

Se você achou este trabalho interessante, talvez se interesse pelos seguintes materiais:

• Diagonalização quântica de Krylov baseada em amostras de um modelo de rede fermiônica - um tutorial relacionado que usa circuitos de evolução temporal em vez de um ansatz variacional.
• Diagonalização quântica baseada em amostras de um Hamiltoniano de química - um tutorial sobre como construir um circuito de Jastrow de cluster unitário local (LUCJ) para simulação de química quântica.
• O artigo SqDRIFT - a literatura na qual este tutorial se baseia. (Observe que algumas das otimizações discutidas nesse artigo ainda estão em desenvolvimento, e este tutorial está sujeito a mudanças no futuro, conforme a evolução das bibliotecas utilizadas.)