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.)

Procurando a versão em C++?

Este tutorial usa Python. Para a implementação em C++, incluindo o código-fonte e as instruções de build, consulte o tutorial SqDRIFT em C++.

Resultados de aprendizagem

Aprenda a criar circuitos de menor profundidade em comparação com a Trotterização
Percorra um fluxo de trabalho completo para estimativa do estado fundamental usando qDRIFT e SQD
Aprenda a usar qiskit-fermions em conjunto com outros addons do Qiskit para implementar tal 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 de Diagonalização Quântica de Krylov Baseada em Amostras (SKQD)

Contexto

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

Este tutorial faz uso do Qiskit Fermions para criar os circuitos fermiônicos mais naturais para o algoritmo qDRIFT, seguido pelo uso de passes de layout e síntese fermiônicos antes de conectar os circuitos ao 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 pré-fator 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 é distinta 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 representa o $k$ -ésimo 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, 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, então cada uma das $n$ etapas evolui pela mesma duração $\lambda t / n$ , independentemente de qual termo foi sorteado. A uniformidade do ângulo de etapa é a característica distintiva do qDRIFT: um coeficiente influencia o resultado por meio de com que frequência seu termo é sorteado, não por quanto esse termo é rotacionado. Os índices são amostrados a partir da distribuição:

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

então a série $(k_1, \ldots, k_n)$ é uma sequência aleatória de índices de termos sorteados a partir dessa distribuição. Como os $c_i$ são positivos e somam $\lambda$ , esta é uma distribuição de probabilidade normalizada, e a expectativa do canal resultante sobre os sorteios aleatórios aproxima a 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$ .

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

Este tutorial mostra como gerar um conjunto de tais circuitos aleatorizados. Depois de criarmos esses circuitos, de maneira semelhante a como criamos um subespaço de Krylov para diferentes operadores, amostramos bitstrings de múltiplos operadores desse tipo com diferentes parâmetros de tempo. Isso garante uma maior sobreposição entre os vetores de estado fundamental e as bitstrings amostradas.

Requisitos

Antes de iniciar 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 — required packages for this notebook
!pip install -q numpy pyscf qiskit qiskit-addon-sqd qiskit-aer qiskit-fermions qiskit-ibm-runtime

# 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 de simulador

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

Lendo e preparando o FCIDump

Para este tutorial, carregaremos o Hamiltoniano de estrutura eletrônica para nitrogênio (N2). Também existem outras formas de criar operadores fermiônicos. 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 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 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, ou seja, 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 suficiente 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 arquivo original na fase ou na ordenação 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 buscá-lo no 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 declaramos norb e nelec, que são o número de orbitais e o número de elétrons, respectivamente. Em seguida, declaramos h1e e h2e, que são as integrais de um e dois elétrons, respectivamente. Todos esses também serão usados posteriormente para o 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 = "assets/sqdrift/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 = "assets/sqdrift/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 a partir 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 o qiskit-fermions, que fornece passes de transpilador e portas específicos para circuitos fermiônicos. Estes serão usados posteriormente antes dos passes tradicionais do transpilador 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 é, portanto, fixa. Isso garante a reprodutibilidade dos operadores criados, pois o passe QDriftTrotterization que usaremos futuramente amostra índices aleatórios para criar os operadores qDRIFT.

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

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

Observe que a group_terms_by_electronic_structure assume termos em ordem normal.

Filtragem de termos diagonais

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

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

o deslocamento de energia constante, 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$ .

Isoladamente, nenhum desses termos move população entre configurações de números de ocupação; eles atuam apenas sobre as fases das configurações já presentes. No entanto, eles não são inertes: essas fases relativas alimentam a interferência gerada pelos termos de excitação mais adiante no circuito, portanto, filtrá-los altera a evolução que é de fato gerada e pode mudar a distribuição de amostragem. Essa é uma aproximação deliberada na etapa de geração do circuito, feita para concentrar a amostragem nos termos de excitação, e não uma etapa que deixa a distribuição amostrada intacta. Diferentemente do agrupamento por simetria acima, que mantém 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, já que 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 permanece 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 in place. Ela os identifica a partir de sua estrutura em ordem normal — o multiconjunto de modos de criação correspondendo ao multiconjunto de modos de aniquilação —, portanto, ela só é válida em um operador que já esteja em ordem 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))

Agora que agrupamos os termos no 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 anteriormente. O operador de evolução é o Hamiltoniano. Mais tarde, executamos passes de transpilador nesses circuitos para criar 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 no qiskit-fermions para realizar otimizações em nível fermiônico, seguidas de transpilar nosso circuito para o backend escolhido. Como este é um experimento de simulador, faremos isso primeiro para o AerSimulator. Cálculo do peso para cada grupo

Nesta etapa, realizamos a amostragem qDRIFT dos termos de forma estocástica, com probabilidades proporcionais aos seus coeficientes no Hamiltoniano. O passe de transpilador qDRIFT faz isso por nós. Agora podemos criar circuitos mais rasos, que podem ser executados no hardware de forma mais eficiente apesar da conectividade limitada de qubits, mesmo quando o Hamiltoniano contém acoplamentos de longo alcance e termos de ordem superior a 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

Otimizações fermiônicas e nativas de 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 passe QDriftTrotterization:

O passe QDriftTrotterization usa internamente o cálculo de pesos e a amostragem para gerar os circuitos que usaremos para amostragem
O passe RelabelModes é outro passe de otimização que pode ser usado para permutar os modos fermiônicos, a fim de otimizar a conectividade entre qubits e reduzir a profundidade de portas; leia mais na referência da API

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

F2QLayout: O gerenciador de passes predefinido aplica o passe TrivialF2QLayout, que mapeia trivialmente $n$ bits fermiônicos para $n$ qubits.
F2QSynth: Um passe de transpilação para mapear 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))

Agora que terminamos as otimizações em 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. Combinaremos todas as contagens dos diferentes circuitos. Nós as convertemos em vetores booleanos antes de, finalmente, pós-processá-las com o 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 sobre as 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, finalmente, 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, adicionamos a nuclear_repulsion_energy à energia resultante.

Observação: A dimensão do subespaço não é fixa entre iterações, mesmo no simulador sem ruído — cada subamostra extrai um conjunto diferente de configurações, e a etapa de recuperação remodela o conjunto entre iterações, de modo que a dimensão relatada varia de uma subamostra para a outra. A amostragem sem ruído não fixa, por si só, a dimensão do subespaço selecionado. A execução em hardware, no entanto, tende a produzir subespaços sistematicamente maiores, porque as medições ruidosas quebram a simetria do número de partículas, e a recuperação de configurações as transforma em vetores de base adicionais. Por isso, também introduziremos 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 de hardware

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

O custo da etapa clássica não é determinado diretamente pelo número de qubits. O SQD diagonaliza o Hamiltoniano projetado sobre o subespaço gerado pelas configurações amostradas, então o que determina o custo clássico é a dimensão desse subespaço selecionado — governado aqui por samples_per_batch, num_batches e quantas configurações distintas os circuitos realmente produzem — juntamente com a álgebra linear esparsa necessária para aplicar o Hamiltoniano projetado. O espaço CI completo cresce combinatoriamente com orbitais e 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 um tanto independente: um espaço orbital mais amplo amostrado em um subespaço modesto pode ser mais barato do que um sistema menor diagonalizado sobre um subespaço muito grande.

Na prática, então, 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 autossolver. Espaços orbitais maiores normalmente exigem um subespaço maior para atingir a precisão química, e é isso que eventualmente motiva o uso de recursos distribuídos — veja qiskit-addon-sqd-hpc para escalar essa etapa. Em vez de assumir um corte fixo, a abordagem prática é observar 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 você esgote a memória disponível.

Observação: Devido ao erro de amostragem causado pelo ruído do hardware, o subespaço criado para diagonalização na execução em hardware será maior do que o obtido ao usar o simulador. Embora isso aumente a dimensão do subespaço que queremos diagonalizar, o fluxo de trabalho ainda nos fornece uma resposta precisa devido à robustez do SQD em relação 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 prosseguir sem podar. Ignorar a poda é geralmente preferível para execuções em hardware, porque isso deixa as medições com simetria quebrada disponíveis para a recuperação de configurações, que pode repará-las em configurações válidas e, assim, ampliar o subespaço em vez de simplesmente descartar essas medições.

Como o nitrogênio só pode ter sete elétrons $\alpha$ e sete $\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, caso não sejam, as descarta. Depois de filtrar as bitstrings espúrias, o restante é enviado ao esquema de diagonalização. Use a flag PRUNE abaixo para alternar entre os dois comportamentos.

Tenha em mente que a poda é apenas uma das várias escolhas que moldam o subespaço final, junto 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 tudo o mais for mantido fixo; o companheiro em C++ discute isso com mais detalhes, já que ele pós-seleciona em vez de recuperar e também difere nesses outros parâmetros.

name = "assets/sqdrift/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 pelo seguinte material:

Diagonalização quântica de Krylov baseada em amostras de um modelo de rede fermiônico - 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 cluster Jastrow unitário local (LUCJ) para simulação de química quântica.
O artigo SqDRIFT - a literatura na qual este tutorial é baseado. (Observe que algumas das otimizações discutidas neste artigo estão atualmente em desenvolvimento, e este tutorial está sujeito a alterações futuras com base na evolução das bibliotecas usadas.)

Resultados de aprendizagem​

Pré-requisitos​

Contexto​

Requisitos​

Configuração​

Exemplo de simulador​

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

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

Etapa 3: Executar usando primitivas do Qiskit​

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

Exemplo de hardware​

Próximos passos​