Algorithme SqDRIFT pour l'estimation de l'état fondamental

Estimation d'utilisation : 180 secondes sur un processeur Heron r3 (REMARQUE : il s'agit uniquement d'une estimation. Ton temps d'exécution peut varier.)

Objectifs d'apprentissage

• Apprendre à créer des circuits de plus faible profondeur que la trotterisation

• Parcourir un workflow de bout en bout pour l'estimation de l'état fondamental avec qDRIFT et SQD

• Apprendre à utiliser qiskit-fermions conjointement avec d'autres addons Qiskit pour mettre en œuvre un tel workflow

Ce tutoriel est présenté sous forme de notebook Python à des fins pédagogiques.

Prérequis

• Lis la présentation de la diagonalisation quantique par échantillonnage (SQD)

• Lis la leçon sur la diagonalisation quantique de Krylov par échantillonnage (SKQD)

Contexte

SqDRIFT est une variante de SKQD qui remplace la nécessité de choisir un ansatz à partir duquel échantillonner des chaînes de bits par un ensemble de circuits d'évolution temporelle construits directement à partir du hamiltonien cible. Pour cela, on sous-échantillonne des opérateurs d'évolution temporelle plus petits à partir du hamiltonien selon ses coefficients, ce qui est connu sous le nom de méthode de trotterisation qDRIFT.

Ce tutoriel utilise Qiskit Fermions pour créer les circuits fermioniques plus naturels pour l'algorithme qDRIFT, puis des passes de placement et de synthèse fermioniques avant d'intégrer les circuits dans le pipeline Qiskit traditionnel pour l'exécution sur le matériel.

Soit un hamiltonien de la forme :

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

où, sans perte de généralité, nous exigeons $c_i > 0$ et que la plus grande valeur propre de $h_i$ soit égale, en valeur absolue, à $1$. Tout préfacteur signé ou complexe est absorbé dans $h_i$, de sorte que les coefficients $c_i$ sont des poids strictement positifs tandis que les $h_i$ portent la direction de chaque terme. Ici $N$ est le nombre de termes (ou, après regroupement, le nombre de groupes) du hamiltonien ; c'est une propriété du hamiltonien, distincte du nombre d'opérateurs échantillonnés dans un seul circuit, noté $n$ ci-dessous.

L'algorithme qDRIFT réalise alors, pour le temps cible $t$, un opérateur $V_k$, où $k$ va de $1 \cdots K$ et désigne le $k_{th}$ circuit SqDRIFT, défini par :

$$V_k = \prod_{j=1}^{n} e^{-i h_{k_j} \lambda t / n }$$

Ici $n$ est le nombre d'opérateurs échantillonnés par circuit et $K$ est le nombre de circuits de l'ensemble. Le produit porte sur les $n$ tirages, et non sur les $N$ termes du hamiltonien, et comme les termes sont tirés avec remise, le même $h_i$ peut apparaître plusieurs fois dans un même $V_k$.

La quantité :

$$\lambda = \sum_{i=1}^{N} c_i$$

est la norme $L_1$ des coefficients, de sorte que chacune des $n$ étapes évolue pendant la même durée $\lambda t / n$ quel que soit le terme tiré. L'uniformité de l'angle de pas est la caractéristique de qDRIFT : un coefficient influence le résultat par la fréquence à laquelle son terme est tiré, et non par l'ampleur de la rotation de ce terme. Les indices sont échantillonnés selon la distribution :

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

de sorte que la série $(k_1, \ldots, k_n)$ est une séquence aléatoire d'indices de termes tirés selon cette distribution. Comme les $c_i$ sont positifs et de somme $\lambda$, il s'agit d'une distribution de probabilité normalisée, et l'espérance du canal obtenu sur les tirages aléatoires approche l'évolution sous $H$, avec une erreur qui diminue quand $n$ augmente. Note que l'erreur d'approximation dépend de $\lambda$ et non du nombre de termes $N$.

(L'article SqDRIFT note le nombre de termes $\mathcal{N}$ et la longueur de la séquence $N$ ; nous utilisons ici $N$ et $n$ pour bien distinguer les deux.)

Ce tutoriel montre comment générer un ensemble de tels circuits aléatoires. Une fois ces circuits créés, de manière analogue à la création d'un sous-espace de Krylov pour différents opérateurs, nous échantillonnons des chaînes de bits à partir de plusieurs de ces opérateurs avec différents paramètres de temps. Cela garantit un meilleur recouvrement entre les vecteurs de l'état fondamental et les chaînes de bits échantillonnées.

Exigences

Avant de commencer ce tutoriel, assure-toi d'avoir installé

• Un environnement virtuel Python (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (note que le nom est au pluriel)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Tu peux installer tous les paquets requis avec :

pip install "qiskit~=2.5" "qiskit-fermions==0.1.0" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy

Configuration

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "qiskit_aer": "qiskit-aer", "qiskit_fermions": "qiskit-fermions", "qiskit_ibm_runtime": "qiskit-ibm-runtime"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

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

# Qiskit Aer
from qiskit_aer import AerSimulator

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

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

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

Exemple de simulateur

Étape 1 : Convertir les entrées classiques en un problème quantique

Lecture et préparation du FCIDump

Pour ce tutoriel, nous allons charger le hamiltonien de structure électronique de l'azote (N2). Il existe d'autres moyens de créer des opérateurs fermioniques. Consulte la documentation de qiskit_fermions.operators.library.

À propos de ce FCIDump. Le fichier N2_sto_3g décrit une molécule d'azote ($N_2$) dans la base minimale STO-3G à une séparation interatomique de 1,09 $\AA$, la longueur de liaison d'équilibre expérimentale. Son en-tête déclare NORB=10, NELEC=14 et MS2=0 : 10 orbitales spatiales (donc 20 spin-orbitales, et 20 qubits sous Jordan-Wigner), 14 électrons dans un singulet de spin, soit sept électrons $\alpha$ et sept électrons $\beta$. Toutes les orbitales portent l'étiquette de symétrie 1, c'est-à-dire qu'aucune symétrie de groupe ponctuel n'est exploitée. S'agissant d'un dump STO-3G de l'espace complet, aucune orbitale n'est gelée et l'espace de corrélation est assez petit pour qu'une énergie de référence FCI exacte puisse être calculée classiquement à titre de comparaison, comme indiqué dans la cellule suivante.

Un fichier équivalent peut être régénéré avec PySCF :

from pyscf import gto, scf, tools

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

Comme les intégrales dépendent des orbitales SCF convergées, un fichier régénéré peut différer de celui fourni par la phase ou l'ordre des orbitales ; les énergies totales ne sont pas affectées.

Obtenir le fichier. Trouve le FCIDump dans ce dépôt GitHub. Tu peux exécuter la cellule ci-dessous pour le récupérer à l'emplacement attendu par le reste du tutoriel.

Nous utilisons d'abord le cisolver fourni par pyscf pour obtenir l'énergie de référence. C'est la véritable énergie de l'état fondamental de la molécule avec laquelle nous travaillons. Pour cela, nous déclarons d'abord norb et nelec, qui sont respectivement le nombre d'orbitales et le nombre d'électrons. Nous déclarons ensuite h1e et h2e, qui sont respectivement les intégrales à un et à deux électrons. Tout cela servira aussi plus tard pour SQD.

import os
from urllib.request import urlopen

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

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

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

Chargement du hamiltonien

Les données nécessaires étant prêtes, nous lisons le hamiltonien à partir du fichier FCI dans un format compatible avec qiskit-fermions

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

Workflows fermioniques avec qiskit-fermions

Nous allons d'abord convertir le hamiltonien en un modèle de circuit fermionique avec qiskit-fermions, qui fournit des passes de transpilation et des portes spécifiques aux circuits fermioniques. Elles seront utilisées plus tard avant les passes de transpilation traditionnelles de Qiskit pour ce workflow.

Regroupement des termes

Pour garantir la reproductibilité des résultats, nous utilisons d'abord canonical_order pour trier les termes uniquement selon leur structure. L'ordre des opérateurs dans la liste canon est donc fixe. Cela garantit la reproductibilité des opérateurs créés, car la passe QDriftTrotterization que nous utiliserons plus tard échantillonne des indices aléatoires pour créer les opérateurs qDRIFT.

Dans cette étape, nous exploitons les nombreuses symétries présentes dans le hamiltonien de structure électronique en regroupant les termes apparentés ayant des coefficients identiques. Bien que cela modifie la distribution des coefficients des opérateurs dont le protocole qDRIFT tire ses échantillons, cela n'affecte pas ses garanties de convergence. Surtout, regrouper les termes apparentés par symétrie entraîne une annulation favorable de termes de Pauli et une profondeur de circuit globalement plus faible lors de l'évolution temporelle d'un état sous leur action.

qiskit-fermions fournit la fonction group_terms_by_electronic_structure qui effectue ce regroupement pour nous.

Note que group_terms_by_electronic_structure suppose des termes en ordre normal.

Filtrage des termes diagonaux

Nous retirons les termes diagonaux du hamiltonien utilisé pour générer les circuits, afin que les $n$ emplacements d'échantillonnage qDRIFT soient consacrés aux termes qui déplacent la population entre configurations. Il vaut mieux filtrer ces termes du hamiltonien à ce stade, avant la construction de la porte Evolution à l'étape suivante.

Les termes en question sont ceux qui sont diagonaux dans la base des nombres d'occupation, c'est-à-dire les produits d'opérateurs nombre $a^\dagger_i a_i$. Trois types de termes relèvent de cette description :

• le décalage d'énergie constant, un produit de zéro opérateur nombre, dont l'évolution temporelle ne contribue qu'à une phase globale ;

• les opérateurs nombre individuels $n_i$, dont l'évolution temporelle se réduit à des rotations $Z$ à un qubit ;

• les produits d'ordre supérieur tels que $n_i n_j$.

Pris isolément, aucun de ces termes ne déplace la population entre configurations de nombres d'occupation ; ils n'agissent que sur les phases des configurations déjà présentes. Ils ne sont cependant pas inertes : ces phases relatives alimentent l'interférence générée par les termes d'excitation plus loin dans le circuit, donc les filtrer modifie l'évolution effectivement générée et peut modifier la distribution d'échantillonnage. C'est une approximation délibérée dans l'étape de génération des circuits, faite pour concentrer l'échantillonnage sur les termes d'excitation, et non une étape qui laisse la distribution échantillonnée intacte. Contrairement au regroupement par symétrie ci-dessus, qui laisse intactes les garanties de convergence de qDRIFT, ce filtre modifie l'opérateur évolué. Les circuits n'approchent donc plus l'évolution sous le hamiltonien complet, et les bornes d'erreur de qDRIFT s'appliquent à l'opérateur filtré et non à l'original. C'est acceptable ici car les circuits ne sont qu'une heuristique d'échantillonnage servant à proposer des configurations : aucun terme n'est perdu dans l'estimation de l'énergie elle-même, puisque le filtre ne s'applique qu'au hamiltonien utilisé pour construire les circuits, tandis que la diagonalisation classique ultérieure utilise le hamiltonien complet, termes diagonaux inclus. La précision de SQD dépend de cette étape classique, qui reste variationnelle dans le sous-espace échantillonné quelle que soit la manière dont les configurations ont été proposées.

La fonction filter_diagonal_terms() retire ces termes d'un opérateur sur place. Elle les identifie à partir de leur structure en ordre normal — le multiensemble des modes de création correspondant au multiensemble des modes d'annihilation — elle n'est donc valide que sur un opérateur déjà en ordre normal. Cette hypothèse n'est pas vérifiée à l'exécution.

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

print(len(canon.groups))

5060

Maintenant que nous avons regroupé les termes du hamiltonien, nous allons décider des paramètres suivants pour générer l'ensemble de circuits :

• Le nombre de circuits à générer : num_circuits
• La longueur de chaque circuit en groupes d'excitation : num_exc
• Le facteur pour les différents temps d'évolution : times

Création des circuits fermioniques

Nous allons maintenant créer des circuits fermioniques pour chacun des pas de temps. Chaque circuit sera constitué d'une seule porte d'évolution, avec le temps d'évolution que nous avons déclaré plus tôt. L'opérateur d'évolution est le hamiltonien. Plus tard, nous exécuterons des passes de transpilation sur ces circuits pour créer des circuits qDRIFT.

Préparation de l'ansatz

Nous préparons l'état de Hartree-Fock avec la classe InitializeModes. Pour l'azote, le processus consiste simplement à appliquer des portes X aux num_elec_a premiers qubits, puis aux num_elec_b qubits, les deux valant sept pour l'azote. Cet état représente les sept électrons $\alpha$ et les sept électrons $\beta$ de l'azote.

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

init_circuits = []

hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

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

Étape 2 : Optimiser le problème pour l'exécution sur matériel quantique

Maintenant que nous avons nos circuits, nous allons d'abord utiliser les passes disponibles dans qiskit-fermions pour effectuer des optimisations au niveau fermionique, puis transpiler notre circuit pour le backend choisi. Comme il s'agit d'une expérience sur simulateur, nous le ferons d'abord pour l'AerSimulator. Calcul du poids de chaque groupe

Dans cette étape, nous effectuons l'échantillonnage qDRIFT des termes de façon stochastique, avec des probabilités proportionnelles à leurs coefficients dans l'hamiltonien. La passe de transpilation qDRIFT s'en charge pour nous. Nous pouvons maintenant créer des circuits moins profonds, qui s'exécutent plus efficacement sur le matériel malgré une connectivité limitée des qubits, même lorsque l'hamiltonien contient des couplages à longue portée et des termes d'ordre supérieur à deux. Après le regroupement des termes, elle échantillonne les opérateurs selon leurs poids. Pour chaque opérateur $h_i$, le poids $W_{h_i}$ est défini comme suit :

$$W_{h_i} = |c_i| / \lambda$$

Comme les termes ont été regroupés à l'étape 1, chaque $h_i$ désigne ici un groupe entier : $c_i$ est le coefficient absolu moyen des termes du groupe $i$, et chaque terme du groupe est évolué avec son coefficient réduit à son signe.

Optimisations fermioniques et natives du matériel

La fonction generate_preset_jw_pass_manager() renvoie un MultiStagePassManager qui prend un FermionicCircuit et produit un circuit final optimisé que nous pouvons transpiler pour l'exécuter sur notre matériel. Nous remplaçons son étape d'optimisation par défaut par un FermionicPassManager contenant notre passe QDriftTrotterization :

• La passe QDriftTrotterization utilise en interne le calcul des poids et l'échantillonnage pour générer les circuits que nous utiliserons pour l'échantillonnage

• La passe RelabelModes est une autre passe d'optimisation qui permet de permuter les modes fermioniques afin d'optimiser la connectivité entre les qubits et de réduire la profondeur en portes ; pour en savoir plus, consulte la référence de l'API

Les étapes restantes du MultiStagePassManager s'exécutent automatiquement et gèrent l'ensemble de la correspondance fermion-qubit :

• F2QLayout : le pass manager prédéfini applique la passe TrivialF2QLayout, qui associe trivialement $n$ bits fermioniques à $n$ qubits.

• F2QSynth : une passe de transpilation qui convertit les instructions de circuit fermioniques en instructions à base de qubits.

qdrift = QDriftTrotterization(num_exc, rng=19)

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

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

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

400

Maintenant que les optimisations au niveau fermionique sont terminées, nous pouvons transpiler les circuits pour les exécuter sur le simulateur.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Étape 3 : exécuter avec les primitives Qiskit

Maintenant que nous avons nos circuits, nous pouvons les exécuter avec les primitives Qiskit sur l'AerSimulator. Nous allons combiner tous les comptages issus des différents circuits. Nous les convertissons en vecteurs booléens avant le post-traitement final avec SQD.

print(
    f"Executing {len(transpiled_circuits)} circuits with {shots} shots each..."
)

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

all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]

print(len(all_counts), "length before post processing")

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

Étape 4 : post-traiter et renvoyer le résultat dans le format classique souhaité

Utilisation des chaînes de bits pour SQD

Nous pouvons maintenant exécuter le schéma de diagonalisation sur les chaînes de bits sélectionnées pour trouver la plus petite valeur propre, qui correspond à l'énergie de l'état fondamental de la molécule. Nous créons une fonction de rappel, déclarons les occupations initiales et fixons les paramètres avant d'exécuter finalement le schéma de diagonalisation. La fonction de rappel sert à afficher l'itération courante et l'estimation courante de la valeur propre à chaque itération.

Enfin, pour obtenir l'estimation de l'état fondamental, nous ajoutons nuclear_repulsion_energy à l'énergie obtenue.

Remarque : la dimension du sous-espace n'est pas fixe d'une itération à l'autre, même sur le simulateur sans bruit : chaque sous-échantillon tire un ensemble différent de configurations, et l'étape de récupération remodèle le pool entre les itérations, de sorte que la dimension rapportée varie d'un sous-échantillon à l'autre. Un échantillonnage sans bruit ne fixe donc pas, à lui seul, la dimension du sous-espace sélectionné. L'exécution sur le matériel tend en revanche à donner des sous-espaces systématiquement plus grands, car les shots bruités brisent la symétrie du nombre de particules et la récupération de configurations les transforme en vecteurs de base supplémentaires. C'est pourquoi nous introduirons aussi une étape supplémentaire d'élagage des chaînes de bits dans la section consacrée au matériel.

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

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

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

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

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

result_history = []

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

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

computed_energy = result.energy + nuclear_repulsion_energy

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

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

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

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

Exemple sur le matériel

Cet exemple utilise 20 qubits (10 orbitales spatiales). Ce choix est une commodité pour un tutoriel qui doit s'exécuter rapidement, et non une limite stricte de la méthode.

Le coût de l'étape classique n'est pas fixé directement par le nombre de qubits. SQD diagonalise l'hamiltonien projeté sur le sous-espace engendré par les configurations échantillonnées ; ce qui détermine le coût classique est donc la dimension de ce sous-espace sélectionné — gouvernée ici par samples_per_batch, num_batches et le nombre de configurations distinctes que les circuits produisent réellement — ainsi que l'algèbre linéaire creuse nécessaire pour appliquer l'hamiltonien projeté. L'espace CI complet croît de façon combinatoire avec les orbitales et les électrons, mais le sous-espace sélectionné en est une petite tranche ajustable, dont nous contrôlons directement la taille. Par conséquent, le nombre de qubits et la difficulté classique peuvent varier de façon assez indépendante : un espace orbital plus large échantillonné dans un sous-espace modeste peut coûter moins cher qu'un système plus petit diagonalisé sur un sous-espace très grand.

En pratique, la taille de système réalisable dépend donc de la dimension de sous-espace dont tu as besoin pour la précision souhaitée, ainsi que de la mémoire et des cœurs disponibles pour le solveur de valeurs propres. Les espaces orbitaux plus grands exigent généralement un sous-espace plus grand pour atteindre la précision chimique, et c'est ce qui motive finalement le recours à des ressources distribuées — consulte qiskit-addon-sqd-hpc pour déployer cette étape à plus grande échelle. Plutôt que de supposer un seuil fixe, l'approche pratique consiste à surveiller la dimension de sous-espace rapportée et la convergence de l'énergie au fil des itérations, et à augmenter la taille du sous-espace jusqu'à ce que l'énergie cesse de s'améliorer ou que tu épuises la mémoire disponible.

Remarque : en raison de l'erreur d'échantillonnage due au bruit du matériel, le sous-espace créé pour la diagonalisation lors de l'exécution sur le matériel sera plus grand que celui obtenu avec le simulateur. Même si cela augmente la dimension du sous-espace que nous voulons diagonaliser, le flux de travail nous donne quand même une réponse précise grâce à la robustesse de SQD face au bruit.

Élagage des chaînes parasites

Ici, nous pouvons choisir d'effectuer une étape supplémentaire. Lorsque nous disposons de toutes les chaînes de bits issues des exécutions des circuits, nous pouvons soit filtrer les chaînes de bits invalides avant d'exécuter SQD, soit continuer sans élagage. Ne pas élaguer est généralement préférable pour les exécutions sur le matériel, car cela laisse les shots à symétrie brisée disponibles pour la récupération de configurations, qui peut les réparer en configurations valides et ainsi élargir le sous-espace au lieu de rejeter purement et simplement ces shots.

Comme l'azote ne peut avoir que sept électrons $\alpha$ et sept électrons $\beta$, toute chaîne de bits qui compte plus ou moins de sept 1 dans la première et la seconde moitié de la sortie peut être écartée. Nous définissons une fonction qui vérifie si les chaînes de bits sont valides et, sinon, les écarte. Une fois les chaînes de bits parasites filtrées, les autres sont envoyées au schéma de diagonalisation. Utilise l'indicateur PRUNE ci-dessous pour basculer entre les deux comportements.

Garde à l'esprit que l'élagage n'est qu'un des nombreux choix qui façonnent le sous-espace final, avec le nombre de circuits, l'ensemble des temps d'évolution et le filtrage des termes diagonaux. Comparer une exécution avec élagage à une exécution sans élagage n'est instructif que si tout le reste est maintenu fixe ; la version C++ de ce tutoriel en discute plus en détail, car elle post-sélectionne au lieu de récupérer et diffère aussi sur ces autres paramètres.

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

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

print(len(canon.groups))

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

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

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

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

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

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

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

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

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

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

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

shots = 100

sampler = Sampler(mode=backend)

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

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

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

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

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

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

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

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

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

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

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

result_history = []

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

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

computed_energy = result.energy + nuclear_repulsion_energy

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

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

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

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

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

Prochaines étapes

Recommandations

Si ce travail t'a intéressé, les ressources suivantes pourraient te plaire :

• Diagonalisation quantique de Krylov basée sur l'échantillonnage d'un modèle de réseau fermionique - un tutoriel connexe qui utilise des circuits d'évolution temporelle au lieu d'un ansatz variationnel.
• Diagonalisation quantique basée sur l'échantillonnage d'un hamiltonien de chimie - un tutoriel sur la construction d'un circuit LUCJ (local unitary cluster Jastrow) pour la simulation de chimie quantique.
• L'article SqDRIFT - la publication sur laquelle repose ce tutoriel. (Note que certaines des optimisations décrites dans cet article sont actuellement en cours de développement, et que ce tutoriel est susceptible d'évoluer en fonction de l'évolution des bibliothèques utilisées.)