SqDRIFT algoritmo para estimativa do estado fundamental

Estimativa de tempo de execução: 180 segundos em um processador Heron r3 (OBSERVAÇÃO: trata-se apenas de uma estimativa. (O tempo de execução pode variar.)

Está procurando a versão em C++?

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

Resultados do aprendizado

Aprenda a criar circuitos com menor profundidade em comparação com a técnica de Trotterização
Conheça um fluxo de trabalho completo para estimativa do estado fundamental utilizando o “ qDRIFT ” e o SQD
Aprenda a usar qiskit-fermions em conjunto com outros complementos do Qiskit para implementar esse fluxo de trabalho

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

Pré-requisitos

Leia a visão geral sobre a diagonalização quântica baseada em amostras (SQD)
Leia a aula sobre a Diagonalização Quântica de Krylov Baseada em Amostras (SKQD)

Segundo plano

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

Este tutorial utiliza o Qiskit Fermions para criar circuitos fermiónicos mais naturais para o algoritmo “ qDRIFT ”, seguido pela aplicação de etapas de layout e síntese fermiónicas antes de integrar os circuitos ao pipeline tradicional do Qiskit para execução em hardware.

Seja o hamiltoniano da forma:

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

onde, sem perda de generalidade, exigimos que $c_i > 0$ e que o maior valor próprio de $h_i$ seja igual, em valor absoluto, a $1$ . Qualquer prefator com sinal ou complexo é absorvido por $h_i$ , de modo que os coeficientes $c_i$ são pesos estritamente positivos, enquanto os $h_i$ determinam a direção de cada termo. Aqui, $N$ é o número de termos (ou, após o agrupamento, o número de grupos) no hamiltoniano; trata-se de uma propriedade do hamiltoniano e é diferente do número de operadores incluídos em um único circuito, denotado por $n$ a seguir.

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

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 é calculado sobre os sorteios $n$ , e não sobre todos os termos hamiltonianos $N$ ; e, como os termos são sorteados com reposição, o mesmo $h_i$ pode aparecer mais de uma vez em um único $V_k$ .

A quantidade:

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

é a norma de $L_1$ dos coeficientes; assim, cada uma das etapas $n$ evolui durante o mesmo intervalo de tempo $\lambda t / n$ , independentemente do termo que tenha sido sorteado. A uniformidade do ângulo de passo é a característica marcante de qDRIFT: : um coeficiente influencia o resultado pela frequência com que seu termo é extraído, e não pelo grau de rotação desse termo. Os índices são extraídos da distribuição:

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

portanto, a série $(k_1, \ldots, k_n)$ é uma sequência aleatória de índices de termos extraídos dessa distribuição. Como os $c_i$ são positivos e somam $\lambda$ , trata-se de uma distribuição de probabilidade normalizada, e a esperança do canal resultante ao longo dos sorteios aleatórios se aproxima da evolução sob $H$ , com um erro que diminui à medida que $n$ cresce. Observe que o erro de aproximação depende de $\lambda$ e não do número de termos $N$ .

(O artigo “ SqDRIFT ” refere-se ao número de termos como “ $\mathcal{N}$ ” e ao comprimento da sequência como “ $N$ ”; usamos aqui “ $N$ ” e “ $n$ ” para manter as duas noções claramente distintas.)

Este tutorial mostra como gerar um conjunto desses circuitos aleatórios. Depois de criarmos esses circuitos, da mesma forma que criamos um subespaço de Krylov para diferentes operadores, amostramos sequências de bits a partir de vários desses operadores com diferentes parâmetros de tempo. Isso garante uma maior sobreposição entre os vetores do estado fundamental e as sequências de bits amostradas.

Requisitos

Antes de iniciar este tutorial, certifique-se de ter instalado

Um ambiente virtual do 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 o seguinte comando:

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

Instalaçã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,
)

Exemplo de simulador

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

Como ler e preparar o FCIDump

Para este tutorial, vamos carregar o hamiltoniano da estrutura eletrônica do nitrogênio ( N2 ). Existem também outras maneiras de criar operadores fermiónicos. Consulte a documentação em qiskit_fermions.operators.library.

Sobre este FCIDump. O arquivo descreve N2_sto_3g uma molécula de nitrogênio ( $N_2$ ) na base mínima STO-3G, com uma separação interatômica de 1.09 $\AA$ , o comprimento de ligação de equilíbrio experimental. Seu cabeçalho declara NORB=10, NELEC=14, e MS2=0: 10 orbitais espaciais (portanto, 20 orbitais de spin e 20 qubits na representação de Jordan-Wigner), 14 elétrons em um singuleto de spin, ou seja, sete elétrons em $\alpha$ e sete em $\beta$ . Todos os orbitais recebem o índice de simetria 1, ou seja, não é utilizada nenhuma simetria de grupo pontual. Por se tratar de um dump de STO-3G em espaço completo, nenhum orbital está congelado e o espaço de correlação é pequeno o suficiente para que uma energia de referência FCI exata possa ser calculada classicamente para fins de comparação, conforme mostrado na próxima célula.

É possível regenerar um arquivo equivalente com o comando 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 convergentes, um arquivo regenerado pode diferir do arquivo original quanto à fase ou à ordem dos orbitais; as energias totais não são afetadas.

Obtenção do arquivo. Encontre o FCIDump neste repositório: GitHub. Você pode executar a célula abaixo para importá-la para o local esperado pelo restante do tutorial.

Primeiro, usamos a função fornecida cisolver pelo pyscf para obter a energia de referência. Essa é a verdadeira energia do estado fundamental da molécula com a qual estamos trabalhando. Para isso, vamos primeiro declarar e norb nelec, que correspondem ao número de orbitais e ao número de elétrons, respectivamente. Em seguida, declaramos e h1e h2e, que são, respectivamente, os integrais de um e de dois elétrons. Todos esses dados também serão utilizados posteriormente para o SQD.

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

Output:

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

Output:

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 preparados, lemos o hamiltoniano do arquivo FCI em um formato compatível com qiskit-fermions

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

Fluxos de trabalho de férmions com qiskit-fermions

Primeiramente, mapearemos o hamiltoniano para um modelo de circuito fermiónico utilizando qiskit-fermions, que oferece passagens de transpiler e portas específicas para circuitos fermiónicos. Esses elementos serão utilizados posteriormente, antes das etapas tradicionais de transpilagem do Qiskit para este fluxo de trabalho.

Agrupamento de termos

Para garantir a reprodutibilidade dos resultados, primeiro utilizamos canonical_order para ordenar os termos com base apenas em sua estrutura. A ordem dos operadores na lista canon é, portanto, fixa. Isso garante a reprodutibilidade dos operadores criados, pois a função pass QDriftTrotterization que utilizaremos nas próximas amostras seleciona índices aleatórios para criar os operadores qDRIFT.

Nesta etapa, aproveitamos as diversas simetrias presentes no hamiltoniano da estrutura eletrônica, agrupando termos relacionados com coeficientes idênticos. Embora isso altere a distribuição dos coeficientes dos operadores a partir da qual o protocolo de qDRIFT e realiza suas amostragens, isso não afeta suas garantias de convergência. Fundamentalmente, o agrupamento de termos relacionados por simetria resulta em uma compensação favorável dos termos de Pauli e em uma profundidade de circuito globalmente menor ao se realizar a evolução temporal de um estado sob a ação desses termos.

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

Observe que a fórmula group_terms_by_electronic_structure pressupõe termos com ordem normal.

Filtragem de termos diagonais

Removemos os termos diagonais do hamiltoniano usado para gerar os circuitos, de modo que $n$ qDRIFT os intervalos de amostragem sejam dedicados aos termos que movimentam a população entre as configurações. É melhor filtrar esses termos do hamiltoniano neste momento, antes que o Evolution portão seja construído na próxima etapa.

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

o deslocamento de energia constante, um produto de operadores de número zero, cuja evolução temporal contribui apenas com uma fase global;
os operadores de número individual $n_i$ , cuja evolução temporal se reduz a rotações de um único qubit $Z$ ;
os produtos de ordem superior, como $n_i n_j$ .

Por si só, nenhuma dessas ações transfere população entre configurações de ocupação e número; elas atuam apenas nas fases das configurações já existentes. No entanto, elas não são inertes: essas fases relativas contribuem para a interferência gerada pelos termos de excitação mais adiante no circuito; portanto, filtrá-las altera a evolução que é efetivamente gerada e pode alterar a distribuição amostral. Trata-se de uma aproximação deliberada na etapa de geração do circuito, feita para concentrar a amostragem nos termos de excitação, e não de uma etapa que deixe a distribuição amostrada inalterada. Ao contrário do agrupamento de simetria acima, que mantém intactas as garantias de convergência do método “ qDRIFT ”, esse filtro altera o operador que está sendo evoluído. Portanto, os circuitos não se aproximam mais da evolução sob o hamiltoniano completo, e os limites de erro de “ qDRIFT ” se aplicam ao operador filtrado, e não ao original. Isso é aceitável neste contexto porque os circuitos são apenas uma heurística de amostragem usada para propor configurações: nenhum termo é perdido na própria estimativa de energia, uma vez que o filtro se aplica apenas ao hamiltoniano utilizado para construir os circuitos, enquanto a diagonalização clássica, realizada posteriormente, utiliza o hamiltoniano completo, incluindo os termos diagonais. A precisão do SQD depende desse passo clássico, que permanece variacional no subespaço amostrado, independentemente de como as configurações foram propostas.

A função filter_diagonal_terms() remove esses termos de um operador diretamente. Ela os identifica a partir de sua estrutura ordenada normalmente — o multiconjunto de modos de criação correspondendo ao multiconjunto de modos de aniquilação —, portanto, só é válida para um operador que já esteja ordenado normalmente. Essa suposição não é verificada durante a execução.

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

Output:

Agora que agrupamos os termos no hamiltoniano, vamos definir os seguintes parâmetros para gerar o conjunto de circuitos:

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

Criação de circuitos fermiónicos

Agora, vamos criar circuitos fermiónicos para cada um dos intervalos de tempo. Cada circuito será composto por um único portão de evolução, com o tempo de evolução que definimos anteriormente. O operador de evolução é o hamiltoniano. Posteriormente, executamos etapas de transpilagem nesses circuitos para criar circuitos d qDRIFT.

Preparação do Ansatz

Preparamos o estado de Hartree-Fock utilizando a classe InitializeModes . Para o nitrogênio, o processo consiste simplesmente em aplicar portas X aos primeiros qubits num_elec_a e, em seguida, aos num_elec_b qubits, sendo que, no caso do nitrogênio, ambos os valores são iguais a sete. Esse estado representa os sete elétrons $\alpha$ e os sete elétrons $\beta$ do nitrogênio.

# 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, vamos primeiro usar as etapas disponíveis em qiskit-fermions para realizar otimizações no nível fermionico e, em seguida, fazer a transpilagem do nosso circuito para o backend de nossa escolha. Como se trata de um experimento em simulador, vamos fazer isso primeiro para o AerSimulator.

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. A etapa do transpiler qDRIFT faz isso por nós. Agora podemos criar circuitos mais simples que podem ser executados no hardware de forma mais eficiente, apesar da conectividade limitada dos qubits, mesmo quando o hamiltoniano contém acoplamentos de longo alcance e termos superiores ao quadrático. Após o agrupamento de termos, ele seleciona os operadores com base em seus pesos. Para cada operador $h_i$ , o peso $W_{h_i}$ é definido da seguinte forma:

W_{h_i} = |c_i| / \lambda

Otimizações fermiónicas e nativas do hardware

A função retorna generate_preset_jw_pass_manager() um objeto MultiStagePassManager que recebe um argumento e FermionicCircuit gera um circuito final otimizado, que podemos transpilá-lo para ser executado em nosso hardware. Substituímos sua etapa de otimização padrão por uma que FermionicPassManager contenha nossa passagem QDriftTrotterization :

A etapa QDriftTrotterization utiliza internamente o cálculo de pesos e a amostragem para gerar os circuitos que usaremos para a amostragem
A etapa RelabelModes é mais uma etapa de otimização que pode ser usada para permutar os modos fermiónicos, a fim de otimizar a conectividade entre os qubits e reduzir a profundidade das portas; leia mais na referência da API

As etapas restantes são executadas MultiStagePassManager automaticamente e cuidam de todo o mapeamento de férmions para qubits:

F2QLayout : O gerenciador de passagens predefinido aplica a TrivialF2QLayout passagem, que mapeia de forma trivial os bits fermiônicos $n$ para os qubits $n$ .
F2QSynth : Uma etapa de transpilagem para mapear instruções de circuitos baseados 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))

Output:

Agora que concluímos as otimizações no nível dos férmions, podemos transpilá-los para execução no simulador.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Etapa 3: Executar usando o comando `Qiskit primitives`

Agora que já temos nossos circuitos, podemos executá-los usando o Qiskit primitives no site AerSimulator. Vamos somar todas as contagens dos diferentes circuitos. Nós os convertemos em vetores booleanos antes de, por fim, realizar o pós-processamento com o SQD.

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

Output:

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

Etapa 4: Realizar o pós-processamento e apresentar o resultado no formato clássico desejado

Utilização de cadeias de bits para SQD

Agora podemos aplicar o esquema de diagonalização às sequências de bits selecionadas para encontrar o menor valor próprio que corresponderá à energia do estado fundamental da molécula. Criamos uma função de retorno de chamada, declaramos as ocupações iniciais e definimos os parâmetros antes de, finalmente, executar o esquema de diagonalização. A função de retorno é usada para exibir a iteração atual e a estimativa atual do valor próprio a cada iteração.

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

Observação : a dimensão do subespaço não é fixa entre as iterações, mesmo no simulador sem ruído — cada subamostra gera um conjunto diferente de configurações, e a etapa de recuperação reestrutura o conjunto entre as iterações; portanto, a dimensão relatada varia de uma subamostra para outra. A amostragem silenciosa, por si só, não determina a dimensão do subespaço selecionado. A execução no hardware, no entanto, tende a gerar subespaços sistematicamente maiores, pois as imagens com ruído quebram a simetria do número de partículas e a recuperação da configuração as transforma em vetores de base adicionais. Por isso, também apresentaremos outra etapa para a poda de cadeias de bits na seção de hardware.

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

Output:

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 utiliza 20 qubits (10 orbitais espaciais). Essa escolha é uma medida de conveniência para um tutorial que deve ser executado rapidamente, e não um limite rígido para o método.

O custo da etapa clássica não é determinado diretamente pelo número de qubits. O SQD diagonaliza o hamiltoniano projetado no subespaço gerado pelas configurações amostradas; portanto, o que determina o custo clássico é a dimensão desse subespaço selecionado — determinada, neste caso, por samples_per_batch, num_batches, e pelo número de configurações distintas que os circuitos realmente produzem — juntamente com a álgebra linear esparsa necessária para aplicar o hamiltoniano projetado. O espaço completo de CI cresce combinatoriamente com os orbitais e os elétrons, mas o subespaço selecionado é uma pequena parcela ajustável desse espaço, e controlamos seu tamanho diretamente. Consequentemente, o número de qubits e a dificuldade clássica podem ser variados de forma relativamente independente: um espaço orbital mais amplo, amostrado em um subespaço modesto, pode ser mais econômico do que um sistema menor diagonalizado em um espaço orbital muito grande.

Na prática, portanto, o tamanho viável do sistema depende da dimensão do subespaço necessária para a precisão desejada e da memória e dos núcleos disponíveis para o solucionador de valores próprios. Espaços orbitais maiores geralmente exigem um subespaço maior para se alcançar precisão química, e é isso que, em última instância, motiva o uso de recursos distribuídos — consulte o qiskit-addon-sqd-hpc para saber como ampliar essa etapa. Em vez de definir um limite fixo, a abordagem prática consiste em acompanhar a dimensão do subespaço relatada e a convergência da energia ao longo das iterações, aumentando o tamanho do subespaço até que a energia pare de melhorar ou até esgotar a memória disponível.

Observação: Devido ao erro de amostragem causado pelo ruído no hardware, o subespaço criado para a diagonalização na execução no hardware será maior do que o obtido ao usar o simulador. Embora aumente a dimensão do subespaço que desejamos diagonalizar, o fluxo de trabalho ainda nos fornece uma resposta precisa devido à robustez do SQD em relação ao ruído.

Eliminação de strings espúrias

Aqui, podemos optar por realizar uma etapa adicional. Quando tivermos todas as sequências de bits resultantes das execuções do circuito, podemos filtrar as sequências inválidas antes de executar o SQD ou prosseguir sem fazer a poda. Geralmente, é preferível não realizar a poda em execuções de hardware, pois isso mantém as tentativas com simetria quebrada disponíveis para a recuperação da configuração, que pode repará-las, transformando-as em configurações válidas e, assim, ampliar o subespaço, em vez de descartar essas tentativas imediatamente.

Como o nitrogênio só pode ter sete $\alpha$ e sete $\beta$ elétrons, quaisquer sequências de bits que tenham mais ou menos do que sete 1s na primeira e na segunda metade da saída podem ser descartadas. Definimos uma função que verifica se as sequências de bits são válidas e, caso não sejam, as descarta. Depois de filtrarmos as sequências de bits espúrias, o restante é enviado para o esquema de diagonalização. Use o sinalizador PRUNE abaixo para alternar entre os dois comportamentos.

Lembre-se de que a poda é apenas uma das várias opções que definem o subespaço final, juntamente com o número de circuitos, o conjunto de tempos de evolução e a filtragem de termos diagonais. Comparar uma execução podada com uma não podada só é informativo se todos os outros fatores forem mantidos constantes; o guia complementar do C++ aborda esse assunto com mais detalhes, já que ele utiliza pós-seleção em vez de recuperação e também difere nesses outros parâmetros.

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

Output:

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óximas etapas

Recomendações

Se você achou este trabalho interessante, talvez se interesse pelo material a seguir:

Diagonalização quântica de Krylov baseada em amostras de um modelo de rede fermiónica — um tutorial relacionado que utiliza 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 em química quântica.
O artigo “ SqDRIFT ” — a literatura na qual este tutorial se baseia. (Observe que algumas das otimizações discutidas neste artigo ainda estão em desenvolvimento, e este tutorial está sujeito a alterações no futuro, de acordo com a evolução das bibliotecas utilizadas.)

Esta página foi útil?

Relate um bug, erro de digitação ou solicite conteúdo no GitHub.

SqDRIFT algoritmo para estimativa do estado fundamental

Resultados do aprendizado

Pré-requisitos

Segundo plano

Requisitos

Instalaçã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 o comando Qiskit primitives

Etapa 4: Realizar o pós-processamento e apresentar o resultado no formato clássico desejado

Exemplo de hardware

Próximas etapas

Etapa 3: Executar usando o comando `Qiskit primitives`