Algoritmo SqDRIFT para estimativa do estado fundamental

Python

Baixe o notebook Python

Baixe o notebook Python para executar este tutorial por conta própria. O notebook baixa a entrada FCIDump caso ela ainda não esteja presente.

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

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.

Ambas as implementações usam o Qiskit Fermions para representar o Hamiltoniano molecular e construir circuitos para o algoritmo qDRIFT. A implementação em Python usa passes fermiônicos do transpilador, enquanto a implementação em C++ amostra grupos explicitamente e os mapeia para operadores de qubit.

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

# 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,
)

Implementação

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

Os passos em Python abaixo usam o simulador Aer. O exemplo em hardware vem em seguida.

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 no próximo bloco de código.

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 o código 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 = "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}")

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

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

5060

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$$

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

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

400

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 eigensolver (resolvedor de autovalores). 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; a aba C++ discute isso com mais detalhes, já que ela pós-seleciona em vez de recuperar 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")

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

C++

Obtenha o código C++

O código-fonte em C++, o projeto CMake e as instruções de compilação estão no diretório docs/tutorials/assets/sqdrift/cpp do repositório da documentação do Qiskit. A entrada FCIDump está no diretório vizinho fcidump_files. Clone o repositório, ou baixe os dois diretórios, para compilar e executar este tutorial por conta própria.

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 aprendizagem (C++)

• 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

Pré-requisitos (C++)

• 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 (C++)

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.

Ambas as implementações usam o Qiskit Fermions para representar o Hamiltoniano molecular e construir circuitos para o algoritmo qDRIFT. A implementação em Python usa passes fermiônicos do transpilador, enquanto a implementação em C++ amostra grupos explicitamente e os mapeia para operadores de qubit.

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 (C++)

Esta implementação é executada em hardware IBM Quantum® e usa um executável SBD externo para a diagonalização clássica. Ela não inclui o exemplo com simulador em Python.

Siga as instruções de compilação em C++ para instalar o compilador, o CMake, o Rust, o Boost e as dependências do Qiskit. Instale o Rust com o rustup: as dependências do Qiskit fixam sua própria versão do Rust, que o rustup baixa automaticamente, mas um Rust instalado pelo Homebrew ou por um gerenciador de pacotes do Linux pode ser antigo demais. A etapa final do SBD requer, além disso, MPI, OpenMP e BLAS/LAPACK. O Windows não é suportado atualmente por esta implementação.

Use o código-fonte em C++ completo e o projeto CMake do repositório. Mantenha o diretório adjacente fcidump_files ao baixar o exemplo. A partir de um checkout do repositório da documentação, compile com os comandos a seguir. A primeira compilação baixa e compila todas as dependências e leva vários minutos.

cd docs/tutorials/assets/sqdrift/cpp
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --config Release

A compilação grava o executável SqDRIFT em docs/tutorials/assets/sqdrift/cpp (a raiz do projeto), não em build/. Execute ./SqDRIFT a partir da raiz do projeto para que ele encontre ../fcidump_files/N2_sto_3g. Todos os outros comandos desta aba também são executados a partir da raiz do projeto, a menos que digam o contrário. Se ./SqDRIFT falhar com um erro de biblioteca compartilhada ausente, siga as instruções de caminho de bibliotecas específicas da plataforma no guia de compilação; no macOS, o executável normalmente roda sem elas.

Configuração (C++)

Os trechos abaixo são excertos do programa SqDRIFT.cpp completo. Compile e execute esse arquivo como um único programa.

Antes de executar o SqDRIFT, certifique-se de que seu token de API do IBM Quantum esteja configurado para a IBM Cloud Platform.

Se você tem Python com qiskit-ibm-runtime instalado, a maneira mais fácil de criar um arquivo de configuração corretamente formatado é:

from qiskit_ibm_runtime import QiskitRuntimeService
QiskitRuntimeService.save_account(channel="ibm_quantum_platform", token="YOUR_IBM_QUANTUM_TOKEN_HERE")

Isso grava ~/.qiskit/qiskit-ibm.json (Linux/macOS) ou %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) com o formato correto. Se preferir criar o arquivo manualmente:

{
  "default-ibm-quantum-platform": {
    "channel": "ibm_quantum_platform",
    "token": "YOUR_IBM_QUANTUM_TOKEN_HERE"
  }
}

Obtenha sua chave de API no IBM Quantum Platform.

Observação: Um arquivo de configuração ausente ou malformado causa um panic unwrap() do Rust pouco útil na inicialização, em vez de um erro claro. Verifique novamente o nome da chave (default-ibm-quantum-platform) e o valor do canal (ibm_quantum_platform) se você tiver um problema.

Implementação (C++)

Etapa 1: Mapeie entradas clássicas para um problema quântico (C++)

Carregar o Hamiltoniano molecular

char filename[] = "../fcidump_files/N2_sto_3g";
QfFCIDump* fcidump = qf_fcidump_from_file(filename);
if (fcidump == nullptr) {
    std::cerr << "Failed to load FCIDump!" << std::endl;
    return 1;
}

uint32_t norb = qf_fcidump_norb(fcidump);
uint32_t nelec = qf_fcidump_nelec(fcidump);
uint32_t n_alpha = nelec / 2;
uint32_t n_beta  = nelec - n_alpha;
uint32_t num_modes = 2 * norb;
uint32_t num_qubits = num_modes;
std::cout << "✓ Loaded N2 molecule (" << norb << " orbitals, "
          << num_qubits << " qubits, " << nelec << " electrons)" << std::endl;

Output:

✓ Loaded N2 molecule (10 orbitals, 20 qubits, 14 electrons)

Carrega o Hamiltoniano da molécula de $N_2$ a partir de um arquivo FCIDump. O número de qubits é igual ao dobro do número de orbitais (spin para cima e spin para baixo).

Criar o Hamiltoniano e colocá-lo em ordem normal

// Create and normal order Hamiltonian
QfFermionOperator* hamiltonian = qf_ferm_op_from_fcidump(fcidump);
QfFermionOperator* normal_ordered = qf_ferm_op_normal_ordered(hamiltonian, nullptr);

Converte os dados do FCIDump em um operador fermiônico e aplica a ordenação normal para simplificar a estrutura do operador.

Agrupar termos por estrutura eletrônica

// Group terms by electronic structure
QfExitCode exit_code = qf_ferm_op_group_terms_by_electronic_structure(
    normal_ordered, num_modes, false
);
if (exit_code != QfExitCode_Success) {
    std::cerr << "Failed to group terms!" << std::endl;
    return 1;
}

uint32_t num_groups = qf_ferm_op_num_groups(normal_ordered);
std::cout << "Grouped into " << num_groups << " groups" << std::endl;

// Split into group operators
QfFermionOperator** group_ops = new QfFermionOperator*[num_groups];
qf_ferm_op_split_out_groups(normal_ordered, nullptr, 0, group_ops);

Output:

Grouped into 1902 groups

Agrupa termos relacionados do Hamiltoniano por estrutura eletrônica antes de construir seus circuitos de evolução.

Calcular os pesos de amostragem, depois normalizar e mapear para operadores de qubit

A ordem importa aqui. A normalização redimensiona todos os coeficientes para magnitude unitária no próprio lugar, portanto os pesos de amostragem devem ser calculados primeiro. Se a normalização for executada primeiro, todo $|c_j|$ volta como exatamente 1.0 e o peso de cada grupo se reduz a uma simples contagem de seus termos — a distribuição de amostragem seria então determinada pelo tamanho do grupo, e não pelos coeficientes físicos.

O peso de cada grupo é o coeficiente absoluto médio de seus termos. O passe QDriftTrotterization da aba Python usa o mesmo peso para um Hamiltoniano agrupado. Usar a soma em vez disso faria $\lambda$ escalar com o tamanho do grupo, o que, por sua vez, distorceria o tempo de evolução $t_k = \lambda\tau/n$.

// Calculate sampling weights BEFORE normalization.
std::vector<double> weights(num_groups);
double total_weight = 0.0;

for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    double abs_sum = 0.0;
    for (uint64_t j = 0; j < num_terms; j++) {
        abs_sum += std::sqrt(coeffs[j].re * coeffs[j].re +
                             coeffs[j].im * coeffs[j].im);
    }
    weights[i] = (num_terms > 0) ? (abs_sum / static_cast<double>(num_terms)) : 0.0;
    total_weight += weights[i];
}

std::cout << "Total weight (λ): " << total_weight << std::endl;

// Normalize the evolution operators term-by-term, now that the sampling
// weights above have already captured the original magnitudes.
for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    for (uint64_t j = 0; j < num_terms; j++) {
        const double magnitude = std::sqrt(coeffs[j].re * coeffs[j].re +
                                           coeffs[j].im * coeffs[j].im);
        if (magnitude > 0.0) {
            coeffs[j].re /= magnitude;
            coeffs[j].im /= magnitude;
        }
    }
}

QkObs** qubit_ops = new QkObs*[num_groups];
for (uint32_t i = 0; i < num_groups; i++) {
    QfExitCode jw_exit = qf_ferm_op_jordan_wigner(group_ops[i], num_qubits, &qubit_ops[i]);
    if (jw_exit != QfExitCode_Success) {
        std::cerr << "Failed to map group " << i << std::endl;
        return 1;
    }
}
std::cout << "Mapped all " << num_groups << " normalized groups to qubit operators" << std::endl;

Output:

Total weight (λ): 337.261
Mapped all 1902 normalized groups to qubit operators

O peso $\lambda_i = \frac{1}{|G_i|}\sum_{j \in G_i} |c_{ij}|$ é o coeficiente absoluto médio sobre os termos do grupo $G_i$, lido a partir dos coeficientes originais. A normalização então redimensiona cada termo fermiônico para magnitude unitária antes da transformação de Jordan-Wigner, preservando apenas sua fase/sinal no operador evoluído. Isso corresponde à evolução de termos normalizados do qDRIFT: as magnitudes dos coeficientes determinam as probabilidades de amostragem por meio de $\lambda$, enquanto a evolução do circuito usa o termo normalizado, de modo que coeficientes grandes não sejam contados duas vezes.

Amostrar conjuntos de operadores

// SqDRIFT Sampling: Create operator sets
    const int num_circuits = 100;    // Number of circuits to create
    const int ops_per_circuit = 10;  // Operators per circuit
    const double time_step = 1;    // Time step for evolution
    
    std::cout << "\n Generating " << num_circuits << " circuits with " 
              << ops_per_circuit << " operators each..." << std::endl;
    
    std::mt19937 gen(42);
    std::discrete_distribution<> dist(weights.begin(), weights.end());
    
    // Store sampled operator indices
    std::vector<std::vector<int>> operator_sets(num_circuits);
    for (int i = 0; i < num_circuits; i++) {
        operator_sets[i].resize(ops_per_circuit);
        for (int j = 0; j < ops_per_circuit; j++) {
            operator_sets[i][j] = dist(gen);
        }
    }

Output:

 Generating 100 circuits with 10 operators each...

Sorteia os conjuntos de operadores para 100 circuitos aleatórios, cada um contendo 10 grupos amostrados da distribuição ponderada. A próxima etapa transforma esses conjuntos no conjunto de circuitos qDRIFT.

Etapa 2: Construir os circuitos qDRIFT (C++)

Na versão em C++, esta etapa constrói os circuitos a partir dos conjuntos de operadores amostrados. A transpilação para o backend de destino acontece na Etapa 3, imediatamente antes de cada circuito ser enviado.

    // Create Suzuki-Trotter circuits by composing all operators
    std::vector<Qiskit::circuit::QuantumCircuit> circuits;
    circuits.reserve(num_circuits);

    const uint32_t trotter_order = 1;
    const uint32_t trotter_reps = 1;
    const bool preserve_order = false;
    const bool insert_barriers = false;

    // Gate name mapping is static — compute once before all loops
    auto name_map = Qiskit::circuit::get_standard_gate_name_mapping();

    for (int circuit_idx = 0; circuit_idx < num_circuits; circuit_idx++) {
        // Create initial circuit with state preparation using C++ API
        Qiskit::circuit::QuantumCircuit qc(num_qubits, num_qubits);

        // State preparation: apply X gates to spin-up qubits 0..n_alpha-1
        // and spin-down qubits norb..norb+n_beta-1 (derived from FCIDump NELEC)
        for (uint32_t q = 0; q < n_alpha && q < num_qubits; q++) {
            qc.x(q);
        }
        for (uint32_t q = norb; q < norb + n_beta && q < num_qubits; q++) {
            qc.x(q);
        }

        // For each operator in this circuit's set
        for (int op_idx = 0; op_idx < ops_per_circuit; op_idx++) {
            int group_idx = operator_sets[circuit_idx][op_idx];
            double evolution_time = (total_weight * time_step) / ops_per_circuit;
            
            // Create Suzuki-Trotter evolution circuit
            QkCircuit* evolution_circuit = qk_circuit_library_suzuki_trotter(
                qubit_ops[group_idx], trotter_order, trotter_reps,
                evolution_time, preserve_order, insert_barriers
            );
            
            if (evolution_circuit != nullptr) {
                // Get number of instructions
                uint32_t num_ops = qk_circuit_num_instructions(evolution_circuit);
                
                // Manually append each instruction
                for (uint32_t i = 0; i < num_ops; i++) {
                    QkCircuitInstruction inst;
                    qk_circuit_get_instruction(evolution_circuit, i, &inst);

                    // Prepare qubit vector
                    std::vector<uint32_t> qubits(inst.num_qubits);
                    for (uint32_t j = 0; j < inst.num_qubits; j++) {
                        qubits[j] = inst.qubits[j];
                    }

                    // Get operation kind
                    QkOperationKind kind = qk_circuit_instruction_kind(evolution_circuit, i);

                    // Get mutable pointer to circuit
                    QkCircuit* mutable_circuit = qc.get_rust_circuit().get();

                    // Append based on operation type
                    if (kind == QkOperationKind_Gate) {
                        // Convert gate name string to QkGate enum
                        std::string gate_name(inst.name);
                        auto gate_it = name_map.find(gate_name);
                        if (gate_it != name_map.end()) {
                            QkGate gate_enum = gate_it->second.gate_map();
                            qk_circuit_parameterized_gate(
                                mutable_circuit,
                                gate_enum,           // Use QkGate enum, not string
                                qubits.data(),
                                inst.params
                            );
                        } else {
                            std::cerr << "Unknown gate in evolution circuit: " << gate_name << std::endl;
                            qk_circuit_instruction_clear(&inst);
                            qk_circuit_free(evolution_circuit);
                            return 1;
                        }
                    } else if (kind == QkOperationKind_Barrier) {
                        qk_circuit_barrier(
                            mutable_circuit,
                            qubits.data(), 
                            inst.num_qubits
                        );
                    } else if (kind == QkOperationKind_Reset) {
                        qk_circuit_reset(
                            mutable_circuit,
                            qubits[0]
                        );
                    }

                    qk_circuit_instruction_clear(&inst);
                }

                qk_circuit_free(evolution_circuit);
            }
        }

        circuits.push_back(qc);

        // Display progress
        if ((circuit_idx + 1) % 10 == 0) {
            std::cout << "  Created " << (circuit_idx + 1) << "/" 
                    << num_circuits << " circuits" << std::endl;
        }
    }

    std::cout << "✓ Created all " << num_circuits << " Suzuki-Trotter circuits" << std::endl;

Output:

  Created 10/100 circuits
  Created 20/100 circuits
  Created 30/100 circuits
  Created 40/100 circuits
  Created 50/100 circuits
  Created 60/100 circuits
  Created 70/100 circuits
  Created 80/100 circuits
  Created 90/100 circuits
  Created 100/100 circuits
✓ Created all 100 Suzuki-Trotter circuits

O programa completo então imprime estatísticas dos circuitos e os índices dos operadores do circuito 0, que são omitidos aqui.

Para cada operador amostrado, cria um circuito de evolução de Suzuki-Trotter com tempo $t_k = \frac{\lambda \tau}{n}$, em que $n$ é o número de operadores amostrados por circuito e $\tau$ é o passo de tempo. Como cada grupo amostrado foi normalizado termo a termo antes do mapeamento, essa evolução aplica apenas a fase/sinal de cada termo durante a rotação; as magnitudes contribuem somente por meio da distribuição de amostragem. O estado inicial de Hartree-Fock é preparado aplicando portas X aos qubits 0-6 e 10-16 (sete elétrons em cada setor de spin). As instruções de cada circuito de evolução são anexadas manualmente ao circuito principal usando a API C++.

Etapa 3: Executar usando primitivas do Qiskit (C++)

std::cout << "\n Connecting to IBM Quantum Cloud..." << std::endl;

// Initialize IBM Quantum Compute Service using qiskit-cpp
Qiskit::service::QiskitRuntimeService service;
std::cout << "✓ Connected to IBM Quantum" << std::endl;

// Get backend
const std::string backend_name = "ibm_fez";
auto backend = service.backend(backend_name);

if (backend.name().empty()) {
    std::cerr << "Backend " << backend_name
              << " is not available to this account." << std::endl;
    return 1;
}

std::cout << "✓ Selected backend: " << backend.name() << std::endl;

// Create sampler primitive
const int32_t shots = 100;  // Number of shots per circuit
Qiskit::primitives::BackendSamplerV2 sampler(backend, shots);

std::cout << "\n Submitting " << num_circuits << " circuits ("
            << shots << " shots each)..." << std::endl;

// Submit one job per circuit.
std::vector<std::shared_ptr<Qiskit::primitives::BasePrimitiveJob>> jobs;

for (int i = 0; i < num_circuits; i++) {
    for (uint32_t q = 0; q < num_qubits; q++) {
        circuits[i].measure(q, q);
    }
    Qiskit::circuit::QuantumCircuit transpiled_qc =
        Qiskit::compiler::transpile(circuits[i], backend, 2, 1.0, 42);
    std::vector<Qiskit::primitives::SamplerPub> pubs;
    pubs.push_back(Qiskit::primitives::SamplerPub(transpiled_qc, shots));
    jobs.push_back(sampler.run(pubs));
}

std::cout << "\n⏳ Waiting for " << num_circuits << " jobs to complete..." << std::endl;

// Collect results — poll each job and extract bitstrings
std::vector<boost::dynamic_bitset<>> all_bitstrings;

for (int i = 0; i < num_circuits; i++) {
    while (!jobs[i]->in_final_state()) {
        sleep(5);
    }
    if (jobs[i]->done()) {
        auto result = jobs[i]->result();
        if (result.size() > 0) {
            auto& pub_result = result[0];
            auto& bit_array  = pub_result.data();
            for (size_t s = 0; s < bit_array.num_shots(); s++) {
                auto sample = bit_array[s].to_string();
                boost::dynamic_bitset<> bs(num_qubits);
                for (uint32_t bit = 0; bit < num_qubits; bit++) {
                    bs[bit] = (sample[num_qubits - 1 - bit] == '1');
                }
                all_bitstrings.push_back(bs);
            }
        }
    } else {
        std::cout << "  Circuit " << (i+1) << " did not complete successfully" << std::endl;
    }
}

std::cout << "\n✓ Collected " << all_bitstrings.size()
            << " total bitstrings from all circuits" << std::endl;

Conecta-se ao IBM Quantum, seleciona o backend pelo nome, transpila cada circuito para esse backend e envia um job por circuito: 100 jobs separados do Sampler de 100 shots cada, em um registrador de 20 qubits. Em seguida, consulta cada job até que ele termine e coleta os bitstrings de result[0]. Cada job consome tempo de QPU e espera separadamente na fila, então espere que esta etapa leve muito mais tempo do que um único job em lote.

Os resultados de hardware citados nesta página vêm de uma execução em ibm_pittsburgh (100 circuitos × 100 shots), com backend_name alterado do padrão ibm_fez. As amostras salvas foram pós-processadas com a ordenação de bits de determinantes corrigida mostrada na Etapa 4. O ruído do hardware varia entre execuções, então suas contagens e sua energia final serão diferentes.

Output:

 Connecting to IBM Quantum Cloud...
✓ Connected to IBM Quantum
✓ Selected backend: ibm_pittsburgh

 Submitting 100 circuits (100 shots each)...

⏳ Waiting for 100 jobs to complete...

✓ Collected 10000 total bitstrings from all circuits

Alterar o backend: Altere backend_name em SqDRIFT.cpp para usar outro dispositivo; service.backends() lista o que sua conta consegue acessar.

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

Pós-selecionar bitstrings e gravar strings CI

// Postselect bitstrings by Hamming weight (derived from FCIDump NELEC)
std::vector<double> bitstring_weights(all_bitstrings.size(), 1.0);
auto [filtered_bitstrings, filtered_weights] = Qiskit::addon::sqd::postselect_bitstrings(
    all_bitstrings,
    bitstring_weights,
    Qiskit::addon::sqd::MatchesRightLeftHamming<uint32_t>(n_alpha, n_beta)
);
std::ignore = filtered_weights;

std::cout << "✓ Postselected " << filtered_bitstrings.size()
          << " bitstrings with Hamming weight ("
          << n_alpha << "," << n_beta << ")" << std::endl;

// Convert bitstrings to CI strings using SQD addon
auto ci_strings = Qiskit::addon::sqd::bitstrings_to_ci_strings_symmetrize_spin(
    filtered_bitstrings,
    std::nullopt  // No dimension limit
);

std::cout << "✓ Generated " << ci_strings.size() << " CI strings" << std::endl;

// Write CI strings to file for SBD
std::ofstream alpha_file("alphadets_from_sqd.txt");
for (const auto& ci_string : ci_strings) {
    std::string bitstr;
    // SBD expects the highest orbital index first (bit zero on the right).
    boost::to_string(ci_string, bitstr);
    alpha_file << bitstr << "\n";
}
alpha_file.close();

std::cout << "✓ Wrote " << ci_strings.size()
          << " CI strings to alphadets_from_sqd.txt" << std::endl;
std::cout << " Ready for SBD diagonalization!" << std::endl;

Output:

✓ Postselected 1137 bitstrings with Hamming weight (7,7)
✓ Generated 49 CI strings
✓ Wrote 49 CI strings to alphadets_from_sqd.txt
 Ready for SBD diagonalization!

Pós-seleciona os bitstrings com o peso de Hamming correto (n_alpha elétrons de spin para cima, n_beta de spin para baixo, derivados do campo NELEC do FCIDump), converte-os em strings de Interação de Configurações (CI) usando simetrização de spin e os grava em alphadets_from_sqd.txt para a Diagonalização de Base Selecionada subsequente.

Liberar recursos

// Cleanup
// circuits vector will be automatically cleaned up (RAII)
// Free qubit operators
for (uint32_t i = 0; i < num_groups; i++) {
    qk_obs_free(qubit_ops[i]);
}
delete[] qubit_ops;
// Free fermionic operators
for (uint32_t i = 0; i < num_groups; i++) {
    qf_ferm_op_free(group_ops[i]);
}
delete[] group_ops;
qf_ferm_op_free(normal_ordered);
qf_ferm_op_free(hamiltonian);
qf_fcidump_free(fcidump);
return 0;

Libera corretamente toda a memória e todos os recursos alocados. O std::vector<QuantumCircuit> do C++ é limpo automaticamente via RAII. Libera manualmente os operadores de qubit, os operadores fermiônicos e os dados do FCIDump.

Agora que criamos a base sobre a qual projetar nosso Hamiltoniano, podemos prosseguir com o processo de diagonalização para obter a estimativa do estado fundamental. Execute os comandos a seguir a partir da raiz do projeto, docs/tutorials/assets/sqdrift/cpp.

Compilar o binário diag do SBD

O CMake busca o SBD apenas como código-fonte em deps/sbd/ (a compilação principal intencionalmente não o compila). O app SBD usa um Makefile simples, e o próprio SBD requer MPI, OpenMP e BLAS/LAPACK — portanto, um simples make não basta até que eles estejam instalados e o arquivo Configuration do app corresponda à sua toolchain.

O Makefile lê seus flags de compilador e de link do arquivo Configuration no mesmo diretório. A versão distribuída upstream é ajustada para um layout específico de macOS/Homebrew e não compila como está em uma máquina padrão — principalmente porque -I/opt/homebrew/opt/llvm/include não contém omp.h (ele fica em um diretório de recursos do clang específico da versão), e o clang++ da Apple rejeita -fopenmp de imediato:

clang++: error: unsupported option '-fopenmp'

A receita a seguir foi verificada e funciona no macOS (Apple Silicon, Darwin 25.6, Open MPI 5.x + Homebrew LLVM).

macOS (Apple Silicon) — testado

# 1. Install MPI and an OpenMP-capable compiler.
#    Apple clang cannot do OpenMP, so Homebrew's LLVM provides both
#    clang++ and libomp. BLAS/LAPACK come from Apple's Accelerate
#    framework, which is already part of macOS — no openblas needed.
brew install open-mpi llvm

Em seguida, substitua deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization/Configuration por:

# Path to the SBD library
SBD_PATH=../..

# MPI C++ compiler. OMPI_CXX points Open MPI's wrapper at Homebrew's clang++,
# which (unlike Apple clang) supports -fopenmp.
CCCOM=OMPI_CXX=/opt/homebrew/opt/llvm/bin/clang++ mpicxx

# Build flags. omp.h lives in a version-specific clang resource dir, so let
# Homebrew's clang++ find it itself rather than hardcoding an -I path.
CCFLAGS= -std=c++17 -fopenmp -O3

# Link flags: OpenMP runtime from Homebrew LLVM, BLAS/LAPACK from the
# Accelerate framework that ships with macOS (no openblas install needed).
SYSLIB= -L/opt/homebrew/opt/llvm/lib -lomp -framework Accelerate

# 2. Build
cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization
make
cd ../../../..

Uma compilação bem-sucedida imprime as duas linhas de compilação/link e deixa um binário diag ao lado do Makefile. Um ld: warning: ignoring duplicate libraries: '-lomp' é inofensivo.

macOS Intel: O prefixo do Homebrew é /usr/local em vez de /opt/homebrew, então substitua-o nos dois caminhos anteriores.

Linux (Ubuntu/Debian) — não testado

Os pré-requisitos a seguir cobrem MPI, OpenMP e BLAS/LAPACK. O GCC suporta -fopenmp nativamente, então nenhuma substituição de compilador é necessária. Isto não foi verificado no Linux:

sudo apt install -y libopenmpi-dev libomp-dev libblas-dev liblapack-dev

SBD_PATH=../..
CCCOM=mpicxx
CCFLAGS= -std=c++17 -fopenmp -O3
SYSLIB= -llapack -lblas

Executar o SqDRIFT para gerar as strings CI

Se você já executou ./SqDRIFT nas Etapas 1–4, alphadets_from_sqd.txt já existe na raiz do projeto e você pode pular esta etapa. Caso contrário, execute-o a partir da raiz do projeto:

./SqDRIFT
# Writes alphadets_from_sqd.txt in the project root

Executar a diagonalização

cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization

# The FCIDump lives in the shared sqdrift/fcidump_files/ directory, one level
# above this tutorial's project root, hence five levels up from here.
ln -sf ../../../../../fcidump_files/N2_sto_3g fcidump.txt

# alphadets_from_sqd.txt is written by ./SqDRIFT into the project root itself,
# which is four levels up.
ln -sf ../../../../alphadets_from_sqd.txt alphadets.txt

./diag \
--fcidump fcidump.txt \
--adetfile alphadets.txt \
--method 0 \
--iteration 100 \
--block 10 \
--tolerance 1e-8 \
| tee ../../../../sbd_output.txt

cd ../../../..

Parâmetros:

• --fcidump: Caminho para o arquivo FCIDump que contém o Hamiltoniano
• --adetfile: Caminho para o arquivo de determinantes alfa (strings CI do SqDRIFT)
• --method 0: Método de diagonalização de Davidson
• --iteration 100: Número máximo de iterações de Davidson
• --block 10: Tamanho do bloco para o algoritmo de Davidson
• --tolerance 1e-8: Tolerância de convergência da energia

 Elapsed time for helper construction 0.001512 (sec)
 Elapsed time for init 1e-06 (sec)
 Davidson iteration 0.0 (tol=0.1857493937033264): -107.493531425221
 Davidson iteration 0.1 (tol=0.0455699976495391): -107.5098150493379 -105.93368397386
 Davidson iteration 0.2 (tol=0.009291923190816917): -107.5108356028703 -106.1245977311231 -105.6823366883346
 Davidson iteration 0.3 (tol=0.001488674426419356): -107.5108711084406 -106.3219246016818 -105.7376912083761 -105.4392974334491
 Davidson iteration 0.4 (tol=0.0003188697654737439): -107.5108722237947 -106.5935152393201 -106.1108017566682 -105.4932621297472
 Davidson iteration 0.5 (tol=3.080273116687244e-05): -107.5108722587463 -106.7245930484005 -106.2313070085243 -105.6423237790271
 Davidson iteration 0.6 (tol=2.981303155550973e-06): -107.5108722590151 -106.8008116251279 -106.330411722324 -105.69489630447
 Davidson iteration 0.7 (tol=4.565075388900609e-07): -107.5108722590186 -106.8090569201597 -106.3345525819528 -106.0578325555555
 Davidson iteration 0.8 (tol=4.870689101320012e-08): -107.5108722590186 -106.8098939189322 -106.3973756614146 -106.1617813840531
 Davidson iteration 0.9 (tol=6.224330513865847e-09): -107.5108722590187 -106.8159869231416 -106.4315891597535 -106.2222021964013
 Elapsed time for davidson 0.031276 (sec)
 Elapsed time for diagonalization 0.031286 (sec)
 Elapsed time for mult 0.002511 (sec)
 Energy = -107.5108722590187
 Elapsed time for measurement 4.6e-05 (sec)
 Sample-based diagonalization: Energy = -107.5108722590187
 Sample-based diagonalization: density = [1.999989340109303,1.99999662423008,1.999199121034164,1.981735389885314,1.999825046204946,1.999847567448985,1.995347419300126,0.01136512984757364,0.009844780833945199,0.002849581105572152
 Sample-based diagonalization: carryover bitstrings = [], size = 0

A linha Energy = -107.5108722590187 é a estimativa do estado fundamental sobre o subespaço amostrado.

Como isto difere do exemplo em Python

Os dois exemplos diferem em vários aspectos ao mesmo tempo, não apenas na forma de tratar bitstrings com simetria incorreta, portanto não se espera que produzam energias idênticas:

	Aba C++	Aba Python
Bitstrings com simetria incorreta	Descartados por postselect_bitstrings com base no peso de Hamming	Reparados por recuperação iterativa de configurações
Estrutura do laço	Passagem única: amostrar → pós-selecionar → diagonalizar uma vez	Laço externo: diagonalizar, ler as ocupações orbitais, recuperar configurações, rediagonalizar
Número de circuitos	100 circuitos (num_circuits = 100)	400 circuitos (200 por tempo de evolução)
Tempos de evolução	Um único $\tau$ (time_step = 1)	Dois tempos, times = [1.0, 10.0]
Termos diagonais	Mantidos no operador amostrado	Removidos por filter_diagonal_terms, alterando a distribuição de amostragem e os circuitos
Ordenação dos termos	Apenas agrupados	canonical_order aplicado antes do agrupamento para reprodutibilidade
Solucionador de diagonalização	SBD (diag, binário externo MPI/OpenMP)	Solucionador em processo qiskit-addon-sqd
Backend	Fixo pelo nome (backend_name = "ibm_fez"; a execução desta página usou ibm_pittsburgh)	Backend operacional least_busy com qubits suficientes
Transpilação	Cada circuito separadamente, nível de otimização 2	Todos os circuitos juntos, nível de otimização 3
Envio de jobs	100 jobs, um circuito cada	Um job do Sampler contendo todos os 400 circuitos
Execução em simulador	Nenhuma; apenas hardware	Exemplo com simulador sem ruído antes da execução em hardware

Ambas as versões usam o mesmo Hamiltoniano, agrupam seus termos por estrutura eletrônica, ponderam cada grupo pelo coeficiente absoluto médio de seus termos, partem do estado de Hartree-Fock e amostram 10 grupos e 100 shots por circuito.

Cada diferença afeta o subespaço amostrado, e portanto a energia, de forma independente. O conjunto maior e o segundo tempo de evolução, mais longo, ampliam as configurações que o exemplo em Python explora; amostrar em mais de um tempo dá ao SqDRIFT sua cobertura semelhante à de Krylov do espaço de baixa energia. Filtrar termos diagonais altera quais operadores o amostrador pode sortear e seus pesos, de modo que as duas execuções amostram distribuições diferentes.

A pós-seleção é um filtro estrito, e um shot com o número de partículas errado é descartado por completo. Na execução acima, isso descartou a maior parte dos dados: 1137 de 10000 shots sobreviveram à verificação de peso de Hamming (7,7) (cerca de 11%), reduzindo-se a apenas 49 $\alpha$-determinantes distintos. A recuperação de configurações, por sua vez, usa as ocupações orbitais médias de uma diagonalização anterior para inverter bits e reparar esses shots, reaproveitando o que a pós-seleção descartaria. Esse feedback a torna inerentemente iterativa: recover_configurations em qiskit-addon-sqd-hpc requer um argumento avg_occupancies que só existe depois que uma diagonalização foi executada.

Esta execução em C++ resulta em um subespaço de 49 $\alpha$-determinantes e uma energia de -107.5109 Ha. O espaço completo de $\alpha$-determinantes para sete elétrons em 10 orbitais tem apenas $\binom{10}{7} = 120$ strings, então fornecer todas as 120 ao mesmo binário diag dá a energia FCI exata -107.6482 Ha. A execução, portanto, fica cerca de 0.1373 Ha acima do exato, abrangendo 49 dos 120 $\alpha$-determinantes. Como o SQD é variacional no subespaço selecionado, um subespaço menor só pode elevar a energia, de modo que o subespaço incompleto é a leitura mais direta da diferença.

Essa verificação com FCI exato só é possível porque o exemplo é pequeno; ela não faz parte do fluxo de trabalho do SqDRIFT e, para tamanhos de sistema em que o SQD é realmente necessário, não existe tal referência. Aqui ela serve apenas para quantificar quanto do espaço esta execução recuperou.

Uma comparação com o exemplo em Python não é controlada, portanto uma diferença não deve ser atribuída apenas à ausência de recuperação de configurações; observe que o número de circuitos, o tempo de evolução único, a filtragem de termos diagonais e o ruído do hardware diferem todos ao mesmo tempo e todos moldam o subespaço.

Para fechar essa diferença, aqui estão alguns próximos passos possíveis para tentar depois de concluir este tutorial: aumente shots, num_circuits ou ops_per_circuit; amostre em mais de um time_step; e adicione um laço de recuperação sobre Qiskit::addon::sqd::recover_configurations.

Próximos passos

Recomendações

• 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 químico — um tutorial sobre como construir um circuito Jastrow de cluster unitário local (LUCJ) para simulação de química quântica.

• O artigo do SqDRIFT — a literatura em que 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 com base na evolução das bibliotecas utilizadas.)