SqDRIFT algoritmo para la estimación del estado fundamental

Estimación de tiempo de ejecución: 180 segundos en un procesador Heron r3 (NOTA: Se trata únicamente de una estimación. (El tiempo de ejecución puede variar.)

¿Buscas la versión en C++?

En este tutorial se utiliza Python. Para consultar la implementación en C++, incluido el código fuente y las instrucciones de compilación, consulta el tutorial « SqDRIFT » en C++.

Resultados del aprendizaje

Descubre cómo crear circuitos con menor profundidad en comparación con la técnica de Trotter
Recorre un flujo de trabajo completo para la estimación del estado fundamental utilizando « qDRIFT » y SQD
Descubre cómo utilizarlo qiskit-fermions junto con otros complementos de Qiskit para implementar un flujo de trabajo de este tipo

Este tutorial se presenta en formato de cuaderno de « Python » con fines didácticos.

Requisitos previos

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

En segundo plano

SqDRIFT es una variante del SKQD que sustituye la necesidad de elegir un ansatz a partir del cual muestrear cadenas de bits por un conjunto de circuitos de evolución temporal construidos directamente a partir del hamiltoniano objetivo. Esto se consigue mediante el submuestreo de operadores de evolución temporal más pequeños a partir del hamiltoniano, basándose en sus coeficientes, lo que se conoce como el método de « qDRIFT » de Trotterización.

Este tutorial utiliza Qiskit Fermions para crear circuitos fermiónicos más naturales para el algoritmo « qDRIFT », tras lo cual se aplican las fases de diseño y síntesis fermiónicas antes de integrar los circuitos en el proceso tradicional de Qiskit para su ejecución en hardware.

Supongamos que el hamiltoniano tiene la forma:

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

donde, sin pérdida de generalidad, se requiere que $c_i > 0$ y que el mayor valor propio de $h_i$ sea igual, en valor absoluto, a $1$ . Cualquier prefactor con signo o complejo se absorbe en $h_i$ , por lo que los coeficientes $c_i$ son pesos estrictamente positivos, mientras que los $h_i$ determinan 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) del hamiltoniano; se trata de una propiedad del hamiltoniano y es distinto del número de operadores muestreados en un único circuito, que a continuación se denota c $n$ .

A continuación, el algoritmo « qDRIFT » realiza, para el tiempo objetivo $t$ , un operador $V_k$ , donde $k$ va desde $1 \cdots K$ y representa el circuito $k_{th}$ 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 se calcula sobre los sorteos de $n$ , no sobre todos los términos hamiltonianos de $N$ , y, dado que los términos se extraen con repetición, 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, por lo que cada uno de los pasos « $n$ » evoluciona durante el mismo intervalo de tiempo « $\lambda t / n$ », independientemente del término que se haya extraído. La uniformidad del ángulo de paso es la característica distintiva de qDRIFT: : un coeficiente influye en el resultado en función de la frecuencia con la que se extrae su término, y no de cuánto se gira dicho término. Los índices se extraen de la distribución:

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

Por lo tanto, la serie « $(k_1, \ldots, k_n)$ » es una secuencia aleatoria de índices de términos extraídos de esta distribución. Dado que los valores de « $c_i$ » son positivos y su suma es « $\lambda$ », se trata de una distribución de probabilidad normalizada, y la esperanza del canal resultante a partir de los muestreos aleatorios se aproxima a la evolución bajo « $H$ », con un error que disminuye a medida que aumenta « $n$ ». Ten en cuenta que el error de aproximación depende de $\lambda$ y no del número de términos $N$ .

(En el artículo « SqDRIFT » se escribe el número de términos como « $\mathcal{N}$ » y la longitud de la sucesión como « $N$ »; aquí utilizamos « $N$ » y « $n$ » para diferenciar claramente ambos conceptos.)

Este tutorial muestra cómo generar un conjunto de circuitos aleatorios de este tipo. Una vez creados estos circuitos, de forma similar a como se crea un subespacio de Krylov para diferentes operadores, tomamos muestras de cadenas de bits de varios de esos operadores con distintos parámetros temporales. Esto garantiza una mayor superposición entre los vectores del estado fundamental y las cadenas de bits muestreadas.

Requisitos

Antes de empezar este tutorial, asegúrate de que tienes 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,
)

Ejemplo de simulador

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

Cómo leer y preparar el FCIDump

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

Acerca de este FCIDump. El archivo describe N2_sto_3g 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$ , que es la longitud de enlace de equilibrio experimental. En su encabezado se declaran NORB=10, NELEC=14, y MS2=0: 10 orbitales espaciales (por lo tanto, 20 orbitales de espín y 20 qubits según la representación de Jordan-Wigner), 14 electrones en un singulete 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 » en todo el espacio, no hay orbitales congelados y el espacio de correlación es lo suficientemente pequeño como para poder calcular de forma clásica una energía de referencia FCI exacta a efectos comparativos, tal y como se muestra en la siguiente celda.

Se puede volver a generar un archivo equivalente con el comando « PySCF: »

from pyscf import gto, scf, tools

mol = gto.M(atom="N 0 0 0; N 0 0 1.09", basis="sto-3g", symmetry=False)
mf = scf.RHF(mol).run()
tools.fcidump.from_scf(mf, "N2_sto_3g")

Dado que las integrales dependen de los orbitales SCF convergentes, un archivo regenerado podría diferir del original en cuanto a la fase o el orden de los orbitales; las energías totales no se ven afectadas.

Cómo obtener el archivo. Encuentra el FCIDump en este repositorio de GitHub. Puedes ejecutar la celda que aparece a continuación para importarlo a la ubicación que indica el resto del tutorial.

En primer lugar, utilizamos la función proporcionada cisolver por pyscf para obtener la energía de referencia. Esta es la energía real del estado fundamental de la molécula con la que estamos trabajando. Para ello, primero definiremos norb y nelec, que son el número de orbitales y el número de electrones, respectivamente. A continuación, definimos y h1e h2e, que son las integrales de un electrón y de dos electrones, respectivamente. Todo esto se utilizará más adelante también para el SQD.

import os
from urllib.request import urlopen

# The FCIDump is stored with this tutorial in the Qiskit documentation repository.
FCIDUMP_URL = "https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g"
FCIDUMP_PATH = "assets/sqdrift/fcidump_files/N2_sto_3g"

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

Output:

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

Output:

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

Carga del hamiltoniano

Una vez que disponemos de los datos necesarios, 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 sobre fermiones con qiskit-fermions

En primer lugar, representaremos el hamiltoniano en un modelo de circuito fermiónico utilizando qiskit-fermions, que proporciona pasadas de transpilador y puertas específicas para circuitos fermiónicos. Estos se utilizarán posteriormente, antes de que se ejecuten las etapas tradicionales del transpilador de Qiskit para este flujo de trabajo.

Agrupación de términos

Para garantizar la reproducibilidad de los resultados, primero utilizamos canonical_order para ordenar los términos basándonos únicamente en su estructura. Por lo tanto, el orden de los operadores en la lista canon es fijo. Esto garantiza la reproducibilidad de los operadores creados, ya que la función QDriftTrotterization «pass» que utilizaremos en las próximas ejemplos selecciona índices aleatorios para crear los operadores « qDRIFT ».

En este paso, aprovechamos las numerosas simetrías presentes en el hamiltoniano de la estructura electrónica agrupando los términos relacionados que tienen coeficientes idénticos. Aunque al hacerlo se modifica la distribución de los coeficientes de los operadores de la que toma muestras el protocolo « qDRIFT », esto no afecta a sus garantías de convergencia. Es fundamental señalar que agrupar los términos relacionados por simetría da lugar a una cancelación favorable de los términos de Pauli y a una profundidad de circuito globalmente menor al hacer evolucionar un estado en el tiempo bajo su acción.

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

Ten en cuenta que se group_terms_by_electronic_structure supone que los términos siguen un orden normal.

Filtrado de términos diagonales

Eliminamos los términos diagonales del hamiltoniano utilizado para generar los circuitos, de modo que $n$ qDRIFT las ranuras de muestreo se dediquen a los términos que transfieren la población entre configuraciones. Lo mejor es eliminar esos términos del hamiltoniano en este momento, antes de construir la puerta Evolution en el siguiente paso.

Los términos en cuestión son aquellos que son diagonales en la base de números de ocupación, es decir, los productos de los operadores numéricos $a^\dagger_i a_i$ . Hay tres tipos de términos que se ajustan a esta descripción:

el desplazamiento energético constante, un producto de operadores de número cero, cuya evolución temporal solo aporta una fase global;
los operadores de número individual $n_i$ , cuya evolución temporal se reduce a rotaciones de un solo qubit $Z$ ;
los productos de orden superior, como $n_i n_j$ .

Por sí mismas, ninguna de ellas provoca cambios en la distribución de la población entre las configuraciones de ocupación y número; solo actúan sobre las fases de las configuraciones ya existentes. Sin embargo, no son inertes: esas fases relativas influyen en la interferencia generada por los términos de excitación más adelante en el circuito, por lo que filtrarlas modifica la evolución que se genera realmente y puede alterar la distribución de muestreo. Se trata de una aproximación deliberada en la etapa de generación del circuito, realizada para centrar el muestreo en los términos de excitación, y no de una etapa que deje intacta la distribución muestreada. A diferencia de la agrupación de simetría anterior, que mantiene intactas las garantías de convergencia de « qDRIFT », este filtro modifica el operador que se somete a la evolución. Por lo tanto, los circuitos ya no se aproximan a 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 en este caso porque los circuitos no son más que una heurística de muestreo utilizada para proponer configuraciones: no se pierde ningún término de la estimación de la energía en sí, ya que el filtro solo se aplica 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 del 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 sin necesidad de reescribirlo. Los identifica a partir de su estructura de 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álido para un operador que ya tenga 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))

Output:

Ahora que hemos agrupado los términos del hamiltoniano, vamos a fijar los siguientes parámetros para generar el conjunto de circuitos:

El número de circuitos que hay que generar: num_circuits
La longitud de cada circuito en términos de grupos de excitación: num_exc
El factor que explica las diferentes duraciones de la evolución: times

Creación de circuitos fermiónicos

Ahora crearemos circuitos fermiónicos para cada uno de los intervalos de tiempo. Cada circuito estará formado por una única puerta de evolución, con el tiempo de evolución que hemos establecido anteriormente. El operador de evolución es el hamiltoniano. Posteriormente, aplicamos varias pasadas del transpilador a estos circuitos para crear circuitos de tipo « qDRIFT ».

Preparación de Ansatz

Preparamos el estado de Hartree-Fock utilizando la clase InitializeModes . En el caso del nitrógeno, el proceso consiste simplemente en aplicar X puertas a los primeros num_elec_a qubits y, a continuación, a los qubits num_elec_b ; en ambos casos, el número es igual a siete para el nitrógeno. Este estado representa los siete electrones $\alpha$ y los 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 su ejecución en hardware cuántico

Ahora que ya tenemos nuestros circuitos, lo primero que haremos será utilizar las etapas disponibles en qiskit-fermions para llevar a cabo optimizaciones a nivel fermiónico, y a continuación compilaremos nuestro circuito para el backend que elijamos. Dado que se trata de un experimento en simulador, lo haremos primero con el « AerSimulator ».

Cálculo del peso de cada grupo

En este paso, realizamos el muestreo « qDRIFT » de los términos de forma estocástica, con probabilidades proporcionales a sus coeficientes en el hamiltoniano. La etapa del transpilador « qDRIFT » se encarga de ello por nosotros. Ahora podemos crear circuitos menos profundos que se pueden ejecutar en el hardware de forma 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 agrupar los términos, realiza un muestreo de los operadores en función de sus ponderaciones. 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 devuelve generate_preset_jw_pass_manager() un objeto MultiStagePassManager que toma un y FermionicCircuit genera un circuito final optimizado que podemos transpilar para ejecutarlo en nuestro hardware. Sustituimos su etapa de optimización predeterminada por una que FermionicPassManager contenga nuestra pasada QDriftTrotterization :

El paso QDriftTrotterization utiliza internamente el cálculo de pesos y el muestreo para generar los circuitos que utilizaremos para el muestreo
Esta pasada RelabelModes es otra pasada de optimización que se puede utilizar para permutar los modos fermiónicos con el fin de optimizar la conectividad entre los qubits y reducir la profundidad de las puertas; para más información, consulta la referencia de la API

Las etapas restantes se ejecutan MultiStagePassManager automáticamente y se encargan de la asignación completa de fermiones a qubits:

F2QLayout : El gestor de pasadas predefinidas aplica la pasada TrivialF2QLayout , que asigna de forma trivial los bits fermiónicos de $n$ a los qubits de $n$ .
F2QSynth : Una fase de transpilación para convertir las instrucciones de circuitos basados 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))

Output:

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 utilizando Qiskit primitives

Ahora que ya tenemos nuestros circuitos, podemos ejecutarlos utilizando Qiskit primitives en AerSimulator. Sumaremos todos los recuentos de los distintos circuitos. Los convertimos en vectores booleanos antes de someterlos finalmente a un 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")

Output:

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

Paso 4: Realizar el posprocesamiento y obtener el resultado en el formato clásico deseado

Uso de cadenas de bits para SQD

Ahora podemos aplicar el método de diagonalización a las cadenas de bits seleccionadas para hallar el valor propio más bajo, que corresponderá a la energía del estado fundamental de la molécula. Creamos una función de devolución de llamada, declaramos las ocupaciones iniciales y configuramos los parámetros antes de ejecutar finalmente el esquema de diagonalización. La función de llamada de retorno se utiliza para mostrar la iteración actual y la estimación del valor propio actual en cada iteración.

Por último, para obtener la estimación del estado fundamental, sumamos el 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 genera un conjunto diferente de configuraciones, y la etapa de recuperación reestructura el conjunto entre iteraciones, por lo que la dimensión indicada varía de una submuestra a otra. El muestreo silencioso no determina por sí solo la dimensión del subespacio seleccionado. Sin embargo, la ejecución en hardware tiende a generar subespacios sistemáticamente más grandes, ya que las tomas con ruido rompen la simetría del número de partículas y la recuperación de la configuración las convierte en vectores base adicionales. Por eso, también introduciremos otro paso para recortar cadenas de bits en la sección dedicada al hardware.

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

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

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

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

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

result_history = []


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


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

computed_energy = result.energy + nuclear_repulsion_energy

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


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

Output:

40000
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64767025226178
		Subspace dimension: 5538
	Subsample 1
		Energy: -107.64772799119115
		Subspace dimension: 5670
	Subsample 2
		Energy: -107.64765512281548
		Subspace dimension: 5767
Iteration 2
	Subsample 0
		Energy: -107.64795948524682
		Subspace dimension: 6080
	Subsample 1
		Energy: -107.64806617355072
		Subspace dimension: 6300
	Subsample 2
		Energy: -107.64802260640258
		Subspace dimension: 6308
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712
 0.99301927 0.0373282  0.0373266  0.00944508]
Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737
 0.99302145 0.03733536 0.03733399 0.0094585 ]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6480661736 Ha
Error:            1.1811817564e-04 Ha

Ejemplo de hardware

En este ejemplo se utilizan 20 qubits (10 orbitales espaciales). Esa elección es una medida de 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, por lo que lo que determina el coste clásico es la dimensión de ese subespacio seleccionado —que en este caso viene determinada por samples_per_batch, num_batches, y el número de configuraciones distintas que 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 pequeña parte ajustable del mismo, y controlamos su tamaño directamente. Por consiguiente, el número de qubits y la dificultad clásica pueden variar de forma relativamente independiente: un espacio orbital más amplio, muestreado en un subespacio modesto, puede resultar más económico que un sistema más pequeño diagonalizado sobre uno muy grande.

En la práctica, pues, el tamaño viable del sistema depende de la dimensión del subespacio que se necesite para alcanzar la precisión deseada y de la memoria y los núcleos de que disponga el solucionador de valores propios. Los espacios orbitales más amplios suelen requerir un subespacio mayor para alcanzar la precisión química, y eso es lo que, en última instancia, justifica el uso de recursos distribuidos; véase qiskit-addon-sqd-hpc para ampliar este paso. En lugar de fijar un límite arbitrario, lo más práctico es observar la dimensión del subespacio indicada 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 provocado por el ruido del hardware, el subespacio creado para la diagonalización en la ejecución en hardware será mayor que el que obtenemos al utilizar el simulador. Aunque aumenta la dimensión del subespacio que queremos diagonalizar, el procedimiento nos sigue proporcionando una respuesta precisa gracias a la robustez del SQD frente al ruido.

Eliminación de cadenas espurias

Aquí podemos optar por realizar un paso adicional. Cuando dispongamos de todas las cadenas de bits resultantes de las ejecuciones del circuito, podemos filtrar las cadenas de bits no válidas antes de ejecutar el SQD, o bien continuar sin realizar ninguna poda. Por lo general, es preferible omitir la poda en las ejecuciones de hardware, ya que así se mantienen disponibles los resultados con simetría rota para la recuperación de la configuración, lo que permite repararlos y convertirlos en configuraciones válidas, ampliando así el subespacio en lugar de descartar dichos resultados directamente.

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

Ten en cuenta que la poda es solo una de las varias opciones que determinan el 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 no podada solo resulta informativo si se mantienen fijos todos los demás parámetros; el complemento de C++ aborda este tema con más detalle, ya que utiliza la poselección en lugar de la 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")

Output:

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

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

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

Pró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 utiliza circuitos de evolución temporal en lugar de un enfoque variacional.
Diagonalización cuántica basada en muestras de un hamiltoniano químico : un tutorial sobre cómo construir un circuito Jastrow de clúster unitario local (LUCJ) para la simulación de química cuántica.
El artículo « SqDRIFT » (La evolución de la tecnología de la información y la comunicación), en el que se basa este tutorial. (Ten en cuenta que algunas de las optimizaciones que se tratan en este artículo se encuentran actualmente en fase de desarrollo, y que este tutorial está sujeto a cambios en el futuro en función de la evolución de las bibliotecas utilizadas.)

¿Le ha resultado útil esta página?

Informe de un error, de una errata o solicite contenido en GitHub.