Algoritmo SqDRIFT para la estimación del estado fundamental

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

Objetivos de aprendizaje

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

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

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

Este tutorial se presenta como un notebook de Python con fines didácticos.

Requisitos previos

• Lee la descripción general de Diagonalización cuántica basada en muestras (SQD)

• Lee la lección de 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 cadenas de bits 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.

Este tutorial usa Qiskit Fermions para crear los circuitos fermiónicos, más naturales, para el algoritmo qDRIFT, seguido del uso de pases de disposición y síntesis fermiónicos antes de introducir los circuitos en el pipeline tradicional de Qiskit para su ejecución en hardware.

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 mayor autovalor de $h_i$ sea igual, en valor absoluto, a $1$. Cualquier prefactor con signo o complejo se absorbe en $h_i$, de modo 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 agrupar, el número de grupos) del Hamiltoniano; es una propiedad del Hamiltoniano y es distinto del número de operadores muestreados en un solo circuito, escrito $n$ más abajo.

El algoritmo qDRIFT realiza entonces, para el tiempo objetivo $t$, algún operador $V_k$, donde $k$ va de $1 \cdots K$ y designa el $k_{th}$ 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 del conjunto. El producto recorre las $n$ extracciones, no los $N$ términos del Hamiltoniano, y como 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, de modo que cada uno de los $n$ pasos evoluciona durante la misma duración $\lambda t / n$ independientemente de qué término se haya extraído. La uniformidad del ángulo de paso es la característica distintiva de qDRIFT: un coeficiente influye en el resultado por la frecuencia con que se extrae su término, no por cuánto se rota ese término. Los índices se muestrean de la distribución:

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

de modo que la serie $(k_1, \ldots, k_n)$ es una secuencia aleatoria de índices de términos extraídos de esta distribución. Como 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 al crecer $n$. Ten en cuenta que el error de aproximación depende de $\lambda$ y no 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 ambos claramente distintos.)

Este tutorial muestra cómo generar un conjunto de tales circuitos aleatorizados. Una vez creados estos circuitos, de forma similar a cómo creamos un subespacio de Krylov para distintos operadores, muestreamos cadenas de bits de varios de estos operadores con distintos parámetros de tiempo. Esto garantiza un mayor solapamiento entre los vectores del estado fundamental y las cadenas de bits muestreadas.

Requisitos

Antes de comenzar este tutorial, asegúrate de tener 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 es 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

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

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

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

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

# Qiskit Aer
from qiskit_aer import AerSimulator

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

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

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

Ejemplo con simulador

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

Lectura y preparación del FCIDump

Para este tutorial cargaremos el Hamiltoniano de estructura electrónica del nitrógeno (N2). Hay otras formas de crear operadores fermiónicos también. 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 a una separación interatómica de 1.09 $\AA$, la longitud de enlace de equilibrio experimental. Su cabecera declara NORB=10, NELEC=14 y MS2=0: 10 orbitales espaciales (por 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 da la etiqueta de simetría 1, es decir, no se explota ninguna simetría de grupo puntual. Al ser un volcado STO-3G de espacio completo, no hay orbitales congelados y el espacio de correlación es lo bastante pequeño como para calcular clásicamente una energía FCI exacta de referencia para comparar, como se muestra en la siguiente celda.

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

Como las integrales dependen de los orbitales SCF convergidos, un archivo regenerado puede diferir del suministrado en la fase o el 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 la celda siguiente para descargarlo en la ubicación que espera el resto del tutorial.

Primero usamos el cisolver que proporciona pyscf para obtener la energía de referencia. Esta es la energía real del estado fundamental de la molécula con la que trabajamos. Para ello declararemos primero norb y nelec, que son el número de orbitales y el número de electrones, respectivamente. Después declaramos h1e y h2e, que son las integrales de uno y dos electrones respectivamente. Todos estos se usarán más adelante también 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}")

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

Carga del 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 asignaremos el Hamiltoniano a un modelo de circuito fermiónico usando qiskit-fermions, que proporciona pases del transpilador y puertas específicos para circuitos fermiónicos. Estos se usarán más adelante antes de los pases tradicionales del transpilador 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 solo en su estructura. El orden de los operadores en la lista canon queda por tanto fijado. Esto garantiza la reproducibilidad de los operadores creados porque el pase QDriftTrotterization que usaremos más adelante muestrea índices aleatorios para crear los operadores qDRIFT.

En este paso aprovechamos las muchas simetrías presentes en el Hamiltoniano de estructura electrónica agrupando términos relacionados con coeficientes idénticos. Aunque al hacerlo cambia la distribución de coeficientes del operador de la que muestrea el protocolo qDRIFT, esto no afecta a sus garantías de convergencia. Es crucial que agrupar términos relacionados por simetría da lugar a una cancelación favorable de términos de Pauli y a una menor profundidad global del circuito al evolucionar temporalmente un estado bajo su acción.

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

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

Filtrado de términos diagonales

Eliminamos los términos diagonales del Hamiltoniano usado para generar los circuitos, de modo que las $n$ ranuras de muestreo de qDRIFT se gasten en términos que mueven población entre configuraciones. Conviene filtrar estos términos del Hamiltoniano en este punto, antes de construir la puerta Evolution en el siguiente paso.

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

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

• los operadores 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 ellos mueve población entre configuraciones de números de ocupación; actúan únicamente 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, de modo que filtrarlos cambia la evolución que realmente se genera y puede cambiar la distribución de muestreo. 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, y no un paso que deje 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 evoluciona. Por tanto, los circuitos ya no aproximan la evolución bajo el Hamiltoniano completo, y las cotas de error de qDRIFT se aplican al operador filtrado y no al original. Esto es aceptable aquí porque los circuitos son solo una heurística de muestreo usada para proponer configuraciones: no se pierde ningún término en la propia estimación de la energía, ya que el filtro se aplica solo al Hamiltoniano usado para construir los circuitos, mientras que la diagonalización clásica posterior usa el Hamiltoniano completo, términos diagonales incluidos. La precisión de SQD depende de ese paso clásico, que sigue siendo variacional en el subespacio muestreado con independencia de cómo se hayan propuesto las configuraciones.

La función filter_diagonal_terms() elimina tales términos de un operador in situ. Los identifica a partir de su estructura en orden normal — el multiconjunto de modos de creación 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 comprueba 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 del 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 distintos tiempos de evolución: times

Creación de circuitos fermiónicos

Ahora crearemos circuitos fermiónicos para cada uno de los pasos de tiempo. Cada circuito constará de una única puerta de evolución, con el tiempo de evolución que declaramos antes. El operador de evolución es el Hamiltoniano. Más adelante ejecutamos pases del transpilador sobre 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 puertas X a los primeros num_elec_a qubits y después 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, usaremos primero los pases disponibles en qiskit-fermions para realizar optimizaciones a nivel fermiónico, seguidos de la transpilación de nuestro circuito para el backend elegido. Como este es un experimento con simulador, lo haremos primero para el AerSimulator. Cálculo de pesos para cada grupo

En este paso, realizamos el muestreo qDRIFT de términos de forma estocástica, con probabilidades proporcionales a sus coeficientes en el hamiltoniano. El pase del Transpiler qDRIFT lo hace por nosotros. Ahora podemos crear circuitos menos profundos que se pueden ejecutar 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. Tras la agrupación de términos, se muestrean 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 ejecutarlo en nuestro hardware. Sustituimos su etapa de optimización predeterminada por un FermionicPassManager que contiene nuestro pase QDriftTrotterization:

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

• El pase RelabelModes es otro pase de optimización que se puede usar para permutar los modos fermiónicos con el fin de optimizar la conectividad entre qubits y reducir la profundidad de las puertas; consulta más información en la referencia de la API

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

• F2QLayout: El gestor de pases preestablecido aplica el pase TrivialF2QLayout, que asigna de forma trivial $n$ bits fermiónicos a $n$ qubits.

• F2QSynth: Un pase de transpilación para convertir instrucciones de circuito basadas en fermiones en 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 ejecutarlos en el simulador.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Paso 3: Ejecutar con las primitivas de Qiskit

Ahora que tenemos nuestros circuitos, podemos ejecutarlos con las primitivas de Qiskit en el AerSimulator. Combinaremos todos los recuentos de los distintos circuitos. Los convertimos en vectores booleanos antes de realizar finalmente el posprocesamiento 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

Uso de cadenas de bits para SQD

Ahora podemos ejecutar el esquema de diagonalización sobre las cadenas de bits seleccionadas para encontrar el menor autovalor, que corresponderá a la energía del estado fundamental de la molécula. Creamos una función de callback, declaramos las ocupaciones iniciales y establecemos los parámetros antes de ejecutar finalmente 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.

Por último, para obtener la estimación del estado fundamental, sumamos nuclear_repulsion_energy a la energía resultante.

Nota: La dimensión del subespacio no es fija entre iteraciones, ni siquiera en el simulador sin ruido: cada submuestra extrae un conjunto distinto de configuraciones y el paso de recuperación remodela 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í solo la dimensión del subespacio seleccionado. Sin embargo, la ejecución en hardware tiende a dar subespacios sistemáticamente mayores, porque los shots con ruido rompen la simetría del número de partículas y la recuperación de configuraciones los convierte en vectores de base adicionales. Por ello, también introduciremos otro paso para podar cadenas de bits 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 en hardware

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

El coste del paso clásico no viene determinado directamente por el número de qubits. SQD diagonaliza el hamiltoniano proyectado sobre el subespacio generado por las configuraciones muestreadas, de modo que lo que determina el coste clásico es la dimensión de ese subespacio seleccionado —regida aquí por samples_per_batch, num_batches y por 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 de forma combinatoria 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 variar de forma algo independiente: un espacio orbital más amplio muestreado en un subespacio modesto puede resultar más barato que un sistema más pequeño diagonalizado en uno muy grande.

En la práctica, entonces, el tamaño de sistema viable depende de la dimensión del subespacio que necesites para la precisión que desees y de la memoria y los núcleos disponibles para el solucionador de autovalores. Los espacios orbitales más grandes suelen requerir un subespacio mayor para alcanzar la precisión química, y eso es lo que acaba motivando el uso de recursos distribuidos; consulta qiskit-addon-sqd-hpc para escalar este paso. En lugar de suponer un límite fijo, el enfoque práctico consiste en 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 agotes la memoria disponible.

Nota: Debido al error de muestreo causado por el ruido del hardware, el subespacio creado para la diagonalización en la ejecución en hardware será mayor que el que obtenemos con el simulador. Aunque aumenta la dimensión del subespacio que queremos diagonalizar, el flujo de trabajo sigue dándonos una respuesta precisa gracias a la robustez de SQD frente al ruido.

Poda de cadenas espurias

Aquí podemos optar por realizar un paso adicional. Cuando tenemos todas las cadenas de bits de las ejecuciones de los circuitos, podemos filtrar las cadenas de bits no válidas antes de ejecutar SQD, o continuar sin podar. Omitir la poda suele ser preferible en las ejecuciones en hardware, porque deja los shots con simetría rota disponibles para la recuperación de configuraciones, que puede repararlos hasta convertirlos en configuraciones válidas y, de ese modo, ampliar el subespacio en lugar de descartar directamente esos shots.

Dado que el nitrógeno solo puede tener siete electrones $\alpha$ y siete $\beta$, se pueden descartar las cadenas de bits que tengan más o menos de siete 1 en la primera y en la segunda mitad de la salida. Definimos una función que comprueba si las cadenas de bits son válidas y, si no lo son, las descarta. Una vez filtradas las cadenas de bits espurias, el resto se envía al esquema de diagonalización. Usa el indicador PRUNE de abajo para alternar entre los dos comportamientos.

Ten en cuenta que la poda es solo una de varias decisiones 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 con poda con otra sin poda solo es informativo si todo lo demás se mantiene fijo; la versión en C++ de este tutorial lo analiza con más detalle, ya que posselecciona en lugar de recuperar y además 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")

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

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

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

Próximos pasos

Recomendaciones

Si este trabajo te ha parecido interesante, quizá te interese el siguiente material:

• 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 químico: un tutorial sobre cómo construir un circuito de Jastrow de clúster unitario local (LUCJ) para la simulación de química cuántica.
• El artículo SqDRIFT: la literatura en la que se basa este tutorial. (Ten en cuenta que algunas de las optimizaciones que se analizan en este artículo son actualmente un trabajo en curso, y este tutorial puede cambiar en el futuro según evolucionen las bibliotecas utilizadas.)