SQD para la estimación energética de un hamiltoniano químico

En esta lección, aplicaremos SQD para estimar la energía de estado fundamental de una molécula.

En concreto, trataremos los siguientes temas utilizando el enfoque del patrón Qiskit $4$ -step:

1. Paso 1: Asignar el problema a circuitos y operadores cuánticos
   • Configurar el Hamiltoniano molecular para $N_2$.
   • Explicar el clúster unitario local Jastrow (LUCJ), inspirado en la química y fácil de usar en hardware [1]
2. Paso 2: Optimización para el hardware de destino
   • Optimizar el número de puertas y el diseño del ansatz para su ejecución en hardware
3. Paso 3: Ejecutar en el hardware de destino
   • Ejecuta el circuito optimizado en una QPU real para generar muestras del subespacio.
4. Paso 4: Tratamiento posterior de los resultados
   • Introducir el bucle de recuperación de configuración autoconsistente [2]
     • Post-procesar el conjunto completo de muestras de cadenas de bits, utilizando el conocimiento previo del número de partículas y la ocupación orbital media calculada en la iteración más reciente.
     • Crear probabilísticamente lotes de submuestras a partir de las cadenas de bits recuperadas.
     • Proyectar y diagonalizar el Hamiltoniano molecular sobre cada subespacio muestreado.
     • Guarda la energía mínima del estado base encontrada en todos los lotes y actualiza la ocupación orbital media.

Utilizaremos varios paquetes de software a lo largo de la lección.

• PySCF para definir la molécula y configurar el Hamiltoniano.
• ffsim para construir el ansatz LUCJ.
• Qiskit para transpilar el ansatz para su ejecución en hardware.
• Qiskit IBM Runtime para ejecutar el circuito en una QPU y recoger muestras.
• Qiskit addon SQD recuperación de la configuración y estimación de la energía del estado base mediante proyección de subespacios y diagonalización de matrices.

1. Asignar el problema a circuitos y operadores cuánticos

Hamiltoniano molecular

Un Hamiltoniano molecular toma la forma genérica:

$\hat{a}^\dagger_{p\sigma}$ / $\hat{a}_{p\sigma}$ son los operadores fermiónicos de creación/aniquilación asociados al $p$ -ésimo elemento del conjunto de bases y al espín $\sigma$. $h_{pr}$ y $(pr|qs)$ son las integrales electrónicas de uno y dos cuerpos. Utilizando pySCF, definiremos la molécula y calcularemos las integrales de uno y dos cuerpos del Hamiltoniano para el conjunto de bases 6-31g.

import warnings
import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# Specify molecule properties
open_shell = False
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],  # Two N atoms 1 angstrom apart
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)  # hcore: one-body integrals
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)  # eri: two-body integrals

# Compute exact energy for comparison
exact_energy = cas.run().e_tot

Output:

converged SCF energy = -108.835236570774
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

En esta lección, utilizaremos la transformación Jordan-Wigner (JW) para mapear una función de onda fermiónica a una función de onda qubit de forma que pueda ser preparada utilizando un circuito cuántico. La transformación JW mapea el espacio de Fock de fermiones en M orbitales espaciales en el espacio de Hilbert de 2M qubits, es decir, un orbital espacial se divide en dos orbitales de espín, uno asociado con un electrón de espín arriba ( $\alpha$ ) y otro con espín abajo ( $\beta$ ). Un orbital de espín puede estar ocupado o desocupado. Normalmente, cuando nos referimos al número de orbitales, estaremos utilizando el número de orbitales espaciales. El número de orbitales de espín será el doble. En los circuitos cuánticos, representaremos cada orbital de espín con un qubit. Así, un conjunto de qubits representará los orbitales de espín hacia arriba o $\alpha$, y otro conjunto representará los orbitales de espín hacia abajo o $\beta$. Por ejemplo, la molécula $N_2$ para el conjunto de bases 6-31g tiene orbitales espaciales $16$ (es decir, $16$ $\alpha$ + $16$ $\beta$ = $32$ orbitales de espín). Por lo tanto, necesitaremos un circuito cuántico $32$ -qubit (es posible que necesitemos qubits ancilla adicionales, como se discutirá más adelante). Los qubits se miden en base computacional para generar cadenas de bits, que representan configuraciones electrónicas o determinantes (de Slater). A lo largo de esta lección, utilizaremos indistintamente los términos cadenas de bits, configuraciones y determinantes. Las cadenas de bits nos indican la ocupación de electrones en orbitales de espín: un $1$ en una posición de bit significa que el orbital de espín correspondiente está ocupado, mientras que un $0$ significa que el orbital de espín está vacío. Dado que los problemas de estructura electrónica preservan las partículas, sólo debe ocuparse un número fijo de orbitales de espín. La molécula $N_2$ tiene $5$ electrones de espín ascendente ( $\alpha$ ) y $5$ electrones de espín descendente ( $\beta$ ). Así, cualquier cadena de bits que represente los orbitales $\alpha$ y $\beta$ debe tener cinco $1\text{s}$ cada uno para la molécula $N_2$.

1.1 Circuito cuántico para la generación de muestras: el enfoque LUCJ

En esta lección, utilizaremos el ansatz local unitario acoplado de Jastrow (LUCJ) \ [1] para la preparación del estado cuántico y el posterior muestreo. En primer lugar, explicaremos los distintos componentes del ansatz UCJ completo y las aproximaciones realizadas en la versión local del mismo. A continuación, utilizando el paquete ffsim, construiremos el ansatz LUCJ y lo optimizaremos utilizando el transpilador Qiskit para su ejecución en hardware.

El ansatz UCJ tiene la siguiente forma (para un producto de $L$ capas o repeticiones del operador UCJ.)

donde, $\vert \Phi_{0} \rangle$ es un estado de referencia, típicamente tomado como el estado Hartree-Fock (HF). Como el estado Hartree-Fock se define por tener ocupados los orbitales de menor número, la preparación del estado HF implicará aplicar puertas X para poner a uno los qubits correspondientes a los orbitales ocupados. Por ejemplo, el bloque de preparación del estado HF para 4 orbitales espaciales y 2 up- y 2 down-spin puede tener el siguiente aspecto:

Una sola repetición del operador UCJ ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ consiste en una evolución diagonal de Coulomb ( $e^{iJ^{(\mu)}}$ ) intercalada por rotaciones orbitales ( $e^{K^{(\mu)}}$ y $e^{-K^{(\mu)}}$ ).

Los bloques de rotación orbital funcionan con una sola especie de espín ( $\alpha$ (up-spin)/ $\beta$ (down-spin)). Para cada especie de electrón, la rotación orbital consiste en una capa de compuertas de un solo qubit $R_{z}$ seguidas de una secuencia de compuertas de rotación de 2 qubits Given ( $XX + YY$ gates).

Las puertas de 2 qubits actúan sobre orbitales de espín adyacentes (qubits vecinos más cercanos) y, por tanto, pueden implementarse en las QPU de IBM® sin necesidad de puertas SWAP.

Esquema

El $e^{iJ^{(\mu)}}$, también conocido como operador diagonal de Coulomb, consta de tres bloques. Dos de ellos funcionan en los mismos sectores de espín ( $e^{iJ_{\alpha \alpha}^{(\mu)}}$ y $e^{iJ_{\beta \beta}^{(\mu)}}$ ), y uno funciona entre dos sectores de espín ( $e^{iJ_{\alpha \beta}^{(\mu)}}$ ).

Todos los bloques en $e^{iJ^{(\mu)}}$ consisten en puertas número-número $U_{nn}(\phi)$ [1]. Una puerta $U_{nn}(\phi)$ puede descomponerse a su vez en una puerta $R_{ZZ}(\frac{\phi}{2})$ seguida de dos puertas $Rz(-\frac{\phi}{2})$ de un solo qubit que actúan sobre dos qubits separados.

Los componentes del mismo espín ( $J_{\alpha \alpha}$ y $J_{\beta \beta}$ ) tienen $U_{nn}$ puertas entre todos los pares posibles de qubits. Sin embargo, como las QPU superconductoras tienen una conectividad restrictiva, los qubits deben intercambiarse para realizar puertas entre qubits no adyacentes.

Por ejemplo, considere el siguiente bloque $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (o $e^{iJ_{\beta \beta}^{(\mu)}}$ ) para orbitales espaciales $N = 4$. Para una conectividad lineal de qubits, las tres últimas puertas no son directamente implementables, ya que funcionan entre qubits no adyacentes (por ejemplo, Q0 y Q2 no están directamente conectados). Por lo tanto, necesitamos puertas SWAP para hacerlas adyacentes (la siguiente figura muestra un ejemplo con $3$ puertas SWAP).

A continuación, $J_{\alpha \beta}$ implementa puertas entre los mismos orbitales indexados de diferentes sectores de espín (por ejemplo, entre $0\alpha$ y $0\beta$ ). Del mismo modo, si los qubits no son físicamente adyacentes en una QPU, estas puertas también requerirán SWAPs.

A partir de la discusión anterior, el ansatz UCJ se enfrenta a algunos obstáculos para la ejecución HW, ya que necesita puertas SWAP debido a las interacciones qubit no adyacentes. La variante local del ansatz UCJ, LUCJ, aborda este reto eliminando algunos $U_{nn}$ del operador de Coulomb diagonal.

En los mismos bloques de especies de electrones, $J_{\alpha \alpha}$ y $J_{\beta \beta}$ ), sólo mantenemos las puertas $U_{nn}$ compatibles con la conectividad de vecino más cercano y eliminamos las puertas entre qubits no adyacentes en la versión LUCJ. La siguiente figura muestra el bloque LUCJ después de eliminar las puertas no adyacentes.

A continuación, la versión LUCJ del bloque $J_{\alpha \beta}$ que funciona entre diferentes especies de electrones puede adoptar diferentes formas en función de la topología del dispositivo.

También en este caso, la versión LUCJ se deshace de las puertas no compatibles. La siguiente figura muestra variantes del bloque $J_{\alpha \beta}$ para diferentes topologías de qubits, incluyendo rejilla, hexagonal, heavy-hex y lineal.

• Cuadrado : podemos tener puertas $U_{nn}$ entre todos los orbitales $\alpha$ y $\beta$ sin ningún SWAP, y por lo tanto, no necesitamos eliminar ninguna puerta $U_{nn}$.
• Heavy-hex : Las interacciones $\alpha$ - $\beta$ se mantienen entre cada $4$ -ésimo orbital de espín indexado (como el 0º, 4º y 8º) y están mediadas por ancilla, es decir, necesitamos qubits ancilla entre las cadenas lineales que representan los orbitales $\alpha$ y $\beta$. Este acuerdo necesita un número limitado de SWAPs.
• Hexagonal : Todos los demás orbitales, como los orbitales 0, 2 y 4 indexados, se convierten en vecinos más cercanos cuando $\alpha$ y $\beta$ se disponen en dos cadenas lineales adyacentes.
• Lineal : Sólo se conectan un orbital $\alpha$ y otro $\beta$, lo que significa que el bloque $J_{\alpha \beta}$ sólo tendrá una puerta.

Diagramas de

Aunque la eliminación de puertas del ansatz UCJ para construir la versión LUCJ lo hace más compatible con HW, el ansatz pierde algo de expresividad. Por lo tanto, pueden ser necesarias más repeticiones ( $L$ ) del operador UCJ modificado cuando se utiliza el ansatz LUCJ.

1.2 Inicialización del ansatz LUCJ

El LUCJ es un ansatz parametrizado, y necesitamos inicializar los parámetros antes de la ejecución del hardware. Una forma de inicializar el ansatz es utilizando las amplitudes t1 y t2 del método clásico de clústeres simples y dobles acoplados (CCSD), donde las amplitudes t1 son el coeficiente de los operadores de excitación simple y las amplitudes t2 son para los operadores de excitación doble.

Obsérvese que aunque la inicialización del ansatz LUCJ con t1 y t2 amplitudes genera resultados decentes, los parámetros del ansatz pueden necesitar una mayor optimización.

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
)
ccsd.run()

t1 = ccsd.t1
t2 = ccsd.t2

Output:

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 Construcción del enfoque LUCJ utilizando ffsim

Utilizaremos el paquete ffsim para crear e inicializar el ansatz con t1 y t2 amplitudes calculadas anteriormente. Dado que nuestra molécula tiene un estado Hartree-Fock de cáscara cerrada, utilizaremos la variante de espín equilibrado del ansatz UCJ, UCJOpSpinBalanced.

Como el hardware de IBM tiene una topología heavy-hex, adoptaremos el patrón zig-zag utilizado en [1] y explicado anteriormente para las interacciones qubit. En este patrón, los orbitales (qubits) con el mismo espín están conectados con una topología de línea (círculos rojos y azules). Debido a la topología hexagonal pesada, los orbitales para diferentes espines tienen conexiones entre cada 4 orbitales, es decir, el 0, el 4, el 8, etc. (círculos morados).

import ffsim
from qiskit import QuantumCircuit, QuantumRegister

n_reps = 2
alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]
alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]

ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),
)

nelec = (num_elec_a, num_elec_b)

# create an empty quantum circuit
qubits = QuantumRegister(2 * num_orbitals, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()
# circuit.decompose().draw("mpl", scale=0.5, fold=-1)

El ansatz LUCJ con capas repetidas puede optimizarse fusionando algunos bloques adyacentes. Consideremos un caso para n_reps=2. Los dos bloques de rotación orbital del centro pueden fusionarse en un único bloque de rotación orbital. El paquete ffsim dispone de un gestor de pasadas denominado ffsim.qiskit.PRE_INIT para optimizar el circuito fusionando dichos bloques adyacentes.

2. Optimizar para el hardware de destino

En primer lugar, buscamos un backend de nuestra elección. Optimizaremos nuestro circuito para el backend, y luego ejecutaremos el circuito optimizado en el mismo backend para generar muestras para el subespacio.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

A continuación, recomendamos los siguientes pasos para optimizar el ansatz y hacerlo compatible con el hardware.

• Seleccione qubits físicos (initial_layout) del hardware de destino que se adhieran al patrón en zig-zag (dos cadenas lineales con un qubit ancilla entre ellas) descrito anteriormente. La disposición de los qubits según este patrón da lugar a un circuito eficiente compatible con hardware y con menos puertas.
• Genere un gestor de pases por etapas utilizando la función generate_preset_pass_manager de Qiskit con su elección de backend y initial_layout.
• Establezca la etapa pre_init de su gestor de pases escalonados en ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT incluye pases de transpilador Qiskit que descomponen las puertas en rotaciones orbitales y luego fusionan las rotaciones orbitales, lo que resulta en menos puertas en el circuito final.
• Ejecuta el gestor de pases en tu circuito.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]
spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]

initial_layout = spin_a_layout + spin_b_layout

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend, initial_layout=initial_layout
)

# without PRE_INIT passes
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}")

# with PRE_INIT passes
# We will use the circuit generated by this pass manager for hardware execution
pass_manager.pre_init = ffsim.qiskit.PRE_INIT
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}")

Output:

Gate counts (w/o pre-init passes): OrderedDict({'rz': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})

3. Ejecutar en el hardware de destino

Tras optimizar el circuito para su ejecución en hardware, estamos listos para ejecutarlo en el hardware de destino y recopilar muestras para la estimación de la energía del estado fundamental. Como solo tenemos un circuito, utilizaremos el y qiskit-ibm-runtimeModo de ejecución de tareas ejecutaremos nuestro circuito.

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
sampler.options.dynamical_decoupling.enable = True

job = sampler.run([isa_circuit], shots=10_000)  # Takes approximately 5sec of QPU time

# Run cell after IQX job completion
primitive_result = job.result()
pub_result = primitive_result[0]
counts = pub_result.data.meas.get_counts()

4. Resultados posteriores al procesamiento

La parte de posprocesamiento del flujo de trabajo de SQD puede resumirse mediante el siguiente diagrama.

El muestreo del ansatz LUCJ en la base de cálculo genera un conjunto de configuraciones ruidosas $\tilde{\mathcal{\chi}}$, que se utilizan en la rutina de posprocesamiento. Se trata de un método denominado recuperación de configuraciones (que se explicará más adelante) para corregir probabilísticamente las configuraciones con números de electrones incorrectos. A continuación, las configuraciones que sólo tienen los números de electrones correctos $\tilde{\mathcal{\chi}}_{R}$ se submuestrean y se distribuyen en varios lotes en función de la frecuencia de aparición de cada configuración única. Cada lote de muestras define un subespacio ( $\mathcal{S^{(k)}}$ ). A continuación, el Hamiltoniano molecular, $H$, se proyecta en subespacios:

Cada Hamiltoniano proyectado $H_{\mathcal{S}^{(k)}}$ se introduce en un Eigensolver, donde se diagonaliza para calcular los valores propios y los vectores propios para reconstruir un estado propio. En esta lección, proyectamos y diagonalizamos el Hamiltoniano utilizando el paquete qiskit-addon-sqd que utiliza el método de Davidson de PySCF para la diagonalización.

A continuación, recogemos el valor propio más bajo (energía) de los lotes y también calculamos la ocupación orbital media, $\text{n}$. La información sobre la ocupación media se utiliza en el paso de recuperación de la configuración para corregir probabilísticamente las configuraciones con ruido.

A continuación, explicamos en detalle el bucle de recuperación de la configuración autoconsistente y mostramos ejemplos concretos de código para implementar los pasos mencionados con el fin de estimar la energía del estado fundamental del Hamiltoniano de $N_2$.

4.1 Recuperación de la configuración: descripción general

Cada bit de una cadena de bits (determinante de Slater) representa un orbital de espín. La mitad derecha de una cadena de bits representa orbitales de espín ascendente y la mitad izquierda orbitales de espín descendente. Un 1 significa que el orbital está ocupado por un electrón, y un 0 significa que el orbital está vacío. Conocemos a priori el número correcto de partículas (electrón de espín ascendente y electrón de espín descendente). Supongamos que tenemos un determinante $x$ con $N_x$ electrones (es decir, hay $N_x$ números de $1$ s en la cadena de bits) en él. El número correcto de partículas es $N$. Si $N_x \neq N$, entonces sabemos que la cadena de bits está corrompida por ruido. La rutina de configuración autoconsistente intenta corregir la cadena de bits volteando probabilísticamente $|N_x - N|$ bits aprovechando la información de ocupación orbital media. La ocupación orbital media ( $n$ ) nos indica la probabilidad de que un orbital esté ocupado por un electrón. Si $N_x < N$, tenemos menos electrones y necesitamos voltear algunos $0$ s a $1$ s y viceversa.

La probabilidad de flipping puede ser $|x[i] - avg\_occupancy[i]|$ para i-th spin orbital. En [2], los autores utilizaron una probabilidad ponderada de volteo utilizando la función ReLU modificada.

Aquí $h$ define la ubicación de la "esquina" de la función ReLU, y el parámetro $\delta$ define el valor de la función ReLU en la esquina. Para $\delta = 0$, $w$ se convierte en la función ReLU verdadera, y para $\delta >0$, se convierte en ReLU* modificada*. En el artículo, los autores utilizaron $\delta = 0.01$ y $h =$ número de partículas alfa (o beta)/número de orbitales de espín alfa (o beta) $= N/M$ (factor de llenado).

La ocupación orbital media ( $n$ ) no se conoce a priori. La primera iteración de la estimación del estado básico comienza con configuraciones que sólo tienen números de partículas correctos en ambas especies de espín. Después de la primera iteración, tenemos una estimación del estado base, y utilizando la estimación, podemos construir la primera conjetura de $n$. Esta conjetura de $n$ se utiliza para recuperar configuraciones, ejecutar la siguiente iteración de estimación del estado base, y auto-consistentemente refinar la conjetura de $n$. El proceso se repite hasta que se cumple un criterio de parada.

Considere el siguiente ejemplo para $N = 2$ y $x = |1000\rangle$ ( $N_x = 1$ ). Tenemos que dar la vuelta a uno de los 0s a 1 para corregirlo para los números de partículas, y las opciones son 1100, 1010, y 1001. En función de la probabilidad de volteo, se seleccionará una de las opciones como configuración recuperada (o la cadena de bits con el número correcto de partículas).

Supongamos que en la primera iteración ejecutamos dos lotes, y los estados básicos estimados a partir de ellos son:

Utilizando los estados base computacionales y sus amplitudes, podemos calcular la probabilidad de las ocupaciones de electrones (en breve ocupaciones ) por orbital de espín (qubit) (nótese que probabilidad = |amplitud| $^2$ ). A continuación tabulamos las ocupaciones por qubit para cada cadena de bits que aparece en el estado básico estimado y calculamos la ocupación orbital total para un lote. Tenga en cuenta que, según la convención de ordenación de Qiskit, el bit situado más a la derecha representa qubit-0 ( Q0 ), y el situado más a la izquierda representa Q3.

Ocupación ( Batch0 ):

Ocupación ( Batch1 )

Ocupación (media de los lotes)

Utilizando la ocupación orbital media calculada anteriormente, podemos hallar las probabilidades de volteo para diferentes orbitales en la configuración $x = \vert 1000 \rangle$. Como el orbital representado por Q3 ya está ocupado y no es necesario voltearlo, fijamos su p(flip) en $0$. Para el resto de orbitales, que están desocupados, la probabilidad de volteo es $\vert x[i] - \text{n}[i] \vert$ para cada uno. Junto con p(flip), también calculamos el peso de probabilidad asociado a flipping utilizando la función ReLU modificada descrita anteriormente.

Probabilidad de flip ( $x = \vert 1000 \rangle$, $\delta = 0.01$, $h = N/M = 2/4 = 0.50$ )

Finalmente, utilizando las probabilidades ponderadas anteriores, podemos voltear uno de los orbitales desocupados Q2, Q1, y Q0. Basándose en los valores anteriores, lo más probable es que se invierta Q0, y una posible configuración recuperada puede ser $\vert \text{1001} \rangle$.

El proceso completo de recuperación de configuraciones autoconsistentes puede resumirse como sigue:

Primera iteración: Supongamos que las cadenas de bits (configuraciones o determinantes de Slater) generadas por el ordenador cuántico forman un conjunto $\widetilde{\chi}$, que incluye tanto configuraciones con número correcto ( $\widetilde{\chi}_{correct}$ ) como incorrecto ( $\widetilde{\chi}_{incorrect}$ ) de partículas en cada sector de espín.

1. Las configuraciones de ( $\widetilde{\chi}_{correct}$ ) se muestrean aleatoriamente para crear lotes $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ de vectores para la proyección del subespacio. El número de lotes y las muestras de cada lote son parámetros definidos por el usuario. Cuanto mayor sea el número de muestras de cada lote, mayor será la dimensión del subespacio y más exigente será la diagonalización desde el punto de vista informático. Por otra parte, un número demasiado pequeño de muestras puede pasar por alto los vectores de apoyo del estado básico y dar lugar a una estimación incorrecta.
2. Ejecute el solucionador de estados propios (es decir, la proyección sobre el subespacio y la diagonalización) en los lotes y obtenga los estados propios aproximados. $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. A partir de los estados propios aproximados, construya la primera conjetura para $n$.

Iteraciones posteriores:

1. Utilizando $n$ corregimos las configuraciones con número de partícula erróneo en $\widetilde{\chi}_{incorrect}$. Supongamos que las llamamos $\widetilde{\chi}_{correct\_new}$. Entonces, $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ forma el nuevo conjunto de configuraciones con número de partículas correcto.
2. $\widetilde{\chi}_{R}$ se muestrea para crear los lotes $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. El solucionador de estados propios se ejecuta con nuevos lotes y genera nuevas estimaciones de los estados básicos $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. A partir de los estados propios aproximados se construyen conjeturas refinadas para $n$.
5. Si no se cumple el criterio de parada, vuelva al paso 2.1.

4.2 Estimación del estado fundamental

En primer lugar, transformaremos los recuentos en una matriz de cadenas de bits y una matriz de probabilidades para su postprocesamiento.

Cada fila de la matriz representa una única cadena de bits. Dado que los qubits se indexan desde la derecha de una cadena de bits en Qiskit, la columna 0 representa el qubit N-1, y la columna N-1 representa el qubit 0, donde N es el número de qubits.

Los orbitales alfa se representan en el rango de índice de columna (N, N/2] (mitad derecha), y los orbitales beta en el rango de columna (N/2, 0] (mitad izquierda).

from qiskit_addon_sqd.counts import counts_to_arrays

# Convert counts into bitstring and probability arrays
bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)

Hay algunas opciones controladas por el usuario que son importantes para esta técnica:

• iterations: Número de iteraciones de recuperación de configuración autoconsistente
• n_batches: Número de lotes de configuraciones utilizados por las diferentes llamadas al solucionador de estados propios
• samples_per_batch: Número de configuraciones únicas a incluir en cada lote
• max_davidson_cycles: Número máximo de ciclos Davidson ejecutados por cada eigensolver

import numpy as np
from qiskit_addon_sqd.configuration_recovery import recover_configurations
from qiskit_addon_sqd.fermion import (
    bitstring_matrix_to_ci_strs,
    solve_fermion,
)
from qiskit_addon_sqd.subsampling import postselect_and_subsample

rng = np.random.default_rng(24)
# SQD options
iterations = 5

# Eigenstate solver options
n_batches = 5
samples_per_batch = 500
max_davidson_cycles = 300

# Self-consistent configuration recovery loop
e_hist = np.zeros((iterations, n_batches))  # energy history
s_hist = np.zeros((iterations, n_batches))  # spin history
occupancy_hist = []
avg_occupancy = None
for i in range(iterations):
    print(f"Starting configuration recovery iteration {i}")
    # On the first iteration, we have no orbital occupancy information from the
    # solver, so we begin with the full set of noisy configurations.
    if avg_occupancy is None:
        bs_mat_tmp = bitstring_matrix_full
        probs_arr_tmp = probs_arr_full

    # If we have average orbital occupancy information, we use it to refine
    # the full set of noisy configurations.
    else:
        bs_mat_tmp, probs_arr_tmp = recover_configurations(
            bitstring_matrix_full,
            probs_arr_full,
            avg_occupancy,
            num_elec_a,
            num_elec_b,
            rand_seed=rng,
        )

    # Create batches of subsamples. We postselect here to remove configurations
    # with incorrect hamming weight during iteration 0, since no config recovery was performed.
    batches = postselect_and_subsample(
        bs_mat_tmp,
        probs_arr_tmp,
        hamming_right=num_elec_a,
        hamming_left=num_elec_b,
        samples_per_batch=samples_per_batch,
        num_batches=n_batches,
        rand_seed=rng,
    )

    # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.
    e_tmp = np.zeros(n_batches)
    s_tmp = np.zeros(n_batches)
    occs_tmp = []
    coeffs = []
    for j in range(n_batches):
        strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])
        print(f"  Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}")
        energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(
            batches[j],
            hcore,
            eri,
            open_shell=open_shell,
            spin_sq=spin_sq,
            max_cycle=max_davidson_cycles,
        )
        energy_sci += nuclear_repulsion_energy
        e_tmp[j] = energy_sci
        s_tmp[j] = spin
        occs_tmp.append(avg_occs)
        coeffs.append(coeffs_sci)

    # Combine batch results
    avg_occupancy = tuple(np.mean(occs_tmp, axis=0))

    # Track optimization history
    e_hist[i, :] = e_tmp
    s_hist[i, :] = s_tmp
    occupancy_hist.append(avg_occupancy)

Output:

Starting configuration recovery iteration 0
  Batch 0 subspace dimension: 21609
  Batch 1 subspace dimension: 21609
  Batch 2 subspace dimension: 21609
  Batch 3 subspace dimension: 21609
  Batch 4 subspace dimension: 21609
Starting configuration recovery iteration 1
  Batch 0 subspace dimension: 609961
  Batch 1 subspace dimension: 616225
  Batch 2 subspace dimension: 627264
  Batch 3 subspace dimension: 633616
  Batch 4 subspace dimension: 624100
Starting configuration recovery iteration 2
  Batch 0 subspace dimension: 564001
  Batch 1 subspace dimension: 605284
  Batch 2 subspace dimension: 582169
  Batch 3 subspace dimension: 559504
  Batch 4 subspace dimension: 591361
Starting configuration recovery iteration 3
  Batch 0 subspace dimension: 550564
  Batch 1 subspace dimension: 549081
  Batch 2 subspace dimension: 531441
  Batch 3 subspace dimension: 527076
  Batch 4 subspace dimension: 531441
Starting configuration recovery iteration 4
  Batch 0 subspace dimension: 544644
  Batch 1 subspace dimension: 580644
  Batch 2 subspace dimension: 527076
  Batch 3 subspace dimension: 531441
  Batch 4 subspace dimension: 537289

4.3 Discusión de los resultados

El primer gráfico muestra que, tras unas pocas iteraciones, estimamos la energía del estado básico con un margen de ~24 mH (la precisión química suele aceptarse en 1 kcal/mol $\approx$ 1.6 mH ). El segundo gráfico muestra la ocupación media de cada orbital espacial tras la iteración final. Podemos ver que tanto los electrones de espín arriba como los de espín abajo ocupan los cinco primeros orbitales con alta probabilidad en nuestras soluciones.

Aunque la energía estimada del estado básico es decente, no está dentro del límite de precisión química ( $\pm \approx 1.6$ mH ). Este desfase puede atribuirse a la pequeña dimensión del subespacio que utilizamos anteriormente para la proyección y la diagonalización. Como utilizamos samples_per_batch=500, el subespacio está abarcado por el máximo de vectores de $500$, que son los vectores que faltan en el soporte del estado base. Aumentar el parámetro samples_per_batch debería mejorar la precisión a costa de más recursos informáticos clásicos y tiempo de ejecución.

# Data for energies plot
x1 = range(iterations)
min_e = [np.min(e) for e in e_hist]
e_diff = [abs(e - exact_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]
x2 = range(len(y2))

import matplotlib.pyplot as plt

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-6)
axs[0].axhline(
    y=chem_accuracy, color="#BF5700", linestyle="--", label="chemical accuracy"
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

print(f"Exact energy: {exact_energy:.5f} Ha")
print(f"SQD energy: {min_e[-1]:.5f} Ha")
print(f"Absolute error: {e_diff[-1]:.5f} Ha")
plt.tight_layout()
plt.show()

Output:

Exact energy: -109.04667 Ha
SQD energy: -109.02234 Ha
Absolute error: 0.02434 Ha

Ejercicio para el lector

Aumente progresivamente el parámetro samples_per_batch (por ejemplo, de $1000$ a $10000$ a un paso de $1000$; permitido mi memoria de su ordenador) y compare las energías estimadas del estado básico.

Referencias

[1] M. Motta et al., "Bridging physical intuition and hardware efficiency for correlated electronic states: the local unitary cluster Jastrow ansatz for electronic structure" (2023). Química. Sci., 2023, 14, 11213.

[2] J. Robledo-Moreno et al., "Química más allá de las soluciones exactas en un superordenador centrado en la cuántica" (2024). arXiv:quant-ph/2405.05068.