SQD pour l'estimation énergétique d'un hamiltonien chimique

Dans cette leçon, nous appliquerons la méthode SQD pour estimer l'énergie de l'état fondamental d'une molécule.

En particulier, nous aborderons les sujets suivants en utilisant l'approche du modèle Qiskit $4$ -step :

1. Étape 1 : Tracer le problème à l'aide de circuits et d'opérateurs quantiques
   • Établir l'hamiltonien moléculaire pour $N_2$.
   • Expliquer le cluster unitaire local de Jastrow (LUCJ) inspiré de la chimie et adapté au matériel [1]
2. Étape 2 : Optimisation 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écution sur le matériel cible
   • Exécuter le circuit optimisé sur une QPU réelle pour générer des échantillons du sous-espace.
4. Étape 4 : Post-traitement des résultats
   • Introduire la boucle de récupération de la configuration autoconsistante [2]
     • Post-traitement de l'ensemble des échantillons de chaînes de bits, en utilisant la connaissance préalable du nombre de particules et de l'occupation orbitale moyenne calculée lors de l'itération la plus récente.
     • Créer de manière probabiliste des lots de sous-échantillons à partir des chaînes de bits récupérées.
     • Projeter et diagonaliser l'hamiltonien moléculaire sur chaque sous-espace échantillonné.
     • Enregistrez l'énergie minimale de l'état fondamental trouvée dans tous les lots et mettez à jour l'occupation orbitale moyenne.

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

• PySCF pour définir la molécule et configurer le hamiltonien.
• ffsim pour construire l'ansatz LUCJ.
• Qiskit pour transposer l'ansatz en vue d'une exécution matérielle.
• Qiskit IBM Runtime pour exécuter le circuit sur une QPU et collecter des échantillons.
• Qiskit addon SQD récupération de la configuration et estimation de l'énergie de l'état fondamental à l'aide de la projection du sous-espace et de la diagonalisation de la matrice.

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{a}^\dagger_{p\sigma}$ / $\hat{a}_{p\sigma}$ sont les opérateurs de création/annihilation fermioniques associés au $p$ -ème élément de base et au spin $\sigma$. $h_{pr}$ et $(pr|qs)$ sont les intégrales électroniques à un et deux corps. En utilisant pySCF,, nous définirons la molécule et calculerons les intégrales à un et deux corps de l'hamiltonien pour l'ensemble de base 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

Dans cette leçon, nous utiliserons la transformation de Jordan-Wigner (JW) pour convertir une fonction d'onde fermionique en une fonction d'onde qubit afin qu'elle puisse être préparée à l'aide d'un circuit quantique. La transformation JW fait correspondre l'espace de Fock des fermions dans M orbitales spatiales à l'espace de Hilbert de 2M qubits, c'est-à-dire qu'une orbitale spatiale est divisée en deux orbitales de spin, l'une associée à un électron de spin haut ( $\alpha$ ) et l'autre à un électron de spin bas ( $\beta$ ). Une orbitale de spin peut être occupée ou inoccupée. En général, lorsque nous parlons du nombre d'orbitales, nous utilisons le nombre d'orbitales spatiales. Le nombre d'orbitales de spin sera doublé. Dans les circuits quantiques, nous représentons chaque orbitale de spin par un qubit. Ainsi, un ensemble de qubits représentera des orbitales de spin-up ou $\alpha$ - et un autre ensemble représentera des orbitales de spin-down ou $\beta$ -. Par exemple, la molécule $N_2$ pour l'ensemble de base 6-31g possède des orbitales spatiales $16$ (c'est-à-dire des orbitales de spin $16$ $\alpha$ + $16$ $\beta$ = $32$ ). Nous aurons donc besoin d'un circuit quantique $32$ -qubit (nous pourrions avoir besoin de qubits ancillaires supplémentaires, comme nous le verrons plus loin). Les qubits sont mesurés dans une base de calcul 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 indiquent l'occupation des électrons dans les orbitales de spin : une adresse $1$ dans une position de bit signifie que l'orbitale de spin correspondante est occupée, tandis qu'une adresse $0$ signifie que l'orbitale de spin est vide. Comme les problèmes de structure électronique préservent les particules, seul un nombre fixe d'orbitales de spin doit être occupé. La molécule $N_2$ possède des électrons $5$ ( $\alpha$ ) et $5$ ( $\beta$ ). Ainsi, toute chaîne de bits représentant les orbitales $\alpha$ et $\beta$ doit comporter cinq $1\text{s}$ pour la molécule $N_2$.

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

Dans cette leçon, nous utiliserons l'anatz de Jastrow (LUCJ) \ [1] pour la préparation d'un état quantique et l'échantillonnage subséquent. Tout d'abord, nous expliquerons les différents éléments constitutifs de l'ansatz UCJ complet et les approximations faites dans la version locale de celui-ci. Ensuite, en utilisant le paquetage ffsim, nous construirons l'ansatz LUCJ et l'optimiserons en utilisant 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)

où $\vert \Phi_{0} \rangle$ est un état de référence, généralement l'état de Hartree-Fock (HF). L'état Hartree-Fock étant défini comme ayant les orbitales les plus basses occupées, la préparation de l'état HF consistera à 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 ascendants et 2 spins descendants peut ressembler à ce qui suit :

Schéma

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)}}$ ) entrecoupée de rotations orbitales ( $e^{K^{(\mu)}}$ et $e^{-K^{(\mu)}}$ ).

Les blocs de rotation orbitale fonctionnent sur une seule espèce de spin ( $\alpha$ (up-spin)/ $\beta$ (down-spin)). Pour chaque espèce d'électron, la rotation orbitale consiste en une couche de portes $R_{z}$ à qubit unique, suivie d'une séquence de portes de rotation de Given à deux qubits (portes $XX + YY$ ).

Les portes à 2 qubits agissent sur des spin-orbitaux adjacents (qubits les plus proches) et peuvent donc être mises en œuvre sur les QPU IBM® sans nécessiter de portes SWAP.

Schéma

Le $e^{iJ^{(\mu)}}$, également connu sous le nom d'opérateur de Coulomb diagonal, se compose de trois blocs. Deux d'entre eux travaillent sur les mêmes secteurs de rotation ( $e^{iJ_{\alpha \alpha}^{(\mu)}}$ et $e^{iJ_{\beta \beta}^{(\mu)}}$ ), et un travaille entre deux secteurs de rotation ( $e^{iJ_{\alpha \beta}^{(\mu)}}$ ).

Tous les blocs de $e^{iJ^{(\mu)}}$ sont constitués de portes numériques $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})$ à qubit unique agissant sur deux qubits distincts.

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

Par exemple, considérons le bloc suivant $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (ou $e^{iJ_{\beta \beta}^{(\mu)}}$ ) pour les orbitales spatiales $N = 4$. Pour une connectivité linéaire des qubits, les trois dernières portes ne sont pas directement réalisables car elles fonctionnent 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 adjacentes (la figure suivante montre un exemple avec $3$ portes SWAP).

Ensuite, le site $J_{\alpha \beta}$ met en œuvre des portes entre les mêmes orbitales indexées 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 une QPU, ces portes nécessiteront également des SWAP.

Il ressort de la discussion ci-dessus que l'ansatz UCJ se heurte à certains obstacles pour l'exécution matérielle, car il nécessite des portes SWAP en raison des interactions entre qubits non adjacents. La variante locale de l'ansatz UCJ, LUCJ, relève ce défi en supprimant certaines $U_{nn}$ de l'opérateur de Coulomb diagonal.

Dans les mêmes blocs d'espèces d'électrons ( $J_{\alpha \alpha}$ et $J_{\beta \beta}$ ), nous ne conservons que les portes $U_{nn}$ compatibles avec la connectivité du plus proche voisin et nous supprimons les portes entre qubits non adjacents dans la version LUCJ. La figure suivante montre le bloc LUCJ après l'élimination des portes non adjacentes.

Ensuite, la version LUCJ du bloc $J_{\alpha \beta}$ qui fonctionne entre différentes espèces d'électrons peut prendre différentes formes en fonction de la topologie du dispositif.

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

• Carré : nous pouvons avoir des portes $U_{nn}$ entre toutes les orbitales $\alpha$ et $\beta$ sans SWAP, et par conséquent, il n'est pas nécessaire de supprimer les portes $U_{nn}$.
• Heavy-hex : Les interactions $\alpha$ - $\beta$ sont maintenues entre chaque $4$ -th indexed (such as the 0th, 4th, and 8th) spin orbitals and are ancilla mediated, that is, we need ancilla qubits between the linear chains representing $\alpha$ and $\beta$ orbitals. Cet arrangement nécessite un nombre limité de SWAP.
• Hexagonale : Toutes les autres orbitales, telles que les orbitales indexées 0e, 2e et 4e, deviennent les plus proches voisines lorsque $\alpha$ et $\beta$ sont disposées en deux chaînes linéaires adjacentes.
• Linéaire : Seules les orbites $\alpha$ et $\beta$ sont connectées, ce qui signifie que le bloc $J_{\alpha \beta}$ n'aura qu'une seule porte.

Schémas de

Si la suppression des portes de l'ansatz UCJ pour construire la version LUCJ le rend plus compatible avec le HW, l'ansatz perd une partie de son 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'approche LUCJ

La 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 grappes couplées simples et doubles (CCSD), où les amplitudes t1 sont le coefficient des opérateurs d'excitation simple et les amplitudes t2 sont celles des opérateurs d'excitation double.

Notez que si l'initialisation de l'ansatz LUCJ avec t1 et t2 amplitudes génère des résultats satisfaisants, les paramètres de l'ansatz peuvent nécessiter une optimisation plus poussée.

# 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 Construction de l'approche LUCJ à l'aide de ffsim

Nous utiliserons le paquetage ffsim pour créer et initialiser l'ansatz avec t1 et t2 amplitudes calculées ci-dessus. Notre molécule ayant un état Hartree-Fock à coquille fermée, nous utiliserons la variante à spin équilibré de l'ansatz UCJ, UCJOpSpinBalanced.

Comme le matériel IBM a une topologie en hexagone lourd, nous adopterons le modèle en zig-zag utilisé dans [1] et expliqué ci-dessus pour les interactions entre qubits. Dans ce schéma, les orbitales (qubits) ayant le même spin sont reliées par une topologie linéaire (cercles rouges et bleus). En raison de la topologie hexagonale lourde, les orbitales de différents spins ont des connexions entre toutes les 4 orbitales, c'est-à-dire la 0e, la 4e, la 8e et ainsi de suite (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. Prenons le cas de n_reps=2. Les deux blocs de rotation orbitale au milieu peuvent être fusionnés en un seul bloc de rotation orbitale. Le paquet ffsim dispose d'un gestionnaire de passage nommé ffsim.qiskit.PRE_INIT pour optimiser le circuit en fusionnant ces blocs adjacents.

2. Optimisation pour le matériel cible

Tout d'abord, nous recherchons un backend de notre choix. Nous optimiserons notre circuit pour le backend, puis nous exécuterons le circuit optimisé sur le même backend afin de générer des échantillons pour le 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 des qubits physiques (initial_layout) dans le matériel cible qui respecte le schéma en zig-zag (deux chaînes linéaires avec un qubit ancilla entre les deux) décrit ci-dessus. La disposition des qubits selon ce schéma permet d'obtenir un circuit efficace compatible avec le matériel et comportant moins de portes.
• Générez un gestionnaire de laissez-passer par étapes en utilisant la fonction generate_preset_pass_manager de Qiskit avec votre choix de backend et initial_layout.
• Réglez l'étape pre_init de votre gestionnaire de passes par étapes sur ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT inclut des passes de transposition Qiskit qui décomposent les portes en rotations orbitales et fusionnent ensuite les rotations orbitales, ce qui permet de réduire le nombre de portes dans le circuit final.
• Lancez le gestionnaire de passes sur votre 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()}")

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

Après avoir optimisé le circuit en vue de son exécution sur le matériel, nous sommes prêts à l'exécuter sur le matériel cible et à collecter des échantillons afin d'estimer l'énergie de l'état fondamental. Comme nous n'avons qu'un seul circuit, nous allons utiliser le et qiskit-ibm-runtimeMode d'exécution des tâches exécuter 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. Résultats post-traitement

La partie post-traitement du flux de travail SQD peut être résumée à l'aide du diagramme suivant.

L'échantillonnage de l'ansatz LUCJ dans la base de calcul génère un ensemble de configurations bruitées $\tilde{\mathcal{\chi}}$, qui sont utilisées dans la routine de post-traitement. Il s'agit d'une méthode appelée (détails discutés plus loin) récupération de configuration pour corriger de manière probabiliste les configurations avec des nombres d'électrons incorrects. Les configurations ne comportant que des numéros d'électrons corrects ( $\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, l'hamiltonien moléculaire, $H$, est projeté sur des sous-espaces :

Chaque hamiltonien projeté $H_{\mathcal{S}^{(k)}}$ est ensuite introduit dans un résolveur d'équations, 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 paquetage qiskit-addon-sqd qui utilise la méthode de Davidson de PySCF pour la diagonalisation.

Nous recueillons ensuite la valeur propre la plus faible (énergie) des lots et calculons également l'occupation orbitale moyenne, $\text{n}$. Les informations relatives à l'occupation moyenne sont utilisées dans l'étape de récupération de la configuration pour corriger de manière probabiliste les configurations bruitées.

Ensuite, nous expliquons en détail la boucle de récupération de la configuration autoconsistante et montrons des exemples de code concrets pour mettre en œuvre les étapes susmentionnées afin d'estimer l'énergie de l'état fondamental de l'hamiltonien $N_2$.

4.1 Récupération de la configuration : aperçu général

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 de spin supérieur et la moitié gauche les orbitales de spin inférieur. 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 exact de particules (l'électron à spin ascendant et l'électron à spin descendant). 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$, nous savons que la chaîne de bits est corrompue par le bruit. La routine de configuration autoconsistante tente de corriger la chaîne de bits en inversant de manière probabiliste les bits $|N_x - N|$ en s'appuyant sur les informations relatives à l'occupation moyenne de l'orbite. L'occupation orbitale moyenne ( $n$ ) indique la probabilité qu'une orbitale soit occupée par un électron. Si $N_x < N$, nous avons moins d'électrons et devons inverser certains $0$ s en $1$ s et vice versa.

La probabilité de retournement peut être $|x[i] - avg\_occupancy[i]|$ pour i-th spin orbital. Dans [2], les auteurs ont utilisé une probabilité pondérée de retournement en utilisant la fonction ReLU modifiée.

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 niveau du coin. Pour $\delta = 0$, $w$ devient la vraie fonction ReLU, et pour $\delta >0$, elle devient ReLU* modifiée*. Dans l'article, les auteurs ont utilisé $\delta = 0.01$ et $h =$ nombre de particules alpha (ou bêta) / nombre d'orbitales de spin alpha (ou bêta) $= 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 par des configurations ne comportant que des nombres corrects de particules dans les deux espèces de spin. Après la première itération, nous disposons d'une estimation de l'état fondamental et, à l'aide de cette estimation, nous pouvons construire la première estimation de $n$. Cette estimation de $n$ est utilisée pour récupérer des configurations, exécuter l'itération suivante de l'estimation de l'état fondamental et affiner de manière autoconsistante l'estimation de $n$. Le processus se répète jusqu'à ce qu'un critère d'arrêt soit rempli.

Prenons l'exemple suivant pour $N = 2$ et $x = |1000\rangle$ ( $N_x = 1$ ). Nous devons faire passer l'un des 0s à 1 pour le corriger en fonction des nombres de particules, et les choix sont 1100, 1010, et 1001. En fonction de la probabilité de renversement, 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 qu'au cours de la première itération, nous exécutions deux lots, et que les états de base estimés à partir de ces lots soient les suivants :

En utilisant les états de base de calcul et leurs amplitudes, nous pouvons calculer la probabilité d'occupation des électrons (en bref, les occupations ) par spin-orbite (qubit) (notez que probabilité = |amplitude| $^2$ ). Nous présentons ci-dessous les occupations par qubit pour chaque chaîne de bits apparaissant dans l'état fondamental estimé et nous calculons l'occupation orbitale totale pour un lot. Notez que, conformément à la convention d'ordonnancement de Qiskit, le bit le plus à droite représente qubit-0 ( Q0 ), et le bit le plus à gauche représente Q3.

Occupation ( Batch0 ) :

Occupation ( Batch1 )

Occupation (moyenne des lots)

En utilisant l'occupation moyenne des orbitales calculée ci-dessus, nous pouvons trouver les probabilités d'inversion pour les différentes orbitales de 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 son p(flip) à $0$. Pour les autres orbitales, qui sont inoccupées, la probabilité de retournement est de $\vert x[i] - \text{n}[i] \vert$ chacune. Avec p(flip), nous calculons également le poids de probabilité associé au flip à l'aide de 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$ )

Enfin, en utilisant les probabilités pondérées ci-dessus, nous pouvons inverser l'une des orbitales inoccupées Q2, Q1 et Q0. Sur la base des valeurs ci-dessus, Q0 sera très probablement inversé, et une configuration récupérée possible peut être $\vert \text{1001} \rangle$.

Le processus complet de récupération de la configuration autoconsistante 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 comprend à la fois des configurations avec un 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 au hasard pour créer des lots $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ de vecteurs pour la projection dans le 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 important, plus la dimension du sous-espace est grande et plus la diagonalisation est exigeante en termes de calcul. D'autre part, un nombre trop faible d'échantillons peut ne pas tenir compte des vecteurs de soutien de l'état fondamental et conduire à une estimation incorrecte.
2. Exécuter le solveur d'états propres (c'est-à-dire la projection sur le sous-espace et la diagonalisation) sur les lots et obtenir des états propres approximatifs. $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. A partir des états propres approximatifs, construire la première estimation de $n$.

Itérations ultérieures :

1. En utilisant $n$, corrigez les configurations dont le nombre de particules est erroné dans $\widetilde{\chi}_{incorrect}$. Supposons que nous les nommions $\widetilde{\chi}_{correct\_new}$. Ensuite, $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ forme le nouvel ensemble de configurations avec les bons numéros 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 fonctionne avec de nouveaux lots et génère de nouvelles estimations des états fondamentaux $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. A partir des états propres approximatifs, construire une estimation affinée pour $n$.
5. Si le critère d'arrêt n'est pas respecté, retourner à l'étape 2.1.

4.2 Estimation de l'état fondamental

Tout d'abord, nous transformerons 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'indices de colonne (N, N/2] (moitié droite), et les orbitales bêta 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)

Quelques options contrôlées par l'utilisateur sont importantes pour cette technique :

• iterations: Nombre d'itérations de récupération de la configuration autoconsistante
• 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 résolveur d'équations

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

Le premier graphique montre qu'après quelques itérations, nous estimons l'énergie de l'état fondamental à ~24 mH (la précision chimique est généralement acceptée comme étant de 1 kcal/mol $\approx$ 1.6 mH ). Le deuxième graphique montre l'occupation moyenne de chaque orbite spatiale après la dernière itération. Nous pouvons constater que les électrons de spin-up et de spin-down occupent les cinq premières orbitales avec une forte probabilité dans nos solutions.

Bien que l'énergie estimée de l'état fondamental soit convenable, elle ne se situe pas dans la limite de la précision chimique ( $\pm \approx 1.6$ mH ). Cet écart peut être attribué à la petite dimension du sous-espace que nous avons utilisée ci-dessus pour la projection et la diagonalisation. Comme nous avons utilisé samples_per_batch=500, le sous-espace est couvert par le maximum de vecteurs $500$, c'est-à-dire les vecteurs manquants du support de l'état fondamental. L'augmentation du paramètre samples_per_batch devrait améliorer la précision au détriment des ressources informatiques classiques et du 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()

Output:

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

Exercice pour le lecteur

Augmentez progressivement le paramètre samples_per_batch (par exemple, de $1000$ à $10000$ avec un pas de $1000$; autorisé par la mémoire de votre ordinateur) et comparez les énergies estimées de l'état fondamental.

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

[2] J. Robledo-Moreno et al, "La chimie au-delà des solutions exactes sur un superordinateur centré sur les quanta" (2024). arXiv:quant-ph/2405.05068.