SQD pour l'estimation de l'énergie d'un Hamiltonien de chimie

Dans cette leçon, nous allons appliquer SQD pour estimer l'énergie de l'état fondamental d'une molécule.

En particulier, nous aborderons les sujets suivants en utilisant l'approche en $4$ étapes de Qiskit :

1. Étape 1 : Mapper le problème sur des circuits quantiques et des opérateurs
   • Configurer le Hamiltonien moléculaire pour $N_2$.
   • Expliquer l'ansatz de Jastrow à cluster unitaire local (LUCJ), compatible avec le matériel et inspiré de la chimie [1]
2. Étape 2 : Optimiser pour le matériel cible
   • Optimiser le nombre de portes et la disposition de l'ansatz pour l'exécution matérielle
3. Étape 3 : Exécuter sur le matériel cible
   • Exécuter le circuit optimisé sur un vrai QPU pour générer des échantillons du sous-espace.
4. Étape 4 : Post-traiter les résultats
   • Présenter la boucle de récupération de configuration auto-cohérente [2]
     • Post-traiter l'ensemble complet des échantillons de chaînes de bits, en utilisant la connaissance a priori du nombre de particules et l'occupation orbitale moyenne calculée lors de la dernière itération.
     • Créer de manière probabiliste des lots de sous-échantillons à partir des chaînes de bits récupérées.
     • Projeter et diagonaliser le Hamiltonien moléculaire sur chaque sous-espace échantillonné.
     • Sauvegarder l'énergie minimale de l'état fondamental trouvée parmi tous les lots et mettre à jour l'occupation orbitale moyenne.

Nous utiliserons plusieurs packages logiciels tout au long de la leçon.

• PySCF pour définir la molécule et configurer le Hamiltonien.
• Le package ffsim pour construire l'ansatz LUCJ.
• Qiskit pour transpiler l'ansatz en vue de l'exécution matérielle.
• Qiskit IBM Runtime pour exécuter le circuit sur un QPU et collecter des échantillons.
• Qiskit addon SQD pour la récupération de configuration et l'estimation de l'énergie de l'état fondamental par projection sur sous-espace et diagonalisation matricielle.

1. Mapper le problème sur des circuits quantiques et des opérateurs

Hamiltonien moléculaire

Un Hamiltonien moléculaire prend la forme générique suivante :

$$\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \sum_{ \substack{prqs\\\sigma\tau} } \frac{(pr|qs)}{2} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma}$$

$\hat{a}^\dagger_{p\sigma}$/$\hat{a}_{p\sigma}$ sont les opérateurs de création/annihilation fermioniques associés au $p$-ième élément de la base et au spin $\sigma$. $h_{pr}$ et $(pr|qs)$ sont les intégrales électroniques à un et deux corps. En utilisant pySCF, nous allons définir la molécule et calculer les intégrales à un et deux corps du Hamiltonien pour la base 6-31g.

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

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

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

Dans cette leçon, nous utiliserons la transformation de Jordan-Wigner (JW) pour mapper une fonction d'onde fermionique sur une fonction d'onde de qubits, afin qu'elle puisse être préparée à l'aide d'un circuit quantique. La transformation JW mappe l'espace de Fock des fermions dans M orbitales spatiales sur l'espace de Hilbert de 2M qubits : une orbitale spatiale est divisée en deux orbitales de spin, l'une associée à un électron de spin haut ($\alpha$) et l'autre au spin bas ($\beta$). Une orbitale de spin peut être occupée ou inoccupée. En général, quand on parle du nombre d'orbitales, on utilise le nombre d'orbitales spatiales. Le nombre d'orbitales de spin sera le double. Dans les circuits quantiques, nous représentons chaque orbitale de spin par un qubit. Ainsi, un ensemble de qubits représente les orbitales spin-haut ou $\alpha$, et un autre ensemble représente les orbitales spin-bas ou $\beta$. Par exemple, la molécule $N_2$ avec la base 6-31g a $16$ orbitales spatiales (soit $16$ $\alpha$ + $16$ $\beta$ = $32$ orbitales de spin). Nous aurons donc besoin d'un circuit quantique à $32$ qubits (nous pourrions avoir besoin de qubits ancilla supplémentaires, comme discuté plus tard). Les qubits sont mesurés dans la base computationnelle pour générer des chaînes de bits, qui représentent des configurations électroniques ou des déterminants (de Slater). Tout au long de cette leçon, nous utiliserons indifféremment les termes chaînes de bits, configurations et déterminants. Les chaînes de bits nous indiquent l'occupation des électrons dans les orbitales de spin : un $1$ à une position de bit signifie que l'orbitale de spin correspondante est occupée, tandis qu'un $0$ signifie que l'orbitale de spin est vide. Comme les problèmes de structure électronique conservent le nombre de particules, seul un nombre fixe d'orbitales de spin doit être occupé. La molécule $N_2$ possède $5$ électrons spin-haut ($\alpha$) et $5$ électrons spin-bas ($\beta$). Ainsi, toute chaîne de bits représentant les orbitales $\alpha$ et $\beta$ doit avoir cinq $1\text{s}$ chacune pour la molécule $N_2$.

1.1 Circuit quantique pour la génération d'échantillons : l'ansatz LUCJ

Dans cette leçon, nous utiliserons l'ansatz de Jastrow à cluster unitaire local (LUCJ) \[1\] pour la préparation d'état quantique et l'échantillonnage ultérieur. Nous expliquerons d'abord les différents blocs constitutifs de l'ansatz UCJ complet et les approximations effectuées dans sa version locale. Ensuite, en utilisant le package ffsim, nous construirons l'ansatz LUCJ et l'optimiserons avec le transpileur Qiskit pour l'exécution matérielle.

L'ansatz UCJ a la forme suivante (pour un produit de $L$ couches ou répétitions de l'opérateur UCJ) :

$$|\psi\rangle = \prod_{\mu=1}^{L}{(e^{K^{\mu}} \times {e^{iJ^{\mu}}} \times {e^{-K^{\mu}}})} |\Phi_{0}\rangle$$

où $\vert \Phi_{0} \rangle$ est un état de référence, typiquement pris comme l'état de Hartree-Fock (HF). Comme l'état Hartree-Fock est défini en occupant les orbitales de plus faible indice, la préparation de l'état HF impliquera d'appliquer des portes X pour mettre à un les qubits correspondant aux orbitales occupées. Par exemple, le bloc de préparation de l'état HF pour 4 orbitales spatiales et 2 spins haut et 2 spins bas peut ressembler à ceci : Une seule répétition de l'opérateur UCJ ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ consiste en une évolution de Coulomb diagonale ($e^{iJ^{(\mu)}}$) encadrée par des rotations orbitales ($e^{K^{(\mu)}}$ et $e^{-K^{(\mu)}}$). Les blocs de rotation orbitale agissent sur une seule espèce de spin ($\alpha$ (spin haut)/$\beta$ (spin bas)). Pour chaque espèce d'électrons, la rotation orbitale consiste en une couche de portes $R_{z}$ à un qubit suivie d'une séquence de portes de rotation de Givens à 2 qubits (portes $XX + YY$).

Les portes à 2 qubits agissent sur des orbitales de spin adjacentes (qubits voisins les plus proches), et sont donc implémentables sur les QPU IBM® sans nécessiter de portes SWAP. Le $e^{iJ^{(\mu)}}$, également connu sous le nom d'opérateur de Coulomb diagonal, est composé de trois blocs. Deux d'entre eux agissent sur les mêmes secteurs de spin ($e^{iJ_{\alpha \alpha}^{(\mu)}}$ et $e^{iJ_{\beta \beta}^{(\mu)}}$), et un agit entre les deux secteurs de spin ($e^{iJ_{\alpha \beta}^{(\mu)}}$).

Tous les blocs de $e^{iJ^{(\mu)}}$ sont composés de portes nombre-nombre $U_{nn}(\phi)$ [1]. Une porte $U_{nn}(\phi)$ peut être décomposée en une porte $R_{ZZ}(\frac{\phi}{2})$ suivie de deux portes $Rz(-\frac{\phi}{2})$ à un qubit agissant sur deux qubits séparés.

Les composantes de même spin ($J_{\alpha \alpha}$ et $J_{\beta \beta}$) ont des portes $U_{nn}$ entre toutes les paires de qubits possibles. Cependant, comme les QPU supraconducteurs ont une connectivité restrictive, les qubits doivent être échangés (swappés) pour réaliser des portes entre des qubits non adjacents.

Par exemple, considérons le bloc $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (ou $e^{iJ_{\beta \beta}^{(\mu)}}$) suivant pour $N = 4$ orbitales spatiales. Pour une connectivité linéaire des qubits, les trois dernières portes ne sont pas directement implémentables car elles agissent entre des qubits non adjacents (par exemple, Q0 et Q2 ne sont pas directement connectés). Nous avons donc besoin de portes SWAP pour les rendre adjacents (la figure suivante montre un exemple avec $3$ portes SWAP). Ensuite, le $J_{\alpha \beta}$ implémente des portes entre des orbitales de même indice issues de différents secteurs de spin (par exemple, entre $0\alpha$ et $0\beta$). De même, si les qubits ne sont pas physiquement adjacents sur un QPU, ces portes nécessiteront également des SWAPs. D'après la discussion ci-dessus, l'ansatz UCJ présente quelques obstacles pour l'exécution sur matériel car il nécessite des portes SWAP dues aux interactions entre qubits non adjacents. La variante locale de l'ansatz UCJ, le LUCJ, relève ce défi en supprimant certains $U_{nn}$ de l'opérateur de Coulomb diagonal.

Dans les blocs de même espèce d'électrons, $J_{\alpha \alpha}$ et $J_{\beta \beta}$), nous ne conservons que les portes $U_{nn}$ compatibles avec la connectivité aux plus proches voisins et supprimons les portes entre qubits non adjacents dans la version LUCJ. La figure suivante montre le bloc LUCJ après suppression des portes non adjacentes. Ensuite, la version LUCJ du bloc $J_{\alpha \beta}$ qui agit entre différentes espèces d'électrons peut prendre différentes formes selon la topologie du dispositif.

Ici aussi, la version LUCJ se débarrasse des portes non compatibles. La figure ci-dessous montre des variantes du bloc $J_{\alpha \beta}$ pour différentes topologies de qubits, notamment grille, hexagonale, heavy-hex et linéaire.

• Carré : nous pouvons avoir des portes $U_{nn}$ entre toutes les orbitales $\alpha$ et $\beta$ sans aucun SWAP, et par conséquent, nous n'avons pas besoin de supprimer de portes $U_{nn}$.
• Heavy-hex : les interactions $\alpha$-$\beta$ sont maintenues entre les orbitales de spin indexées tous les $4$ (comme les 0e, 4e et 8e) et sont médiées par des qubits ancilla, c'est-à-dire que nous avons besoin de qubits ancilla entre les chaînes linéaires représentant les orbitales $\alpha$ et $\beta$. Cette disposition nécessite un nombre limité de SWAPs.
• Hexagonale : une orbitale sur deux, comme les orbitales d'indice 0, 2 et 4, devient voisine la plus proche lorsque $\alpha$ et $\beta$ sont disposés en deux chaînes linéaires adjacentes.
• Linéaire : une seule orbitale $\alpha$ et une orbitale $\beta$ sont connectées, ce qui signifie que le bloc $J_{\alpha \beta}$ n'aura qu'une seule porte. Bien que la suppression de portes de l'ansatz UCJ pour construire la version LUCJ le rende plus compatible avec le matériel, l'ansatz perd en expressivité. Par conséquent, davantage de répétitions ($L$) de l'opérateur UCJ modifié peuvent être nécessaires lors de l'utilisation de l'ansatz LUCJ.

1.2 Initialisation de l'ansatz LUCJ

Le LUCJ est un ansatz paramétré, et nous devons initialiser les paramètres avant l'exécution matérielle. Une façon d'initialiser l'ansatz est d'utiliser les amplitudes t1 et t2 de la méthode classique des clusters couplés simple et double (CCSD), où les amplitudes t1 sont les coefficients des opérateurs d'excitation simple et les amplitudes t2 sont ceux des opérateurs d'excitation double.

Note que même si l'initialisation de l'ansatz LUCJ avec les amplitudes t1 et t2 produit des résultats corrects, les paramètres de l'ansatz peuvent nécessiter une optimisation supplémentaire.

# 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

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 Construction de l'ansatz LUCJ avec ffsim

Nous utiliserons le package ffsim pour créer et initialiser l'ansatz avec les amplitudes t1 et t2 calculées ci-dessus. Comme notre molécule a un état de Hartree-Fock à couche fermée, nous utiliserons la variante équilibrée en spin de l'ansatz UCJ, UCJOpSpinBalanced.

Comme le matériel IBM a une topologie heavy-hex, nous adopterons le motif zig-zag utilisé dans [1] et expliqué ci-dessus pour les interactions entre qubits. Dans ce motif, les orbitales (qubits) de même spin sont connectées avec une topologie linéaire (cercles rouges et bleus). En raison de la topologie heavy-hex, les orbitales de spin différent ont des connexions entre chaque 4e orbitale, c'est-à-dire les 0e, 4e, 8e, etc. (cercles violets).

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)

L'ansatz LUCJ avec des couches répétées peut être optimisé en fusionnant certains blocs adjacents. Considérons le cas pour n_reps=2. Les deux blocs de rotation orbitale au milieu peuvent être fusionnés en un seul bloc de rotation orbitale. Le package ffsim dispose d'un gestionnaire de passes nommé ffsim.qiskit.PRE_INIT pour optimiser le circuit en fusionnant ces blocs adjacents.

2. Optimiser pour le matériel cible

Tout d'abord, nous récupérons un backend de notre choix. Nous optimiserons notre circuit pour ce backend, puis exécuterons le circuit optimisé sur le même backend pour générer des échantillons du sous-espace.

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

Ensuite, nous recommandons les étapes suivantes pour optimiser l'ansatz et le rendre compatible avec le matériel.

• Sélectionner les qubits physiques (initial_layout) du matériel cible qui respectent le motif zig-zag (deux chaînes linéaires avec un qubit ancilla entre elles) décrit ci-dessus. Disposer les qubits selon ce motif conduit à un circuit efficace et compatible avec le matériel, avec moins de portes.
• Générer un gestionnaire de passes par étapes à l'aide de la fonction generate_preset_pass_manager de Qiskit avec le backend et le initial_layout de ton choix.
• Définir l'étape pre_init de ton gestionnaire de passes par étapes sur ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT inclut des passes du transpileur Qiskit qui décomposent les portes en rotations orbitales puis fusionnent les rotations orbitales, ce qui résulte en moins de portes dans le circuit final.
• Exécuter le gestionnaire de passes sur ton circuit.

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()}")

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. Exécuter sur le matériel cible

Après avoir optimisé le circuit pour l'exécution matérielle, nous sommes prêts à l'exécuter sur le matériel cible et à collecter des échantillons pour l'estimation de l'énergie de l'état fondamental. Comme nous n'avons qu'un seul circuit, nous utiliserons le mode d'exécution par Job de Qiskit Runtime et exécuterons notre circuit.

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. Post-traiter les résultats

La partie post-traitement du flux de travail SQD peut être résumée par le diagramme suivant.

L'échantillonnage de l'ansatz LUCJ dans la base computationnelle génère un ensemble de configurations bruitées $\tilde{\mathcal{\chi}}$, qui sont utilisées dans la routine de post-traitement. Celle-ci implique une méthode appelée récupération de configuration (détails discutés plus loin) pour corriger de manière probabiliste les configurations avec un nombre d'électrons incorrect. Les configurations avec uniquement le nombre correct d'électrons $\tilde{\mathcal{\chi}}_{R}$ sont ensuite sous-échantillonnées et réparties en plusieurs lots en fonction de la fréquence d'apparition de chaque configuration unique. Chaque lot d'échantillons définit un sous-espace ($\mathcal{S^{(k)}}$). Ensuite, le Hamiltonien moléculaire, $H$, est projeté sur les sous-espaces :

$$H_{\mathcal{S}^{(k)}} = P_{\mathcal{S}^{(k)}} H _{\mathcal{S}^{(k)}} \text{ with } P_{\mathcal{S}^{(k)}} = \sum_{x \in \mathcal{S}^{(k)}} \vert x \rangle \langle x \vert$$

Chaque Hamiltonien projeté $H_{\mathcal{S}^{(k)}}$ est ensuite transmis à un solveur propre (Eigensolver), où il est diagonalisé pour calculer les valeurs propres et les vecteurs propres afin de reconstruire un état propre. Dans cette leçon, nous projetons et diagonalisons le Hamiltonien à l'aide du package qiskit-addon-sqd, qui utilise la méthode de Davidson de PySCF pour la diagonalisation.

$$H_{\mathcal{S}^{(k)}} \vert \psi^{(k)} \rangle = E^{(k)} \vert \psi^{(k)} \rangle$$

Nous collectons ensuite la valeur propre (énergie) la plus basse parmi les lots, et calculons également l'occupation orbitale moyenne, $\text{n}$. L'information d'occupation moyenne est utilisée dans l'étape de récupération de configuration pour corriger de manière probabiliste les configurations bruitées.

Ensuite, nous expliquons en détail la boucle de récupération de configuration auto-cohérente et montrons des exemples de code concrets pour implémenter les étapes mentionnées ci-dessus afin d'estimer l'énergie de l'état fondamental du Hamiltonien $N_2$.

4.1 Récupération de configuration : aperçu

Chaque bit d'une chaîne de bits (déterminant de Slater) représente une orbitale de spin. La moitié droite d'une chaîne de bits représente les orbitales spin-haut, et la moitié gauche représente les orbitales spin-bas. Un 1 signifie que l'orbitale est occupée par un électron, et un 0 signifie que l'orbitale est vide. Nous connaissons a priori le nombre correct de particules (électrons spin-haut et spin-bas). Supposons que nous ayons un déterminant $x$ avec $N_x$ électrons (c'est-à-dire qu'il y a $N_x$ nombres de $1$s dans la chaîne de bits). Le nombre correct de particules est $N$. Si $N_x \neq N$, alors nous savons que la chaîne de bits est corrompue par le bruit. La routine de configuration auto-cohérente tente de corriger la chaîne de bits en retournant de manière probabiliste $|N_x - N|$ bits en exploitant les informations d'occupation orbitale moyenne. L'occupation orbitale moyenne ($n$) nous indique la probabilité qu'une orbitale soit occupée par un électron. Si $N_x < N$, nous avons moins d'électrons et devons retourner certains $0$s en $1$s et vice versa.

La probabilité de retournement peut être $|x[i] - avg\_occupancy[i]|$ pour la $i$-ème orbitale de spin. Dans [2], les auteurs ont utilisé une probabilité de retournement pondérée à l'aide d'une fonction ReLU modifiée.

$$\begin{align} w(y) = \begin{cases} \delta \frac{y}{h} & \text{if } y \leq h\\ \nonumber \delta + (1 - \delta) \frac{y - h}{1 - h} & \text{if } y > h \end{cases} \end{align}$$

Ici $h$ définit l'emplacement du « coin » de la fonction ReLU, et le paramètre $\delta$ définit la valeur de la fonction ReLU au coin. Pour $\delta = 0$, $w$ devient une vraie fonction ReLU, et pour $\delta >0$, elle devient une ReLU modifiée. Dans l'article, les auteurs ont utilisé $\delta = 0.01$ et $h =$ nombre de particules alpha (ou beta) / nombre d'orbitales de spin alpha (ou beta) $= N/M$ (facteur de remplissage).

L'occupation orbitale moyenne ($n$) n'est pas connue a priori. La première itération de l'estimation de l'état fondamental commence avec les configurations ayant uniquement le nombre correct de particules dans les deux espèces de spin. Après la première itération, nous avons une estimation de l'état fondamental, et à partir de cette estimation, nous pouvons construire la première approximation de $n$. Cette approximation de $n$ est utilisée pour récupérer des configurations, exécuter l'itération suivante d'estimation de l'état fondamental, et affiner l'approximation de $n$ de manière auto-cohérente. Le processus se répète jusqu'à ce qu'un critère d'arrêt soit atteint.

Considérons l'exemple suivant pour $N = 2$ et $x = |1000\rangle$ ($N_x = 1$). Nous devons retourner l'un des 0s en 1 pour corriger le nombre de particules, et les choix sont 1100, 1010 et 1001. En fonction de la probabilité de retournement, l'un des choix sera sélectionné comme configuration récupérée (ou la chaîne de bits avec le nombre correct de particules).

Supposons que lors de la première itération nous exécutions deux lots, et que les états fondamentaux estimés soient :

$$\begin{align}\nonumber \text{Batch0: } \vert \psi \rangle &= 0.8 \times \vert 1001 \rangle + 0.6 \times \vert 0110 \rangle \\ \nonumber \text{Batch1: } \vert \psi \rangle &= \frac{1}{\sqrt{3}} \left( \vert 1001 \rangle + \vert 0101 \rangle + \vert 0110 \rangle \right) \nonumber \end{align}$$

En utilisant les états de la base computationnelle et leurs amplitudes, nous pouvons calculer la probabilité d'occupation des électrons (en abrégé occupancies) par orbitale de spin (qubit) (noter que la probabilité = |amplitude|$^2$). Ci-dessous, nous tabulons les occupations par qubit pour chaque chaîne de bits apparaissant dans l'état fondamental estimé et calculons l'occupation orbitale totale pour un lot. Notez que, selon la convention d'ordre de Qiskit, le bit le plus à droite représente qubit-0 (Q0), et le bit le plus à gauche représente Q3.

Occupation (Batch0) :

	Q3	Q2	Q1	Q0
1001	0.64	0.0	0.0	0.64
0110	0.0	0.36	0.36	0.0
n (Batch0)	0.64	0.36	0.36	0.64

Occupation (Batch1)

	Q3	Q2	Q1	Q0
1001	0.33	0.00	0.00	0.33
0101	0.0	0.33	0.00	0.33
0110	0.0	0.33	0.33	0.00
n (Batch1)	0.33	0.66	0.33	0.66

Occupation (moyenne sur les lots)

	Q3	Q2	Q1	Q0
n (Batch0)	0.64	0.36	0.36	0.64
n (Batch1)	0.33	0.66	0.33	0.66
n (average)	0.49	0.51	0.35	0.65

En utilisant l'occupation orbitale moyenne calculée ci-dessus, nous pouvons trouver les probabilités de retournement pour différentes orbitales dans la configuration $x = \vert 1000 \rangle$. Comme l'orbitale représentée par Q3 est déjà occupée et n'a pas besoin d'être retournée, nous fixons sa p(flip) à $0$. Pour les orbitales restantes, qui sont inoccupées, la probabilité de retournement est $\vert x[i] - \text{n}[i] \vert$ chacune. En plus de p(flip), nous calculons également le poids de probabilité associé au retournement en utilisant la fonction ReLU modifiée décrite ci-dessus.

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

	Q3	Q2	Q1	Q0
p(flip) ($\vert x[i] - \text{n}[i] \vert$)	0	0.51	0.35	0.65
w(p(flip))	0	0.03	0.007	0.31

Finalement, en utilisant les probabilités pondérées ci-dessus, nous pouvons retourner l'une des orbitales inoccupées Q2, Q1 et Q0. D'après les valeurs ci-dessus, Q0 sera retourné avec la plus grande probabilité, et une configuration récupérée possible peut être $\vert \text{1001} \rangle$. Le processus complet de récupération de configuration auto-cohérente peut être résumé comme suit :

Première itération : Supposons que les chaînes de bits (configurations ou déterminants de Slater) générées par l'ordinateur quantique forment un ensemble $\widetilde{\chi}$, qui inclut à la fois les configurations avec le nombre correct ($\widetilde{\chi}_{correct}$) et incorrect ($\widetilde{\chi}_{incorrect}$) de particules dans chaque secteur de spin.

1. Les configurations de ($\widetilde{\chi}_{correct}$) sont échantillonnées de manière aléatoire pour créer des lots $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ de vecteurs pour la projection sur sous-espace. Le nombre de lots et d'échantillons dans chaque lot sont des paramètres définis par l'utilisateur. Plus le nombre d'échantillons dans chaque lot est grand, plus la dimension du sous-espace est grande et plus la diagonalisation est coûteuse en calcul. D'autre part, un nombre d'échantillons trop faible peut manquer les vecteurs de support de l'état fondamental et conduire à une estimation incorrecte.
2. Exécuter le solveur d'états propres (c'est-à-dire la projection sur sous-espace et la diagonalisation) sur les lots et obtenir des états propres approximatifs $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. À partir des états propres approximatifs, construire la première approximation de $n$.

Itérations suivantes :

1. En utilisant $n$, corriger les configurations avec le mauvais nombre de particules dans $\widetilde{\chi}_{incorrect}$. Nommons-les $\widetilde{\chi}_{correct\_new}$. Alors, $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ forme le nouvel ensemble de configurations avec le nombre correct de particules.
2. $\widetilde{\chi}_{R}$ est échantillonné pour créer des lots $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. Le solveur d'états propres s'exécute avec les nouveaux lots et génère de nouvelles estimations des états fondamentaux $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. À partir des états propres approximatifs, construire une approximation affinée de $n$.
5. Si le critère d'arrêt n'est pas atteint, retourner à l'étape 2.1.

4.2 Estimation de l'état fondamental

Tout d'abord, nous allons transformer les comptages en une matrice de chaînes de bits et un tableau de probabilités pour le post-traitement.

Chaque ligne de la matrice représente une chaîne de bits unique. Comme les qubits sont indexés à partir de la droite d'une chaîne de bits dans Qiskit, la colonne 0 représente le qubit N-1, et la colonne N-1 représente le qubit 0, où N est le nombre de qubits.

Les orbitales alpha sont représentées dans la plage d'index de colonne (N, N/2] (moitié droite), et les orbitales beta sont représentées dans la plage de colonnes (N/2, 0] (moitié gauche).

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)

Il existe quelques options contrôlables par l'utilisateur qui sont importantes pour cette technique :

• iterations : Nombre d'itérations de récupération de configuration auto-cohérente
• n_batches : Nombre de lots de configurations utilisés par les différents appels au solveur d'états propres
• samples_per_batch : Nombre de configurations uniques à inclure dans chaque lot
• max_davidson_cycles : Nombre maximum de cycles de Davidson exécutés par chaque solveur propre

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_davidson=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)

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 Discussion des résultats

Le premier graphique montre qu'après quelques itérations, nous estimons l'énergie de l'état fondamental avec une erreur d'environ 24 mH (la précision chimique est généralement acceptée à 1 kcal/mol $\approx$ 1,6 mH). Le second graphique montre l'occupation moyenne de chaque orbitale spatiale après la dernière itération. On peut voir que les électrons spin-haut et spin-bas occupent tous deux les cinq premières orbitales avec une forte probabilité dans nos solutions.

Bien que l'énergie de l'état fondamental estimée soit correcte, elle n'est pas dans la limite de précision chimique ($\pm \approx 1,6$ mH). Cet écart peut être attribué à la faible dimension du sous-espace utilisée ci-dessus pour la projection et la diagonalisation. Comme nous avons utilisé samples_per_batch=500, le sous-espace est engendré par au maximum $500$ vecteurs, ce qui manque des vecteurs du support de l'état fondamental. Augmenter le paramètre samples_per_batch devrait améliorer la précision au détriment de plus de ressources de calcul classique et de temps d'exécution.

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

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

Exercice pour le lecteur

Augmente progressivement le paramètre samples_per_batch (par exemple, de $1000$ à $10000$ par pas de $1000$ ; dans la limite de la mémoire de ton ordinateur) et compare les énergies d'état fondamental estimées.

Références

[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). Chem. Sci., 2023, 14, 11213.

[2] J. Robledo-Moreno et al., "Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer" (2024). arXiv:quant-ph/2405.05068.