Algoritmo SqDRIFT para la estimación del estado fundamental

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

¿Buscas la versión en C++?

Este tutorial usa Python. Para la implementación en C++, incluyendo el código fuente y las instrucciones de compilación, consulta el tutorial de SqDRIFT en C++.

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

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

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.

Este tutorial hace uso de Qiskit Fermions para crear los circuitos fermiónicos más naturales para el algoritmo qDRIFT, seguido del uso de pases de layout y síntesis fermiónicos antes de conectar los circuitos al pipeline tradicional de Qiskit para la 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 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

# Added by doQumentation — required packages for this notebook
!pip install -q numpy pyscf qiskit qiskit-addon-sqd qiskit-aer qiskit-fermions qiskit-ibm-runtime

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

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

# Qiskit Aer
from qiskit_aer import AerSimulator

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

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

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

Ejemplo de simulador

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

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 la siguiente celda.

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 la celda 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 = "assets/sqdrift/fcidump_files/N2_sto_3g"

if not os.path.exists(FCIDUMP_PATH):
    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)
    with urlopen(FCIDUMP_URL) as response:
        contents = response.read()
    with open(FCIDUMP_PATH, "wb") as f:
        f.write(contents)
    print(f"Downloaded FCIDump to {FCIDUMP_PATH}")
else:
    print(f"Using existing FCIDump at {FCIDUMP_PATH}")

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

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 compuertas 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 contribuye 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 compuerta 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 más superficiales 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$$

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 compuertas; 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; el complemento en 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 = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []
hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

# Calculate weights for sampling (one per group)
qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService(channel="ibm_quantum_platform")

# Select backend (choose based on qubit requirements)
backend = service.least_busy(
    operational=True,
    simulator=False,
    min_num_qubits=2 * norb,
)

print(f"Selected backend: {backend.name} ({backend.num_qubits} qubits)")

# Transpile for hardware
transpiled_circuits = transpile(
    sqdrift_circuits,
    backend=backend,
    optimization_level=3,
    seed_transpiler=42,
)

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-SqDRIFT"]

job = sampler.run(transpiled_circuits, shots=shots)
result = job.result()

# Extract counts from SamplerV2 results
all_counts = [pub_result.data.meas.get_counts() for pub_result in result]

# Set to True to filter out bitstrings that violate electron-number conservation
PRUNE = False

def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )

if PRUNE:
    all_counts_filtered = []
    for counts in all_counts:
        filtered_count = {}
        for key in counts:
            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):
                continue
            elif key not in filtered_count.keys():
                filtered_count[key] = counts[key]
            else:
                filtered_count[key] += counts[key]
        all_counts_filtered.append(filtered_count)
    all_counts = all_counts_filtered

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print("Electron configuration:")
print(f"  Total electrons: {nelec}")
print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

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

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

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

Próximos pasos

Recomendaciones

Si te pareció interesante este trabajo, podrían interesarte los siguientes materiales:

• Diagonalización cuántica de Krylov basada en muestras de un modelo de red fermiónico - 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 Jastrow (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 discutidas en este artículo actualmente están en desarrollo, y este tutorial está sujeto a cambios en el futuro según la evolución de las bibliotecas utilizadas.)