Algoritmo SqDRIFT para la estimación del estado fundamental

Python

Descarga el cuaderno de Python

Descarga el cuaderno de Python para ejecutar este tutorial por tu cuenta. El cuaderno descarga la entrada FCIDump si aún no está presente.

Estimación de uso: 180 segundos en un procesador Heron r3 (NOTA: esto es solo una estimación. Tu tiempo de ejecución puede variar.)

Resultados del aprendizaje

• Aprende a crear circuitos de menor profundidad en comparación con la trotterización

• Recorre un flujo de trabajo de extremo a extremo para la estimación del estado fundamental usando qDRIFT y SQD

• Aprende a usar qiskit-fermions junto con otros addons de Qiskit para implementar dicho flujo de trabajo

Prerrequisitos

• Lee la introducción a la diagonalización cuántica basada en muestras (SQD)

• Lee la lección sobre diagonalización cuántica de Krylov basada en muestras (SKQD)

Antecedentes

SqDRIFT es una variante de SKQD que sustituye la necesidad de elegir un ansatz del que muestrear bitstrings por un conjunto de circuitos de evolución temporal construidos directamente a partir del hamiltoniano objetivo. Esto se logra submuestreando operadores de evolución temporal más pequeños del hamiltoniano según sus coeficientes, lo que se conoce como el método de trotterización qDRIFT.

Ambas implementaciones usan Qiskit Fermions para representar el hamiltoniano molecular y construir circuitos para el algoritmo qDRIFT. La implementación en Python usa pasadas de transpilador fermiónicas, mientras que la implementación en C++ muestrea grupos de forma explícita y los asigna a operadores de qubits.

Sea el hamiltoniano de la forma:

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

donde, sin pérdida de generalidad, exigimos $c_i > 0$ y que el autovalor más grande de $h_i$ sea igual, en valor absoluto, a $1$. Cualquier prefactor con signo o complejo se absorbe en $h_i$, así que los coeficientes $c_i$ son pesos estrictamente positivos mientras que los $h_i$ llevan la dirección de cada término. Aquí $N$ es el número de términos (o, tras la agrupación, el número de grupos) en el hamiltoniano; es una propiedad del hamiltoniano y es distinto del número de operadores muestreados en un único circuito, escrito $n$ a continuación.

El algoritmo qDRIFT entonces realiza, para el tiempo objetivo $t$, algún operador $V_k$, donde $k$ va desde $1 \cdots K$ y denota el $k$-ésimo circuito SqDRIFT, definido como:

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

Aquí $n$ es el número de operadores muestreados por circuito y $K$ es el número de circuitos en el conjunto. El producto se realiza sobre las $n$ extracciones, no sobre los $N$ términos del hamiltoniano, y debido a que los términos se extraen con reemplazo, el mismo $h_i$ puede aparecer más de una vez en un único $V_k$.

La cantidad:

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

es la norma $L_1$ de los coeficientes, así que cada uno de los $n$ pasos evoluciona durante la misma duración $\lambda t / n$ sin importar qué término se extrajo. La uniformidad del ángulo de paso es la característica distintiva de qDRIFT: un coeficiente influye en el resultado a través de con qué frecuencia se extrae su término, no a través de cuánto se rota ese término. Los índices se muestrean a partir de la distribución:

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

así que la serie $(k_1, \ldots, k_n)$ es una secuencia aleatoria de índices de términos extraída de esta distribución. Dado que los $c_i$ son positivos y suman $\lambda$, esta es una distribución de probabilidad normalizada, y la esperanza del canal resultante sobre las extracciones aleatorias aproxima la evolución bajo $H$, con un error que disminuye a medida que crece $n$. Ten en cuenta que el error de aproximación depende de $\lambda$ en lugar de depender del número de términos $N$.

(El artículo de SqDRIFT escribe el número de términos como $\mathcal{N}$ y la longitud de la secuencia como $N$; aquí usamos $N$ y $n$ para mantener las dos cosas claramente distintas.)

Este tutorial muestra cómo generar un conjunto de dichos circuitos aleatorizados. Después de haber creado estos circuitos, de forma similar a como creamos un subespacio de Krylov para diferentes operadores, muestreamos bitstrings de múltiples operadores de este tipo con diferentes parámetros temporales. Esto garantiza un mayor solapamiento entre los vectores del estado fundamental y los bitstrings muestreados.

Requisitos

Antes de comenzar este tutorial, asegúrate de haber instalado

• Un entorno virtual de Python (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (ten en cuenta que el nombre está en plural)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Puedes instalar todos los paquetes necesarios con:

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

Configuración

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

Implementación

Paso 1: Mapear las entradas clásicas a un problema cuántico

Los pasos en Python que siguen usan el simulador Aer. El ejemplo en hardware viene después.

Lectura y preparación del FCIDump

Para este tutorial, cargaremos el hamiltoniano de estructura electrónica para el nitrógeno (N2). También hay otras formas de crear operadores fermiónicos. Consulta la documentación en qiskit_fermions.operators.library.

Acerca de este FCIDump. El archivo N2_sto_3g describe una molécula de nitrógeno ($N_2$) en la base mínima STO-3G con una separación interatómica de 1.09 $\AA$, la longitud de enlace de equilibrio experimental. Su encabezado declara NORB=10, NELEC=14 y MS2=0: 10 orbitales espaciales (por lo tanto 20 espín-orbitales, y 20 qubits bajo Jordan-Wigner), 14 electrones en un singlete de espín, es decir, siete electrones $\alpha$ y siete $\beta$. A todos los orbitales se les asigna la etiqueta de simetría 1, es decir, no se aprovecha ninguna simetría de grupo puntual. Al tratarse de un volcado STO-3G de espacio completo, no se congela ningún orbital y el espacio de correlación es lo suficientemente pequeño como para poder calcular clásicamente una energía de referencia FCI exacta para comparación, como se muestra en el siguiente bloque de código.

Se puede regenerar un archivo equivalente con 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")

Dado que las integrales dependen de los orbitales SCF convergidos, un archivo regenerado podría diferir del proporcionado en la fase u orden de los orbitales; las energías totales no se ven afectadas.

Obtención del archivo. Encuentra el FCIDump en este repositorio de GitHub. Puedes ejecutar el código de abajo para obtenerlo en la ubicación que el resto del tutorial espera.

Primero usamos el cisolver proporcionado por pyscf para obtener la energía de referencia. Esta es la verdadera energía del estado fundamental de la molécula con la que estamos trabajando. Para esto, primero declararemos norb y nelec, que son el número de orbitales y el número de electrones, respectivamente. Luego declaramos h1e y h2e, que son las integrales de uno y dos electrones respectivamente. Todas estas también se usarán más adelante para 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

Cargando el Hamiltoniano

Con los datos necesarios listos, leemos el Hamiltoniano del archivo FCI en un formato compatible con qiskit-fermions

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

Flujos de trabajo fermiónicos con qiskit-fermions

Primero mapearemos el hamiltoniano en un modelo de circuito fermiónico usando qiskit-fermions, que proporciona pasadas de transpilador y puertas específicas para circuitos fermiónicos. Estas se usarán más adelante antes de las pasadas de transpilador tradicionales de Qiskit para este flujo de trabajo.

Agrupación de términos

Para garantizar la reproducibilidad de los resultados, primero usamos canonical_order para ordenar los términos basándonos únicamente en su estructura. El orden de los operadores en la lista canon queda, por lo tanto, fijo. Esto garantiza la reproducibilidad de los operadores creados porque la pasada QDriftTrotterization que usaremos más adelante muestrea índices aleatorios para crear los operadores qDRIFT.

En este paso, aprovechamos las numerosas simetrías presentes en el Hamiltoniano de estructura electrónica agrupando términos relacionados con coeficientes idénticos. Si bien hacer esto cambia la distribución de coeficientes del operador de la que muestrea el protocolo qDRIFT, esto no afecta sus garantías de convergencia. Fundamentalmente, agrupar términos relacionados por simetría produce una cancelación favorable de términos de Pauli y, en general, una profundidad de circuito más corta al evolucionar un estado bajo su acción.

qiskit-fermions proporciona la función group_terms_by_electronic_structure que realiza esta agrupación por nosotros.

Ten en cuenta que group_terms_by_electronic_structure asume términos en orden normal.

Filtrado de términos diagonales

Eliminamos los términos diagonales del Hamiltoniano utilizado para generar los circuitos, de modo que las $n$ ranuras de muestreo qDRIFT se dediquen a términos que mueven población entre configuraciones. Es mejor filtrar dichos términos del Hamiltoniano en este punto, antes de que se construya la compuerta Evolution en el siguiente paso.

Los términos en cuestión son los que son diagonales en la base de número de ocupación, es decir, los productos de operadores de número $a^\dagger_i a_i$. Tres tipos de términos entran en esta descripción:

• el desplazamiento de energía constante, un producto de cero operadores de número, cuya evolución temporal solo aporta una fase global;

• los operadores de número individuales $n_i$, cuya evolución temporal se reduce a rotaciones $Z$ de un solo qubit;

• los productos de orden superior como $n_i n_j$.

Por sí solos, ninguno de estos mueve población entre configuraciones de número de ocupación; solo actúan sobre las fases de las configuraciones ya presentes. Sin embargo, no son inertes: esas fases relativas alimentan la interferencia generada por los términos de excitación más adelante en el circuito, por lo que filtrarlos cambia la evolución que realmente se genera y puede cambiar la distribución de muestreo. Esta es una aproximación deliberada en el paso de generación de circuitos, hecha para centrar el muestreo en los términos de excitación, en lugar de un paso que deja intacta la distribución muestreada. A diferencia de la agrupación por simetría anterior, que deja intactas las garantías de convergencia de qDRIFT, este filtro cambia el operador que se está evolucionando. Por lo tanto, los circuitos ya no aproximan la evolución bajo el Hamiltoniano completo, y los límites de error de qDRIFT se aplican al operador filtrado en lugar de al original. Esto es aceptable aquí porque los circuitos son solo una heurística de muestreo utilizada para proponer configuraciones: no se pierde ningún término de la estimación de energía en sí, ya que el filtro se aplica solo al Hamiltoniano utilizado para construir los circuitos, mientras que la diagonalización clásica posterior utiliza el Hamiltoniano completo, incluidos los términos diagonales. La precisión de SQD depende de ese paso clásico, que sigue siendo variacional en el subespacio muestreado independientemente de cómo se hayan propuesto las configuraciones.

La función filter_diagonal_terms() elimina dichos términos de un operador en el mismo lugar (in place). Los identifica a partir de su estructura en orden normal, el multiconjunto de modos de creación que coincide con el multiconjunto de modos de aniquilación, por lo que solo es válida en un operador que ya está en orden normal. Esta suposición no se verifica en tiempo de ejecución.

# 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

Ahora que hemos agrupado los términos en el Hamiltoniano, decidiremos los siguientes parámetros para generar el conjunto de circuitos:

• El número de circuitos a generar: num_circuits
• La longitud de cada circuito en términos de grupos de excitación: num_exc
• El factor para los diferentes tiempos de evolución: times

Creación de circuitos fermiónicos

Ahora crearemos circuitos fermiónicos para cada uno de los pasos temporales. Cada circuito consistirá en una única puerta de evolución, con el tiempo de evolución que declaramos anteriormente. El operador de evolución es el hamiltoniano. Más adelante ejecutaremos pasadas de transpilador en estos circuitos para crear circuitos qDRIFT.

Preparación del ansatz

Preparamos el estado de Hartree-Fock usando la clase InitializeModes. Para el nitrógeno, el proceso consiste simplemente en aplicar compuertas X a los primeros num_elec_a qubits y luego a los num_elec_b qubits, ambos iguales a siete para el nitrógeno. Este estado representa los siete electrones $\alpha$ y siete $\beta$ del nitrógeno.

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

Paso 2: Optimizar el problema para la ejecución en hardware cuántico

Ahora que tenemos nuestros circuitos, primero usaremos las pasadas disponibles en qiskit-fermions para realizar optimizaciones a nivel fermiónico, seguidas de la transpilación de nuestro circuito para el backend elegido. Dado que este es un experimento de simulador, primero haremos esto para el AerSimulator.

Cálculo de peso para cada grupo

En este paso, realizamos el muestreo qDRIFT de términos de manera estocástica con probabilidades proporcionales a sus coeficientes en el hamiltoniano. La pasada de transpilador qDRIFT hace esto por nosotros. Ahora podemos crear circuitos menos profundos que pueden ejecutarse en el hardware de manera más eficiente a pesar de la conectividad limitada de los qubits, incluso cuando el hamiltoniano contiene acoplamientos de largo alcance y términos de orden superior al cuadrático. Después de la agrupación de términos, muestrea los operadores según sus pesos. Para cada operador $h_i$, el peso $W_{h_i}$ se define de la siguiente manera:

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

Como los términos se agruparon en el Paso 1, cada $h_i$ aquí es un grupo completo: $c_i$ es el coeficiente absoluto medio de los términos del grupo $i$, y cada término del grupo evoluciona con su coeficiente reducido a su signo.

Optimizaciones fermiónicas y nativas del hardware

La función generate_preset_jw_pass_manager() devuelve un MultiStagePassManager que toma un FermionicCircuit y produce un circuito final optimizado que podemos transpilar para ejecutar en nuestro hardware. Reemplazamos su etapa de optimización predeterminada por un FermionicPassManager que contiene nuestra pasada QDriftTrotterization:

• La pasada QDriftTrotterization utiliza internamente el cálculo de pesos y el muestreo para generar los circuitos que usaremos para el muestreo

• La pasada RelabelModes es otra pasada de optimización que puede usarse para permutar los modos fermiónicos y optimizar la conectividad entre qubits y reducir la profundidad de puertas; lee más en la referencia de la API

Las etapas restantes del MultiStagePassManager se ejecutan automáticamente y gestionan el mapeo completo de fermiones a qubits:

• F2QLayout: El gestor de pasadas preestablecido aplica la pasada TrivialF2QLayout, que mapea de manera trivial $n$ bits fermiónicos a $n$ qubits.

• F2QSynth: Una pasada de transpilación para mapear instrucciones de circuito basadas en fermiones a instrucciones basadas en 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

Ahora que hemos terminado con las optimizaciones a nivel fermiónico, podemos transpilar los circuitos para su ejecución en el simulador.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Paso 3: Ejecutar usando primitivas de Qiskit

Ahora que tenemos nuestros circuitos, podemos ejecutarlos usando primitivas de Qiskit en el AerSimulator. Combinaremos todos los conteos de los diferentes circuitos. Los convertimos en vectores booleanos antes de finalmente posprocesarlos con 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

Paso 4: Posprocesar y devolver el resultado en el formato clásico deseado

Usando bitstrings para SQD

Ahora podemos ejecutar el esquema de diagonalización en los bitstrings seleccionados para encontrar el autovalor más bajo, que corresponderá a la energía del estado fundamental de la molécula. Creamos una función de callback, declaramos ocupaciones iniciales y establecemos los parámetros antes de finalmente ejecutar el esquema de diagonalización. La función de callback se usa para imprimir la iteración actual y la estimación actual del autovalor en cada iteración.

Finalmente, para obtener la estimación del estado fundamental, agregamos la nuclear_repulsion_energy a la energía resultante.

Nota: La dimensión del subespacio no es fija entre iteraciones, incluso en el simulador sin ruido; cada submuestra extrae un conjunto diferente de configuraciones, y el paso de recuperación reorganiza el conjunto entre iteraciones, por lo que la dimensión reportada varía de una submuestra a otra. El muestreo sin ruido no fija por sí mismo la dimensión del subespacio seleccionado. La ejecución en hardware, sin embargo, tiende a dar subespacios sistemáticamente más grandes, porque las mediciones ruidosas rompen la simetría de número de partículas y la recuperación de configuraciones las convierte en vectores base adicionales. Debido a esto, también introduciremos otro paso para podar bitstrings en la sección 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

Ejemplo de hardware

Este ejemplo utiliza 20 qubits (10 orbitales espaciales). Esa elección es una conveniencia para un tutorial que debe ejecutarse rápidamente, no un límite estricto del método.

El costo del paso clásico no está determinado directamente por el número de qubits. SQD diagonaliza el Hamiltoniano proyectado sobre el subespacio abarcado por las configuraciones muestreadas, por lo que lo que determina el costo clásico es la dimensión de ese subespacio seleccionado, gobernada aquí por samples_per_batch, num_batches, y cuántas configuraciones distintas producen realmente los circuitos, junto con el álgebra lineal dispersa necesaria para aplicar el Hamiltoniano proyectado. El espacio CI completo crece combinatoriamente con los orbitales y los electrones, pero el subespacio seleccionado es una porción pequeña y ajustable de él, y controlamos su tamaño directamente. En consecuencia, el número de qubits y la dificultad clásica pueden variarse hasta cierto punto de forma independiente: un espacio orbital más amplio muestreado en un subespacio modesto puede ser más barato que un sistema más pequeño diagonalizado sobre uno muy grande.

En la práctica, entonces, el tamaño de sistema factible depende de la dimensión del subespacio que necesites para la precisión que deseas y de la memoria y los núcleos disponibles para el resolutor de autovalores. Los espacios orbitales más grandes típicamente sí requieren un subespacio más grande para alcanzar la precisión química, y eso es lo que eventualmente motiva los recursos distribuidos; consulta qiskit-addon-sqd-hpc para escalar este paso. En lugar de suponer un límite fijo, el enfoque práctico es observar la dimensión del subespacio reportada y la convergencia de la energía a lo largo de las iteraciones, y aumentar el tamaño del subespacio hasta que la energía deje de mejorar o se agote la memoria disponible.

Nota: Debido al error de muestreo por el ruido en el hardware, el subespacio creado para la diagonalización en la ejecución de hardware será más grande que el que obtenemos al usar el simulador. Si bien esto aumenta la dimensión del subespacio que queremos diagonalizar, el flujo de trabajo igualmente nos da una respuesta precisa gracias a la robustez de SQD frente al ruido.

Poda de cadenas espurias

Aquí podemos elegir realizar un paso adicional. Cuando tenemos todos los bitstrings de las ejecuciones del circuito, podemos filtrar los bitstrings inválidos antes de ejecutar SQD, o continuar sin podar. Omitir la poda es generalmente preferible para las ejecuciones en hardware, porque deja disponibles las mediciones con simetría rota para la recuperación de configuraciones, que puede repararlas en configuraciones válidas y así ampliar el subespacio en lugar de descartar esas mediciones directamente.

Dado que el nitrógeno solo puede tener siete electrones $\alpha$ y siete $\beta$, se puede descartar cualquier bitstring que tenga más o menos de siete 1s en la primera y la segunda mitad de la salida. Definimos una función que verifica si los bitstrings son válidos, y si no lo son, los descarta. Una vez que filtramos los bitstrings espurios, el resto se envía al esquema de diagonalización. Usa la bandera PRUNE a continuación para alternar entre los dos comportamientos.

Ten en cuenta que la poda es solo una de varias opciones que dan forma al subespacio final, junto con el número de circuitos, el conjunto de tiempos de evolución y el filtrado de términos diagonales. Comparar una ejecución podada con una sin podar solo es informativo si todo lo demás se mantiene fijo; la pestaña de C++ analiza esto con más detalle, ya que realiza una posselección en lugar de una recuperación y también difiere en esos otros 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++

Obtén el código en C++

El código fuente en C++, el proyecto de CMake y las instrucciones de compilación están en el directorio docs/tutorials/assets/sqdrift/cpp del repositorio de documentación de Qiskit. La entrada FCIDump está en el directorio vecino fcidump_files. Clona el repositorio, o descarga ambos directorios, para compilar y ejecutar este tutorial por tu cuenta.

Estimación de uso: 180 segundos en un procesador Heron r3 (NOTA: esto es solo una estimación. Tu tiempo de ejecución puede variar.)

Resultados del aprendizaje (C++)

• Aprende a crear circuitos de menor profundidad en comparación con la trotterización

• Recorre un flujo de trabajo de extremo a extremo para la estimación del estado fundamental usando qDRIFT y SQD

• Aprende a usar qiskit-fermions junto con otros addons de Qiskit para implementar dicho flujo de trabajo

Prerrequisitos (C++)

• Lee la introducción a la diagonalización cuántica basada en muestras (SQD)

• Lee la lección sobre diagonalización cuántica de Krylov basada en muestras (SKQD)

Antecedentes (C++)

SqDRIFT es una variante de SKQD que sustituye la necesidad de elegir un ansatz del que muestrear bitstrings por un conjunto de circuitos de evolución temporal construidos directamente a partir del hamiltoniano objetivo. Esto se logra submuestreando operadores de evolución temporal más pequeños del hamiltoniano según sus coeficientes, lo que se conoce como el método de trotterización qDRIFT.

Ambas implementaciones usan Qiskit Fermions para representar el hamiltoniano molecular y construir circuitos para el algoritmo qDRIFT. La implementación en Python usa pasadas de transpilador fermiónicas, mientras que la implementación en C++ muestrea grupos de forma explícita y los asigna a operadores de qubits.

Sea el hamiltoniano de la forma:

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

donde, sin pérdida de generalidad, exigimos $c_i > 0$ y que el autovalor más grande de $h_i$ sea igual, en valor absoluto, a $1$. Cualquier prefactor con signo o complejo se absorbe en $h_i$, así que los coeficientes $c_i$ son pesos estrictamente positivos mientras que los $h_i$ llevan la dirección de cada término. Aquí $N$ es el número de términos (o, tras la agrupación, el número de grupos) en el hamiltoniano; es una propiedad del hamiltoniano y es distinto del número de operadores muestreados en un único circuito, escrito $n$ a continuación.

El algoritmo qDRIFT entonces realiza, para el tiempo objetivo $t$, algún operador $V_k$, donde $k$ va desde $1 \cdots K$ y denota el $k$-ésimo circuito SqDRIFT, definido como:

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

Aquí $n$ es el número de operadores muestreados por circuito y $K$ es el número de circuitos en el conjunto. El producto se realiza sobre las $n$ extracciones, no sobre los $N$ términos del hamiltoniano, y debido a que los términos se extraen con reemplazo, el mismo $h_i$ puede aparecer más de una vez en un único $V_k$.

La cantidad:

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

es la norma $L_1$ de los coeficientes, así que cada uno de los $n$ pasos evoluciona durante la misma duración $\lambda t / n$ sin importar qué término se extrajo. La uniformidad del ángulo de paso es la característica distintiva de qDRIFT: un coeficiente influye en el resultado a través de con qué frecuencia se extrae su término, no a través de cuánto se rota ese término. Los índices se muestrean a partir de la distribución:

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

así que la serie $(k_1, \ldots, k_n)$ es una secuencia aleatoria de índices de términos extraída de esta distribución. Dado que los $c_i$ son positivos y suman $\lambda$, esta es una distribución de probabilidad normalizada, y la esperanza del canal resultante sobre las extracciones aleatorias aproxima la evolución bajo $H$, con un error que disminuye a medida que crece $n$. Ten en cuenta que el error de aproximación depende de $\lambda$ en lugar de depender del número de términos $N$.

(El artículo de SqDRIFT escribe el número de términos como $\mathcal{N}$ y la longitud de la secuencia como $N$; aquí usamos $N$ y $n$ para mantener las dos cosas claramente distintas.)

Este tutorial muestra cómo generar un conjunto de dichos circuitos aleatorizados. Después de haber creado estos circuitos, de forma similar a como creamos un subespacio de Krylov para diferentes operadores, muestreamos bitstrings de múltiples operadores de este tipo con diferentes parámetros temporales. Esto garantiza un mayor solapamiento entre los vectores del estado fundamental y los bitstrings muestreados.

Requisitos (C++)

Esta implementación se ejecuta en hardware de IBM Quantum® y usa un ejecutable externo de SBD para la diagonalización clásica. No incluye el ejemplo con simulador en Python.

Sigue las instrucciones de compilación de C++ para instalar el compilador, CMake, Rust, Boost y las dependencias de Qiskit. Instala Rust con rustup: las dependencias de Qiskit fijan su propia versión de Rust, que rustup descarga automáticamente, pero un Rust instalado desde Homebrew o desde un gestor de paquetes de Linux podría ser demasiado antiguo. El paso final de SBD requiere además MPI, OpenMP y BLAS/LAPACK. Windows no es compatible actualmente con esta implementación.

Usa el código fuente en C++ completo y el proyecto de CMake del repositorio. Conserva el directorio adyacente fcidump_files al descargar el ejemplo. Desde una copia del repositorio de documentación, compila con los siguientes comandos. La primera compilación descarga y compila todas las dependencias y tarda varios minutos.

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

La compilación escribe el ejecutable SqDRIFT en docs/tutorials/assets/sqdrift/cpp (la raíz del proyecto), no en build/. Ejecuta ./SqDRIFT desde la raíz del proyecto para que pueda encontrar ../fcidump_files/N2_sto_3g. Todos los demás comandos de esta pestaña también se ejecutan desde la raíz del proyecto, salvo que se indique lo contrario. Si ./SqDRIFT falla con un error de biblioteca compartida ausente, sigue las instrucciones de ruta de bibliotecas específicas de tu plataforma en la guía de compilación; en macOS, el ejecutable normalmente funciona sin ellas.

Configuración (C++)

Los fragmentos siguientes son extractos del programa SqDRIFT.cpp completo. Compila y ejecuta ese archivo como un único programa.

Antes de ejecutar SqDRIFT, asegúrate de que tu token de API de IBM Quantum esté configurado para IBM Cloud Platform.

Si tienes Python con qiskit-ibm-runtime instalado, la forma más sencilla de crear un archivo de configuración con el formato correcto es:

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

Esto escribe ~/.qiskit/qiskit-ibm.json (Linux/macOS) o %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) con la estructura correcta. Si prefieres crear el archivo manualmente:

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

Obtén tu clave de API en IBM Quantum Platform.

Nota: Un archivo de configuración ausente o con formato incorrecto provoca un pánico poco útil de Rust unwrap() al iniciar, en lugar de un error claro. Revisa bien el nombre de la clave (default-ibm-quantum-platform) y el valor del canal (ibm_quantum_platform) si encuentras algún problema.

Implementación (C++)

Paso 1: Mapear las entradas clásicas a un problema cuántico (C++)

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

Carga el hamiltoniano de la molécula $N_2$ desde un archivo FCIDump. El número de qubits es igual al doble del número de orbitales (espín arriba y espín abajo).

Crea el hamiltoniano y aplica el orden normal

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

Convierte los datos de FCIDump en un operador fermiónico y aplica el orden normal para simplificar la estructura del operador.

Agrupa los términos por estructura electró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 los términos relacionados del hamiltoniano por estructura electrónica antes de construir sus circuitos de evolución.

Calcula los pesos de muestreo, luego normaliza y asigna a operadores de qubits

Aquí el orden importa. La normalización reescala cada coeficiente a magnitud unitaria in situ, por lo que los pesos de muestreo deben calcularse primero. Si la normalización se ejecuta primero, cada $|c_j|$ se lee como exactamente 1.0 y el peso de cada grupo se reduce a un simple recuento de sus términos: la distribución de muestreo quedaría entonces determinada por el tamaño del grupo y no por los coeficientes físicos.

El peso de cada grupo es el coeficiente absoluto medio de sus términos. La pasada QDriftTrotterization de la pestaña de Python usa el mismo peso para un hamiltoniano agrupado. Usar la suma haría que $\lambda$ escalara con el tamaño del grupo, lo que a su vez sesgaría el tiempo de evolución $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

El peso $\lambda_i = \frac{1}{|G_i|}\sum_{j \in G_i} |c_{ij}|$ es el coeficiente absoluto medio sobre los términos del grupo $G_i$, leído de los coeficientes originales. La normalización reescala entonces cada término fermiónico a magnitud unitaria antes de la transformación de Jordan-Wigner, conservando solo su fase/signo en el operador evolucionado. Esto coincide con la evolución de términos normalizados de qDRIFT: las magnitudes de los coeficientes determinan las probabilidades de muestreo mediante $\lambda$, mientras que la evolución del circuito usa el término normalizado para que los coeficientes grandes no se cuenten dos veces.

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

Extrae los conjuntos de operadores para 100 circuitos aleatorios, cada uno con 10 grupos muestreados de la distribución ponderada. El siguiente paso convierte estos conjuntos en el conjunto de circuitos qDRIFT.

Paso 2: Construir los circuitos qDRIFT (C++)

En la versión en C++, este paso construye los circuitos a partir de los conjuntos de operadores muestreados. La transpilación para el backend de destino ocurre en el Paso 3, inmediatamente antes de enviar cada circuito.

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

El programa completo imprime después estadísticas de los circuitos y los índices de operadores del circuito 0, que aquí se omiten.

Para cada operador muestreado, crea un circuito de evolución de Suzuki-Trotter con tiempo $t_k = \frac{\lambda \tau}{n}$, donde $n$ es el número de operadores muestreados por circuito y $\tau$ es el paso de tiempo. Como cada grupo muestreado se normalizó término a término antes de la asignación, esta evolución aplica solo la fase/signo de cada término durante la rotación; las magnitudes contribuyen únicamente a través de la distribución de muestreo. El estado inicial de Hartree-Fock se prepara aplicando compuertas X a los qubits 0-6 y 10-16 (siete electrones en cada sector de espín). Las instrucciones de cada circuito de evolución se añaden manualmente al circuito principal con la API de C++.

Paso 3: Ejecutar usando primitivas de 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;

Se conecta a IBM Quantum, selecciona el backend por nombre, transpila cada circuito para ese backend y envía un trabajo por circuito: 100 trabajos de Sampler independientes de 100 mediciones cada uno en un registro de 20 qubits. Luego consulta cada trabajo hasta que termina y recoge las cadenas de bits de result[0]. Cada trabajo consume tiempo de QPU y espera por separado en la cola, así que es de esperar que este paso tarde mucho más que un único trabajo por lotes.

Los resultados de hardware citados en esta página provienen de una ejecución en ibm_pittsburgh (100 circuitos × 100 mediciones), con backend_name cambiado respecto al valor predeterminado ibm_fez. Las muestras guardadas se posprocesaron con el orden de bits de los determinantes corregido que se muestra en el Paso 4. El ruido del hardware varía entre ejecuciones, por lo que tus recuentos y tu energía final serán distintos.

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

Cambiar el backend: Cambia backend_name en SqDRIFT.cpp para apuntar a otro dispositivo; service.backends() enumera los que tu cuenta puede alcanzar.

Paso 4: Posprocesar y devolver el resultado en el formato clásico deseado (C++)

Posselecciona las cadenas de bits y escribe las cadenas 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!

Posselecciona las cadenas de bits con el peso de Hamming correcto (n_alpha electrones con espín arriba, n_beta con espín abajo, derivados del campo NELEC de FCIDump), las convierte en cadenas de Interacción de Configuraciones (CI) mediante simetrización de espín y las escribe en alphadets_from_sqd.txt para la posterior Diagonalización de Base Seleccionada.

Libera los 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 correctamente toda la memoria y los recursos asignados. El std::vector<QuantumCircuit> de C++ se limpia automáticamente mediante RAII. Libera manualmente los operadores de qubits, los operadores fermiónicos y los datos de FCIDump.

Ahora que hemos creado la base sobre la que proyectar nuestro hamiltoniano, podemos continuar con el proceso de diagonalización para obtener la estimación del estado fundamental. Ejecuta los siguientes comandos desde la raíz del proyecto, docs/tutorials/assets/sqdrift/cpp.

Compila el binario diag de SBD

CMake obtiene SBD solo como código fuente en deps/sbd/ (la compilación principal omite compilarlo a propósito). La aplicación SBD usa un Makefile simple, y el propio SBD requiere MPI, OpenMP y BLAS/LAPACK — por lo que un simple make no basta hasta que estén instalados y el archivo Configuration de la aplicación coincida con tu cadena de herramientas.

El Makefile lee sus indicadores de compilación y de enlace del archivo Configuration del mismo directorio. La versión distribuida en el repositorio original está ajustada a una disposición concreta de macOS/Homebrew y no se compilará tal cual en un equipo estándar — sobre todo porque -I/opt/homebrew/opt/llvm/include no contiene omp.h (está en un directorio de recursos de clang específico de cada versión), y el clang++ de Apple rechaza -fopenmp de plano:

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

La siguiente receta está verificada y funciona en macOS (Apple Silicon, Darwin 25.6, Open MPI 5.x + LLVM de Homebrew).

macOS (Apple Silicon) — probado

# 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

Luego reemplaza 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 ../../../..

Una compilación correcta imprime las dos líneas de compilación/enlace y deja un binario diag junto al Makefile. Un ld: warning: ignoring duplicate libraries: '-lomp' es inofensivo.

macOS con Intel: El prefijo de Homebrew es /usr/local en lugar de /opt/homebrew, así que sustitúyelo en las dos rutas anteriores.

Linux (Ubuntu/Debian) — sin probar

Los siguientes requisitos previos cubren MPI, OpenMP y BLAS/LAPACK. GCC admite -fopenmp de forma nativa, por lo que no hace falta sustituir el compilador. Esto no se ha verificado en 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

Ejecuta SqDRIFT para generar las cadenas CI

Si ya ejecutaste ./SqDRIFT durante los Pasos 1–4, alphadets_from_sqd.txt ya existe en la raíz del proyecto y puedes omitir este paso. En caso contrario, ejecútalo desde la raíz del proyecto:

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

Ejecuta la diagonalización

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: ruta al archivo FCIDump que contiene el hamiltoniano
• --adetfile: ruta al archivo de determinantes alfa (cadenas CI de SqDRIFT)
• --method 0: método de diagonalización de Davidson
• --iteration 100: número máximo de iteraciones de Davidson
• --block 10: tamaño de bloque del algoritmo de Davidson
• --tolerance 1e-8: tolerancia de convergencia de la energía

 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

La línea Energy = -107.5108722590187 es la estimación del estado fundamental sobre el subespacio muestreado.

En qué se diferencia del ejemplo en Python

Los dos ejemplos difieren en varios aspectos a la vez, no solo en cómo tratan las cadenas de bits con mala simetría, por lo que no se espera que produzcan energías idénticas:

	Pestaña de C++	Pestaña de Python
Cadenas de bits con mala simetría	Descartadas por postselect_bitstrings según el peso de Hamming	Reparadas mediante recuperación iterativa de configuraciones
Estructura del bucle	Una sola pasada: muestrear → posseleccionar → diagonalizar una vez	Bucle externo: diagonalizar, leer las ocupaciones orbitales, recuperar configuraciones, volver a diagonalizar
Número de circuitos	100 circuitos (num_circuits = 100)	400 circuitos (200 por tiempo de evolución)
Tiempos de evolución	Un único $\tau$ (time_step = 1)	Dos tiempos, times = [1.0, 10.0]
Términos diagonales	Se conservan en el operador muestreado	Eliminados por filter_diagonal_terms, lo que cambia la distribución de muestreo y los circuitos
Orden de los términos	Solo agrupados	canonical_order aplicado antes de agrupar para lograr reproducibilidad
Solucionador de diagonalización	SBD (diag, binario externo MPI/OpenMP)	Solucionador en proceso de qiskit-addon-sqd
Backend	Fijado por nombre (backend_name = "ibm_fez"; la ejecución de esta página usó ibm_pittsburgh)	Backend operativo least_busy con suficientes qubits
Transpilación	Cada circuito por separado, nivel de optimización 2	Todos los circuitos juntos, nivel de optimización 3
Envío de trabajos	100 trabajos, un circuito cada uno	Un trabajo de Sampler con los 400 circuitos
Ejecución en simulador	Ninguna; solo hardware	Ejemplo con simulador sin ruido antes de la ejecución en hardware

Ambas versiones usan el mismo hamiltoniano, agrupan sus términos por estructura electrónica, ponderan cada grupo por el coeficiente absoluto medio de sus términos, parten del estado de Hartree-Fock y muestrean 10 grupos y 100 mediciones por circuito.

Cada diferencia afecta al subespacio muestreado, y por tanto a la energía, de forma independiente. El conjunto más grande y el segundo tiempo de evolución, más largo, amplían las configuraciones que explora el ejemplo en Python; muestrear en más de un tiempo da a SqDRIFT su cobertura del espacio de baja energía de tipo Krylov. Filtrar los términos diagonales cambia qué operadores puede extraer el muestreador y sus pesos, de modo que las dos ejecuciones muestrean distribuciones distintas.

La posselección es un filtro estricto, y una medición con el número de partículas incorrecto se descarta por completo. En la ejecución anterior eso descartó la mayor parte de los datos: sobrevivieron 1137 de 10000 mediciones a la comprobación del peso de Hamming (7,7) (alrededor del 11 %), que se redujeron a solo 49 $\alpha$-determinantes distintos. La recuperación de configuraciones, en cambio, usa las ocupaciones orbitales medias de una diagonalización previa para invertir bits y reparar esas mediciones, reciclando lo que la posselección descartaría. Esa retroalimentación la hace intrínsecamente iterativa: recover_configurations en qiskit-addon-sqd-hpc requiere un argumento avg_occupancies que solo existe después de ejecutar una diagonalización.

Esta ejecución en C++ produce un subespacio de 49 $\alpha$-determinantes y una energía de -107.5109 Ha. El espacio completo de $\alpha$-determinantes para siete electrones en 10 orbitales es de solo $\binom{10}{7} = 120$ cadenas, así que al pasar las 120 al mismo binario diag se obtiene la energía FCI exacta -107.6482 Ha. Por tanto, la ejecución queda unos 0.1373 Ha por encima del valor exacto abarcando 49 de los 120 $\alpha$-determinantes. Dado que SQD es variacional en el subespacio seleccionado, un subespacio más pequeño solo puede elevar la energía, por lo que el subespacio incompleto es la lectura más directa de la diferencia.

Esta comprobación con FCI exacto solo es posible porque el ejemplo es pequeño; no forma parte del flujo de trabajo de SqDRIFT, y para los tamaños de sistema en los que SQD realmente hace falta no existe tal referencia. Aquí sirve solo para cuantificar qué parte del espacio recuperó esta ejecución.

La comparación con el ejemplo en Python no es controlada, por lo que una diferencia no debe atribuirse únicamente a la ausencia de recuperación de configuraciones; ten en cuenta que el número de circuitos, el único tiempo de evolución, el filtrado de términos diagonales y el ruido del hardware difieren a la vez y todos dan forma al subespacio.

Para cerrar la diferencia, estos son algunos posibles pasos siguientes tras completar este tutorial: aumentar shots, num_circuits u ops_per_circuit; muestrear en más de un time_step; y añadir un bucle de recuperación sobre Qiskit::addon::sqd::recover_configurations.

Próximos pasos

Recomendaciones

• Diagonalización cuántica de Krylov basada en muestras de un modelo de red fermiónica — un tutorial relacionado que usa circuitos de evolución temporal en lugar de un ansatz variacional.

• Diagonalización cuántica basada en muestras de un hamiltoniano de química — un tutorial sobre cómo construir un circuito de clúster unitario local de Jastrow (LUCJ) para la simulación de química cuántica.

• El artículo de SqDRIFT — la literatura en la que se basa este tutorial. (Ten en cuenta que algunas de las optimizaciones comentadas en este artículo están actualmente en desarrollo, y este tutorial puede cambiar en el futuro según la evolución de las bibliotecas utilizadas.)