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

Python

Télécharger le notebook Python

Télécharge le notebook Python pour exécuter toi-même ce tutoriel. Le notebook télécharge l'entrée FCIDump si elle n'est pas déjà présente.

Estimation d'utilisation : 180 secondes sur un processeur Heron r3 (REMARQUE : ceci n'est qu'une estimation. Ton temps d'exécution peut varier.)

Résultats d'apprentissage

• Apprends à créer des circuits de plus faible profondeur par rapport à la trotterisation

• Parcours un flux de travail de bout en bout pour l'estimation de l'état fondamental à l'aide de qDRIFT et SQD

• Apprends à utiliser qiskit-fermions en tandem avec d'autres addons Qiskit pour implémenter un tel flux de travail

Prérequis

• Lis l'aperçu de Sample-based quantum diagonalization (SQD)

• Lis la leçon Sample-based Krylov Quantum Diagonalization (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. Ceci est réalisé en sous-échantillonnant des opérateurs d'évolution temporelle plus petits à partir du hamiltonien en fonction de ses coefficients, ce qui est connu sous le nom de méthode de trotterisation qDRIFT.

Les deux implémentations utilisent Qiskit Fermions pour représenter l'hamiltonien moléculaire et construire les circuits de l'algorithme qDRIFT. L'implémentation Python utilise des passes de transpilation fermioniques, tandis que l'implémentation C++ échantillonne explicitement des groupes et les associe à des opérateurs de qubits.

Soit le hamiltonien de la forme :

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

où, sans perte de généralité, nous exigeons que $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) dans le 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 certain opérateur $V_k$, où $k$ va de $1 \cdots K$ et désigne le $k_{th}$ circuit SqDRIFT, défini comme :

$$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 dans l'ensemble. Le produit porte sur les $n$ tirages, et non sur les $N$ termes du hamiltonien au total, et comme les termes sont tirés avec remise, le même $h_i$ peut apparaître plus d'une fois dans un seul $V_k$.

La quantité :

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

est la norme $L_1$ des coefficients, donc 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 des étapes est la caractéristique déterminante 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 à partir de la distribution :

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

donc la série $(k_1, \ldots, k_n)$ est une séquence aléatoire d'indices de termes tirés de cette distribution. Comme les $c_i$ sont positifs et somment à $\lambda$, il s'agit d'une distribution de probabilité normalisée, et l'espérance du canal résultant sur les tirages aléatoires approxime l'évolution sous $H$, avec une erreur qui diminue à mesure que $n$ augmente. Remarque que l'erreur d'approximation dépend de $\lambda$ plutôt que 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 garder les deux clairement distincts.)

Ce tutoriel montre comment générer un ensemble de tels circuits randomisés. Une fois ces circuits créés, de façon similaire à 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 temporels. Cela garantit un meilleur recouvrement entre les vecteurs d'é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 (remarque 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

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

Implémentation

Étape 1 : associer les entrées classiques à un problème quantique

Les étapes Python ci-dessous utilisent le simulateur Aer. L'exemple sur matériel suit ensuite.

Lecture et préparation du FCIDump

Pour ce tutoriel, nous allons charger le hamiltonien de structure électronique pour l'azote (N2). Il existe aussi d'autres façons de créer des opérateurs fermioniques. Consulte la documentation à 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 orbitales de spin, et 20 qubits sous Jordan-Wigner), 14 électrons dans un singulet de spin, donc sept électrons $\alpha$ et sept électrons $\beta$. Toutes les orbitales reçoivent 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 en 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 à des fins de comparaison, comme montré dans le bloc de code suivant.

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 du fichier fourni en phase ou en 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 le code 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. Il s'agit de 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. Ensuite, nous déclarons h1e et h2e, qui sont respectivement les intégrales à un et à deux électrons. Toutes ces valeurs seront également utilisées 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}")

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

Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

Chargement du hamiltonien

Avec les données nécessaires prêtes, nous lisons l'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

Flux de travail fermioniques avec qiskit-fermions

Nous allons d'abord associer l'hamiltonien à un modèle de circuit fermionique à l'aide de qiskit-fermions, qui fournit des passes de transpileur et des portes spécifiques aux circuits fermioniques. Celles-ci seront utilisées plus tard avant les passes traditionnelles du transpileur de Qiskit pour ce flux de travail.

Regroupement de termes

Pour garantir la reproductibilité des résultats, nous utilisons d'abord canonical_order pour trier les termes en fonction uniquement de 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 de l'opérateur dans laquelle le protocole qDRIFT effectue ses tirages, cela n'affecte pas ses garanties de convergence. Fait crucial, le regroupement de termes liés par symétrie entraîne une annulation favorable des termes de Pauli et une profondeur de circuit globalement plus courte 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 de l'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 les configurations. De tels termes sont mieux filtrés hors de l'hamiltonien à ce stade, avant que la porte Evolution ne soit construite à l'étape suivante.

Les termes en question sont ceux qui sont diagonaux dans la base du nombre d'occupation, c'est-à-dire les produits d'opérateurs de 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 de nombre, dont l'évolution temporelle ne contribue qu'à une phase globale ;

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

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

À elles seules, aucune de ces contributions ne déplace la population entre les configurations du nombre d'occupation ; elles agissent uniquement sur les phases des configurations déjà présentes. Elles 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, si bien que les filtrer modifie l'évolution effectivement générée et peut changer la distribution d'échantillonnage. Il s'agit d'une approximation délibérée dans l'étape de génération de circuit, faite pour concentrer l'échantillonnage sur les termes d'excitation, plutôt que d'une étape qui laisse intacte la distribution échantillonnée. Contrairement au regroupement par symétrie ci-dessus, qui laisse intactes les garanties de convergence de qDRIFT, ce filtre modifie l'opérateur en cours d'évolution. Les circuits n'approximent donc plus l'évolution sous l'hamiltonien complet, et les bornes d'erreur de qDRIFT s'appliquent à l'opérateur filtré plutôt qu'à l'original. Ceci est acceptable ici car les circuits ne constituent qu'une heuristique d'échantillonnage utilisée pour proposer des configurations : aucun terme n'est perdu pour l'estimation d'énergie elle-même, puisque le filtre s'applique uniquement à l'hamiltonien utilisé pour construire les circuits, tandis que la diagonalisation classique ultérieure utilise l'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() supprime de tels 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 dans l'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 termes de groupes d'excitation : num_exc
• Le facteur pour les différents temps d'évolution : times

Création de circuits fermioniques

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

Préparation de l'ansatz

Nous préparons l'état de Hartree-Fock à l'aide de la classe InitializeModes. Pour l'azote, le processus consiste simplement à appliquer des portes X aux premiers num_elec_a qubits, puis aux num_elec_b qubits, tous deux égaux à 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 de simulation, nous ferons d'abord cela pour l'AerSimulator.

Calcul du poids pour chaque groupe

Dans cette étape, nous effectuons l'échantillonnage qDRIFT des termes de manière stochastique avec des probabilités proportionnelles à leurs coefficients dans l'hamiltonien. La passe de transpileur qDRIFT fait cela pour nous. Nous pouvons maintenant créer des circuits moins profonds qui peuvent être exécutés sur le matériel plus efficacement malgré une connectivité de qubits limitée, même lorsque l'hamiltonien contient des couplages à longue portée et des termes d'ordre supérieur au quadratique. Après le regroupement des termes, il échantillonne les opérateurs en fonction de 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$ ici est 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 au matériel

La fonction generate_preset_jw_pass_manager() retourne 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 peut être utilisée pour permuter les modes fermioniques afin d'optimiser la connectivité entre les qubits et réduire la profondeur des 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 le mappage complet de fermion à qubit :

• F2QLayout : Le gestionnaire de passes préconfiguré applique la passe TrivialF2QLayout, qui associe trivialement $n$ bits fermioniques à $n$ qubits.

• F2QSynth : Une passe de transpilation pour associer des instructions de circuit basées sur les fermions à des instructions basées sur les 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 nous avons terminé les optimisations au niveau fermionique, nous pouvons transpiler les circuits pour l'exécution sur le simulateur.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Étape 3 : Exécuter à l'aide des primitives Qiskit

Maintenant que nous avons nos circuits, nous pouvons les exécuter à l'aide des primitives Qiskit sur l'AerSimulator. Nous allons combiner tous les comptes 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 retourner 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 valeur propre la plus basse, qui correspondra à l'énergie de l'état fondamental de la molécule. Nous créons une fonction de rappel, déclarons les occupations initiales et définissons les paramètres avant d'exécuter enfin le schéma de diagonalisation. La fonction de rappel est utilisée pour afficher l'itération en cours et l'estimation actuelle de la valeur propre à chaque itération.

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

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 réservoir entre les itérations, si bien que la dimension rapportée varie d'un sous-échantillon à l'autre. L'échantillonnage sans bruit ne fixe pas en soi la dimension du sous-espace sélectionné. L'exécution sur matériel, en revanche, tend à donner des sous-espaces systématiquement plus grands, car les tirages bruités brisent la symétrie du nombre de particules et la récupération de configuration les transforme en vecteurs de base supplémentaires. Pour cette raison, nous introduirons également une autre étape d'élagage des chaînes de bits dans la section matérielle.

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 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 un plafond strict pour la méthode.

Le coût de l'étape classique n'est pas déterminé directement par le nombre de qubits. SQD diagonalise l'hamiltonien projeté sur le sous-espace couvert par les configurations échantillonnées, si bien que ce qui détermine le coût classique est la dimension de ce sous-espace sélectionné — régie 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 manière combinatoire avec les orbitales et les électrons, mais le sous-espace sélectionné en est une tranche petite et ajustable, dont nous contrôlons directement la taille. Par conséquent, le nombre de qubits et la difficulté classique peuvent varier de manière quelque peu indépendante : un espace orbital plus large échantillonné dans un sous-espace modeste peut être moins coûteux 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 nécessitent 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 — voir qiskit-addon-sqd-hpc pour faire évoluer 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 matériel sera plus grand que celui obtenu avec le simulateur. Bien que cela augmente la dimension du sous-espace que nous voulons diagonaliser, le flux de travail nous donne toujours 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 de circuits, nous pouvons soit filtrer les chaînes de bits invalides avant d'exécuter SQD, soit poursuivre sans élagage. Il est généralement préférable de sauter l'élagage pour les exécutions sur matériel, car cela laisse les tirages ayant brisé la symétrie disponibles pour la récupération de configuration, qui peut les réparer en configurations valides et ainsi élargir le sous-espace au lieu de rejeter purement et simplement ces tirages.

Comme l'azote ne peut avoir que sept électrons $\alpha$ et sept électrons $\beta$, toute chaîne de bits dont la première ou la seconde moitié de la sortie ne comporte pas exactement sept 1 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, le reste est envoyé dans le 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, aux côtés du nombre de circuits, de l'ensemble des temps d'évolution et du filtrage des termes diagonaux. Comparer une exécution élaguée à une exécution non élaguée n'est informatif que si tout le reste est maintenu fixe ; l'onglet C++ traite cela plus en détail, puisqu'il effectue une post-sélection plutôt qu'une récupération et diffère également 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")

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

C++

Obtenir le code C++

Le code source C++, le projet CMake et les instructions de compilation se trouvent dans le répertoire docs/tutorials/assets/sqdrift/cpp du dépôt de documentation de Qiskit. L'entrée FCIDump se trouve dans le répertoire voisin fcidump_files. Clone le dépôt, ou télécharge les deux répertoires, pour compiler et exécuter toi-même ce tutoriel.

Estimation d'utilisation : 180 secondes sur un processeur Heron r3 (REMARQUE : ceci n'est qu'une estimation. Ton temps d'exécution peut varier.)

Résultats d'apprentissage (C++)

• Apprends à créer des circuits de plus faible profondeur par rapport à la trotterisation

• Parcours un flux de travail de bout en bout pour l'estimation de l'état fondamental à l'aide de qDRIFT et SQD

• Apprends à utiliser qiskit-fermions en tandem avec d'autres addons Qiskit pour implémenter un tel flux de travail

Prérequis (C++)

• Lis l'aperçu de Sample-based quantum diagonalization (SQD)

• Lis la leçon Sample-based Krylov Quantum Diagonalization (SKQD)

Contexte (C++)

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. Ceci est réalisé en sous-échantillonnant des opérateurs d'évolution temporelle plus petits à partir du hamiltonien en fonction de ses coefficients, ce qui est connu sous le nom de méthode de trotterisation qDRIFT.

Les deux implémentations utilisent Qiskit Fermions pour représenter l'hamiltonien moléculaire et construire les circuits de l'algorithme qDRIFT. L'implémentation Python utilise des passes de transpilation fermioniques, tandis que l'implémentation C++ échantillonne explicitement des groupes et les associe à des opérateurs de qubits.

Soit le hamiltonien de la forme :

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

où, sans perte de généralité, nous exigeons que $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) dans le 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 certain opérateur $V_k$, où $k$ va de $1 \cdots K$ et désigne le $k_{th}$ circuit SqDRIFT, défini comme :

$$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 dans l'ensemble. Le produit porte sur les $n$ tirages, et non sur les $N$ termes du hamiltonien au total, et comme les termes sont tirés avec remise, le même $h_i$ peut apparaître plus d'une fois dans un seul $V_k$.

La quantité :

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

est la norme $L_1$ des coefficients, donc 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 des étapes est la caractéristique déterminante 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 à partir de la distribution :

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

donc la série $(k_1, \ldots, k_n)$ est une séquence aléatoire d'indices de termes tirés de cette distribution. Comme les $c_i$ sont positifs et somment à $\lambda$, il s'agit d'une distribution de probabilité normalisée, et l'espérance du canal résultant sur les tirages aléatoires approxime l'évolution sous $H$, avec une erreur qui diminue à mesure que $n$ augmente. Remarque que l'erreur d'approximation dépend de $\lambda$ plutôt que 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 garder les deux clairement distincts.)

Ce tutoriel montre comment générer un ensemble de tels circuits randomisés. Une fois ces circuits créés, de façon similaire à 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 temporels. Cela garantit un meilleur recouvrement entre les vecteurs d'état fondamental et les chaînes de bits échantillonnées.

Exigences (C++)

Cette implémentation s'exécute sur du matériel IBM Quantum® et utilise un exécutable SBD externe pour la diagonalisation classique. Elle n'inclut pas l'exemple de simulateur Python.

Suis les instructions de compilation C++ pour installer le compilateur, CMake, Rust, Boost et les dépendances de Qiskit. Installe Rust avec rustup : les dépendances de Qiskit épinglent leur propre version de Rust, que rustup télécharge automatiquement, mais un Rust installé depuis Homebrew ou un gestionnaire de paquets Linux peut être trop ancien. L'étape finale SBD nécessite en outre MPI, OpenMP et BLAS/LAPACK. Windows n'est actuellement pas pris en charge par cette implémentation.

Utilise le code source C++ complet et le projet CMake du dépôt. Conserve le répertoire adjacent fcidump_files lorsque tu télécharges l'exemple. Depuis une copie du dépôt de documentation, compile avec les commandes suivantes. La première compilation télécharge et compile toutes les dépendances et prend plusieurs minutes.

cd docs/tutorials/assets/sqdrift/cpp
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --config Release

La compilation écrit l'exécutable SqDRIFT dans docs/tutorials/assets/sqdrift/cpp (la racine du projet), et non dans build/. Exécute ./SqDRIFT depuis la racine du projet afin qu'il puisse trouver ../fcidump_files/N2_sto_3g. Toutes les autres commandes de cet onglet s'exécutent aussi depuis la racine du projet, sauf indication contraire. Si ./SqDRIFT échoue avec une erreur de bibliothèque partagée manquante, suis les instructions de chemin de bibliothèque propres à ta plateforme dans le guide de compilation ; sous macOS, l'exécutable fonctionne généralement sans elles.

Configuration (C++)

Les extraits ci-dessous sont tirés du programme SqDRIFT.cpp complet. Compile et exécute ce fichier comme un seul programme.

Avant d'exécuter SqDRIFT, assure-toi que ton jeton d'API IBM Quantum est configuré pour IBM Cloud Platform.

Si Python est installé avec qiskit-ibm-runtime, le moyen le plus simple de créer un fichier de configuration correctement formaté est :

from qiskit_ibm_runtime import QiskitRuntimeService
QiskitRuntimeService.save_account(channel="ibm_quantum_platform", token="YOUR_IBM_QUANTUM_TOKEN_HERE")

Cela écrit ~/.qiskit/qiskit-ibm.json (Linux/macOS) ou %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) avec la bonne structure. Si tu préfères créer le fichier manuellement :

{
  "default-ibm-quantum-platform": {
    "channel": "ibm_quantum_platform",
    "token": "YOUR_IBM_QUANTUM_TOKEN_HERE"
  }
}

Obtiens ta clé d'API sur IBM Quantum Platform.

Remarque : un fichier de configuration manquant ou mal formé provoque au démarrage une panique Rust unwrap() peu explicite au lieu d'une erreur claire. Vérifie bien le nom de la clé (default-ibm-quantum-platform) et la valeur du canal (ibm_quantum_platform) si tu rencontres un problème.

Implémentation (C++)

Étape 1 : associer les entrées classiques à un problème quantique (C++)

Charger l'hamiltonien moléculaire

char filename[] = "../fcidump_files/N2_sto_3g";
QfFCIDump* fcidump = qf_fcidump_from_file(filename);
if (fcidump == nullptr) {
    std::cerr << "Failed to load FCIDump!" << std::endl;
    return 1;
}

uint32_t norb = qf_fcidump_norb(fcidump);
uint32_t nelec = qf_fcidump_nelec(fcidump);
uint32_t n_alpha = nelec / 2;
uint32_t n_beta  = nelec - n_alpha;
uint32_t num_modes = 2 * norb;
uint32_t num_qubits = num_modes;
std::cout << "✓ Loaded N2 molecule (" << norb << " orbitals, "
          << num_qubits << " qubits, " << nelec << " electrons)" << std::endl;

Output:

✓ Loaded N2 molecule (10 orbitals, 20 qubits, 14 electrons)

Charge l'hamiltonien de la molécule $N_2$ à partir d'un fichier FCIDump. Le nombre de qubits est égal au double du nombre d'orbitales (spin haut et spin bas).

Créer l'hamiltonien et le mettre en ordre normal

// Create and normal order Hamiltonian
QfFermionOperator* hamiltonian = qf_ferm_op_from_fcidump(fcidump);
QfFermionOperator* normal_ordered = qf_ferm_op_normal_ordered(hamiltonian, nullptr);

Convertit les données FCIDump en opérateur fermionique et applique l'ordre normal pour simplifier la structure de l'opérateur.

Regrouper les termes par structure électronique

// Group terms by electronic structure
QfExitCode exit_code = qf_ferm_op_group_terms_by_electronic_structure(
    normal_ordered, num_modes, false
);
if (exit_code != QfExitCode_Success) {
    std::cerr << "Failed to group terms!" << std::endl;
    return 1;
}

uint32_t num_groups = qf_ferm_op_num_groups(normal_ordered);
std::cout << "Grouped into " << num_groups << " groups" << std::endl;

// Split into group operators
QfFermionOperator** group_ops = new QfFermionOperator*[num_groups];
qf_ferm_op_split_out_groups(normal_ordered, nullptr, 0, group_ops);

Output:

Grouped into 1902 groups

Regroupe les termes liés de l'hamiltonien par structure électronique avant de construire leurs circuits d'évolution.

Calculer les poids d'échantillonnage, puis normaliser et associer à des opérateurs de qubits

L'ordre compte ici. La normalisation remet chaque coefficient à une magnitude unitaire sur place, donc les poids d'échantillonnage doivent être calculés d'abord. Si la normalisation s'exécute en premier, chaque $|c_j|$ vaut exactement 1.0 et le poids de chaque groupe se réduit à un simple décompte de ses termes — la distribution d'échantillonnage serait alors pilotée par la taille des groupes plutôt que par les coefficients physiques.

Le poids de chaque groupe est le coefficient absolu moyen de ses termes. La passe QDriftTrotterization de l'onglet Python utilise le même poids pour un hamiltonien groupé. Utiliser la somme à la place ferait croître $\lambda$ avec la taille du groupe, ce qui fausserait à son tour le temps d'évolution $t_k = \lambda\tau/n$.

// Calculate sampling weights BEFORE normalization.
std::vector<double> weights(num_groups);
double total_weight = 0.0;

for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    double abs_sum = 0.0;
    for (uint64_t j = 0; j < num_terms; j++) {
        abs_sum += std::sqrt(coeffs[j].re * coeffs[j].re +
                             coeffs[j].im * coeffs[j].im);
    }
    weights[i] = (num_terms > 0) ? (abs_sum / static_cast<double>(num_terms)) : 0.0;
    total_weight += weights[i];
}

std::cout << "Total weight (λ): " << total_weight << std::endl;

// Normalize the evolution operators term-by-term, now that the sampling
// weights above have already captured the original magnitudes.
for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    for (uint64_t j = 0; j < num_terms; j++) {
        const double magnitude = std::sqrt(coeffs[j].re * coeffs[j].re +
                                           coeffs[j].im * coeffs[j].im);
        if (magnitude > 0.0) {
            coeffs[j].re /= magnitude;
            coeffs[j].im /= magnitude;
        }
    }
}

QkObs** qubit_ops = new QkObs*[num_groups];
for (uint32_t i = 0; i < num_groups; i++) {
    QfExitCode jw_exit = qf_ferm_op_jordan_wigner(group_ops[i], num_qubits, &qubit_ops[i]);
    if (jw_exit != QfExitCode_Success) {
        std::cerr << "Failed to map group " << i << std::endl;
        return 1;
    }
}
std::cout << "Mapped all " << num_groups << " normalized groups to qubit operators" << std::endl;

Output:

Total weight (λ): 337.261
Mapped all 1902 normalized groups to qubit operators

Le poids $\lambda_i = \frac{1}{|G_i|}\sum_{j \in G_i} |c_{ij}|$ est le coefficient absolu moyen sur les termes du groupe $G_i$, lu à partir des coefficients d'origine. La normalisation remet ensuite chaque terme fermionique à une magnitude unitaire avant la transformation de Jordan-Wigner, ne conservant que sa phase/son signe dans l'opérateur évolué. Cela correspond à l'évolution à termes normalisés de qDRIFT : les magnitudes des coefficients déterminent les probabilités d'échantillonnage via $\lambda$, tandis que l'évolution du circuit utilise le terme normalisé afin que les grands coefficients ne soient pas comptés deux fois.

Échantillonner des ensembles d'opérateurs

// SqDRIFT Sampling: Create operator sets
    const int num_circuits = 100;    // Number of circuits to create
    const int ops_per_circuit = 10;  // Operators per circuit
    const double time_step = 1;    // Time step for evolution
    
    std::cout << "\n Generating " << num_circuits << " circuits with " 
              << ops_per_circuit << " operators each..." << std::endl;
    
    std::mt19937 gen(42);
    std::discrete_distribution<> dist(weights.begin(), weights.end());
    
    // Store sampled operator indices
    std::vector<std::vector<int>> operator_sets(num_circuits);
    for (int i = 0; i < num_circuits; i++) {
        operator_sets[i].resize(ops_per_circuit);
        for (int j = 0; j < ops_per_circuit; j++) {
            operator_sets[i][j] = dist(gen);
        }
    }

Output:

 Generating 100 circuits with 10 operators each...

Tire les ensembles d'opérateurs pour 100 circuits aléatoires, chacun contenant 10 groupes échantillonnés dans la distribution pondérée. L'étape suivante transforme ces ensembles en l'ensemble de circuits qDRIFT.

Étape 2 : construire les circuits qDRIFT (C++)

Dans la version C++, cette étape construit les circuits à partir des ensembles d'opérateurs échantillonnés. La transpilation pour le Backend cible a lieu à l'étape 3, juste avant la soumission de chaque circuit.

    // Create Suzuki-Trotter circuits by composing all operators
    std::vector<Qiskit::circuit::QuantumCircuit> circuits;
    circuits.reserve(num_circuits);

    const uint32_t trotter_order = 1;
    const uint32_t trotter_reps = 1;
    const bool preserve_order = false;
    const bool insert_barriers = false;

    // Gate name mapping is static — compute once before all loops
    auto name_map = Qiskit::circuit::get_standard_gate_name_mapping();

    for (int circuit_idx = 0; circuit_idx < num_circuits; circuit_idx++) {
        // Create initial circuit with state preparation using C++ API
        Qiskit::circuit::QuantumCircuit qc(num_qubits, num_qubits);

        // State preparation: apply X gates to spin-up qubits 0..n_alpha-1
        // and spin-down qubits norb..norb+n_beta-1 (derived from FCIDump NELEC)
        for (uint32_t q = 0; q < n_alpha && q < num_qubits; q++) {
            qc.x(q);
        }
        for (uint32_t q = norb; q < norb + n_beta && q < num_qubits; q++) {
            qc.x(q);
        }

        // For each operator in this circuit's set
        for (int op_idx = 0; op_idx < ops_per_circuit; op_idx++) {
            int group_idx = operator_sets[circuit_idx][op_idx];
            double evolution_time = (total_weight * time_step) / ops_per_circuit;
            
            // Create Suzuki-Trotter evolution circuit
            QkCircuit* evolution_circuit = qk_circuit_library_suzuki_trotter(
                qubit_ops[group_idx], trotter_order, trotter_reps,
                evolution_time, preserve_order, insert_barriers
            );
            
            if (evolution_circuit != nullptr) {
                // Get number of instructions
                uint32_t num_ops = qk_circuit_num_instructions(evolution_circuit);
                
                // Manually append each instruction
                for (uint32_t i = 0; i < num_ops; i++) {
                    QkCircuitInstruction inst;
                    qk_circuit_get_instruction(evolution_circuit, i, &inst);

                    // Prepare qubit vector
                    std::vector<uint32_t> qubits(inst.num_qubits);
                    for (uint32_t j = 0; j < inst.num_qubits; j++) {
                        qubits[j] = inst.qubits[j];
                    }

                    // Get operation kind
                    QkOperationKind kind = qk_circuit_instruction_kind(evolution_circuit, i);

                    // Get mutable pointer to circuit
                    QkCircuit* mutable_circuit = qc.get_rust_circuit().get();

                    // Append based on operation type
                    if (kind == QkOperationKind_Gate) {
                        // Convert gate name string to QkGate enum
                        std::string gate_name(inst.name);
                        auto gate_it = name_map.find(gate_name);
                        if (gate_it != name_map.end()) {
                            QkGate gate_enum = gate_it->second.gate_map();
                            qk_circuit_parameterized_gate(
                                mutable_circuit,
                                gate_enum,           // Use QkGate enum, not string
                                qubits.data(),
                                inst.params
                            );
                        } else {
                            std::cerr << "Unknown gate in evolution circuit: " << gate_name << std::endl;
                            qk_circuit_instruction_clear(&inst);
                            qk_circuit_free(evolution_circuit);
                            return 1;
                        }
                    } else if (kind == QkOperationKind_Barrier) {
                        qk_circuit_barrier(
                            mutable_circuit,
                            qubits.data(), 
                            inst.num_qubits
                        );
                    } else if (kind == QkOperationKind_Reset) {
                        qk_circuit_reset(
                            mutable_circuit,
                            qubits[0]
                        );
                    }

                    qk_circuit_instruction_clear(&inst);
                }

                qk_circuit_free(evolution_circuit);
            }
        }

        circuits.push_back(qc);

        // Display progress
        if ((circuit_idx + 1) % 10 == 0) {
            std::cout << "  Created " << (circuit_idx + 1) << "/" 
                    << num_circuits << " circuits" << std::endl;
        }
    }

    std::cout << "✓ Created all " << num_circuits << " Suzuki-Trotter circuits" << std::endl;

Output:

  Created 10/100 circuits
  Created 20/100 circuits
  Created 30/100 circuits
  Created 40/100 circuits
  Created 50/100 circuits
  Created 60/100 circuits
  Created 70/100 circuits
  Created 80/100 circuits
  Created 90/100 circuits
  Created 100/100 circuits
✓ Created all 100 Suzuki-Trotter circuits

Le programme complet affiche ensuite des statistiques sur les circuits et les indices d'opérateurs du circuit 0, qui sont omis ici.

Pour chaque opérateur échantillonné, crée un circuit d'évolution de Suzuki-Trotter avec le temps $t_k = \frac{\lambda \tau}{n}$, où $n$ est le nombre d'opérateurs échantillonnés par circuit et $\tau$ est le pas de temps. Comme chaque groupe échantillonné a été normalisé terme par terme avant l'association, cette évolution n'applique que la phase/le signe de chaque terme pendant la rotation ; les magnitudes ne contribuent que par la distribution d'échantillonnage. L'état initial de Hartree-Fock est préparé en appliquant des portes X aux qubits 0-6 et 10-16 (sept électrons dans chaque secteur de spin). Les instructions de chaque circuit d'évolution sont ajoutées manuellement au circuit principal à l'aide de l'API C++.

Étape 3 : Exécuter à l'aide des primitives Qiskit (C++)

std::cout << "\n Connecting to IBM Quantum Cloud..." << std::endl;

// Initialize IBM Quantum Compute Service using qiskit-cpp
Qiskit::service::QiskitRuntimeService service;
std::cout << "✓ Connected to IBM Quantum" << std::endl;

// Get backend
const std::string backend_name = "ibm_fez";
auto backend = service.backend(backend_name);

if (backend.name().empty()) {
    std::cerr << "Backend " << backend_name
              << " is not available to this account." << std::endl;
    return 1;
}

std::cout << "✓ Selected backend: " << backend.name() << std::endl;

// Create sampler primitive
const int32_t shots = 100;  // Number of shots per circuit
Qiskit::primitives::BackendSamplerV2 sampler(backend, shots);

std::cout << "\n Submitting " << num_circuits << " circuits ("
            << shots << " shots each)..." << std::endl;

// Submit one job per circuit.
std::vector<std::shared_ptr<Qiskit::primitives::BasePrimitiveJob>> jobs;

for (int i = 0; i < num_circuits; i++) {
    for (uint32_t q = 0; q < num_qubits; q++) {
        circuits[i].measure(q, q);
    }
    Qiskit::circuit::QuantumCircuit transpiled_qc =
        Qiskit::compiler::transpile(circuits[i], backend, 2, 1.0, 42);
    std::vector<Qiskit::primitives::SamplerPub> pubs;
    pubs.push_back(Qiskit::primitives::SamplerPub(transpiled_qc, shots));
    jobs.push_back(sampler.run(pubs));
}

std::cout << "\n⏳ Waiting for " << num_circuits << " jobs to complete..." << std::endl;

// Collect results — poll each job and extract bitstrings
std::vector<boost::dynamic_bitset<>> all_bitstrings;

for (int i = 0; i < num_circuits; i++) {
    while (!jobs[i]->in_final_state()) {
        sleep(5);
    }
    if (jobs[i]->done()) {
        auto result = jobs[i]->result();
        if (result.size() > 0) {
            auto& pub_result = result[0];
            auto& bit_array  = pub_result.data();
            for (size_t s = 0; s < bit_array.num_shots(); s++) {
                auto sample = bit_array[s].to_string();
                boost::dynamic_bitset<> bs(num_qubits);
                for (uint32_t bit = 0; bit < num_qubits; bit++) {
                    bs[bit] = (sample[num_qubits - 1 - bit] == '1');
                }
                all_bitstrings.push_back(bs);
            }
        }
    } else {
        std::cout << "  Circuit " << (i+1) << " did not complete successfully" << std::endl;
    }
}

std::cout << "\n✓ Collected " << all_bitstrings.size()
            << " total bitstrings from all circuits" << std::endl;

Se connecte à IBM Quantum, sélectionne le Backend par son nom, transpile chaque circuit pour ce Backend et soumet un job par circuit : 100 jobs Sampler distincts de 100 shots chacun sur un registre de 20 qubits. Il interroge ensuite chaque job jusqu'à sa fin et collecte les bitstrings depuis result[0]. Chaque job consomme du temps QPU et attend séparément dans la file, donc attends-toi à ce que cette étape prenne beaucoup plus de temps qu'un seul job groupé.

Les résultats sur matériel cités sur cette page proviennent d'une exécution sur ibm_pittsburgh (100 circuits × 100 shots), avec backend_name modifié par rapport à la valeur par défaut ibm_fez. Les échantillons enregistrés ont été post-traités avec l'ordre des bits des déterminants corrigé, montré à l'étape 4. Le bruit matériel varie d'une exécution à l'autre, donc tes comptages et ton énergie finale seront différents.

Output:

 Connecting to IBM Quantum Cloud...
✓ Connected to IBM Quantum
✓ Selected backend: ibm_pittsburgh

 Submitting 100 circuits (100 shots each)...

⏳ Waiting for 100 jobs to complete...

✓ Collected 10000 total bitstrings from all circuits

Changer de Backend : modifie backend_name dans SqDRIFT.cpp pour cibler un autre appareil ; service.backends() liste ceux auxquels ton compte a accès.

Étape 4 : Post-traiter et retourner le résultat dans le format classique souhaité (C++)

Post-sélectionner les bitstrings et écrire les chaînes CI

// Postselect bitstrings by Hamming weight (derived from FCIDump NELEC)
std::vector<double> bitstring_weights(all_bitstrings.size(), 1.0);
auto [filtered_bitstrings, filtered_weights] = Qiskit::addon::sqd::postselect_bitstrings(
    all_bitstrings,
    bitstring_weights,
    Qiskit::addon::sqd::MatchesRightLeftHamming<uint32_t>(n_alpha, n_beta)
);
std::ignore = filtered_weights;

std::cout << "✓ Postselected " << filtered_bitstrings.size()
          << " bitstrings with Hamming weight ("
          << n_alpha << "," << n_beta << ")" << std::endl;

// Convert bitstrings to CI strings using SQD addon
auto ci_strings = Qiskit::addon::sqd::bitstrings_to_ci_strings_symmetrize_spin(
    filtered_bitstrings,
    std::nullopt  // No dimension limit
);

std::cout << "✓ Generated " << ci_strings.size() << " CI strings" << std::endl;

// Write CI strings to file for SBD
std::ofstream alpha_file("alphadets_from_sqd.txt");
for (const auto& ci_string : ci_strings) {
    std::string bitstr;
    // SBD expects the highest orbital index first (bit zero on the right).
    boost::to_string(ci_string, bitstr);
    alpha_file << bitstr << "\n";
}
alpha_file.close();

std::cout << "✓ Wrote " << ci_strings.size()
          << " CI strings to alphadets_from_sqd.txt" << std::endl;
std::cout << " Ready for SBD diagonalization!" << std::endl;

Output:

✓ Postselected 1137 bitstrings with Hamming weight (7,7)
✓ Generated 49 CI strings
✓ Wrote 49 CI strings to alphadets_from_sqd.txt
 Ready for SBD diagonalization!

Post-sélectionne les bitstrings ayant le bon poids de Hamming (n_alpha électrons de spin haut, n_beta électrons de spin bas, dérivés du champ NELEC de FCIDump), les convertit en chaînes d'interaction de configurations (CI) par symétrisation de spin, et les écrit dans alphadets_from_sqd.txt pour la diagonalisation en base sélectionnée qui suit.

Libérer les ressources

// Cleanup
// circuits vector will be automatically cleaned up (RAII)
// Free qubit operators
for (uint32_t i = 0; i < num_groups; i++) {
    qk_obs_free(qubit_ops[i]);
}
delete[] qubit_ops;
// Free fermionic operators
for (uint32_t i = 0; i < num_groups; i++) {
    qf_ferm_op_free(group_ops[i]);
}
delete[] group_ops;
qf_ferm_op_free(normal_ordered);
qf_ferm_op_free(hamiltonian);
qf_fcidump_free(fcidump);
return 0;

Libère correctement toute la mémoire et toutes les ressources allouées. Le std::vector<QuantumCircuit> C++ est nettoyé automatiquement via RAII. Libère manuellement les opérateurs de qubits, les opérateurs fermioniques et les données FCIDump.

Maintenant que nous avons créé la base sur laquelle projeter notre hamiltonien, nous pouvons passer au processus de diagonalisation pour obtenir l'estimation de l'état fondamental. Exécute les commandes suivantes depuis la racine du projet, docs/tutorials/assets/sqdrift/cpp.

Compiler le binaire SBD diag

CMake récupère SBD en source uniquement dans deps/sbd/ (la compilation principale évite volontairement de le compiler). L'application SBD utilise un simple Makefile, et SBD lui-même nécessite MPI, OpenMP et BLAS/LAPACK — donc un simple make ne suffit pas tant que ceux-ci ne sont pas installés et que le fichier Configuration de l'application ne correspond pas à ta chaîne d'outils.

Le Makefile lit ses options de compilation et d'édition de liens dans le fichier Configuration du même répertoire. La version fournie en amont est réglée pour une disposition macOS/Homebrew précise et ne se compilera pas telle quelle sur une machine standard — notamment, -I/opt/homebrew/opt/llvm/include ne contient pas omp.h (il se trouve dans un répertoire de ressources clang propre à la version), et le clang++ d'Apple rejette -fopenmp d'emblée :

clang++: error: unsupported option '-fopenmp'

La recette suivante est vérifiée et fonctionnelle sur macOS (Apple Silicon, Darwin 25.6, Open MPI 5.x + Homebrew LLVM).

macOS (Apple Silicon) — testé

# 1. Install MPI and an OpenMP-capable compiler.
#    Apple clang cannot do OpenMP, so Homebrew's LLVM provides both
#    clang++ and libomp. BLAS/LAPACK come from Apple's Accelerate
#    framework, which is already part of macOS — no openblas needed.
brew install open-mpi llvm

Remplace ensuite deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization/Configuration par :

# Path to the SBD library
SBD_PATH=../..

# MPI C++ compiler. OMPI_CXX points Open MPI's wrapper at Homebrew's clang++,
# which (unlike Apple clang) supports -fopenmp.
CCCOM=OMPI_CXX=/opt/homebrew/opt/llvm/bin/clang++ mpicxx

# Build flags. omp.h lives in a version-specific clang resource dir, so let
# Homebrew's clang++ find it itself rather than hardcoding an -I path.
CCFLAGS= -std=c++17 -fopenmp -O3

# Link flags: OpenMP runtime from Homebrew LLVM, BLAS/LAPACK from the
# Accelerate framework that ships with macOS (no openblas install needed).
SYSLIB= -L/opt/homebrew/opt/llvm/lib -lomp -framework Accelerate

# 2. Build
cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization
make
cd ../../../..

Une compilation réussie affiche les deux lignes de compilation/édition de liens et laisse un binaire diag à côté du Makefile. Un message ld: warning: ignoring duplicate libraries: '-lomp' est sans danger.

macOS Intel : le préfixe de Homebrew est /usr/local et non /opt/homebrew, donc remplace-le dans les deux chemins précédents.

Linux (Ubuntu/Debian) — non testé

Les prérequis suivants couvrent MPI, OpenMP et BLAS/LAPACK. GCC prend en charge -fopenmp nativement, donc aucun remplacement de compilateur n'est nécessaire. Cela n'a pas été vérifié sur Linux :

sudo apt install -y libopenmpi-dev libomp-dev libblas-dev liblapack-dev

SBD_PATH=../..
CCCOM=mpicxx
CCFLAGS= -std=c++17 -fopenmp -O3
SYSLIB= -llapack -lblas

Exécuter SqDRIFT pour générer les chaînes CI

Si tu as déjà exécuté ./SqDRIFT pendant les étapes 1 à 4, alphadets_from_sqd.txt existe déjà à la racine du projet et tu peux sauter cette étape. Sinon, exécute-le depuis la racine du projet :

./SqDRIFT
# Writes alphadets_from_sqd.txt in the project root

Exécuter la diagonalisation

cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization

# The FCIDump lives in the shared sqdrift/fcidump_files/ directory, one level
# above this tutorial's project root, hence five levels up from here.
ln -sf ../../../../../fcidump_files/N2_sto_3g fcidump.txt

# alphadets_from_sqd.txt is written by ./SqDRIFT into the project root itself,
# which is four levels up.
ln -sf ../../../../alphadets_from_sqd.txt alphadets.txt

./diag \
--fcidump fcidump.txt \
--adetfile alphadets.txt \
--method 0 \
--iteration 100 \
--block 10 \
--tolerance 1e-8 \
| tee ../../../../sbd_output.txt

cd ../../../..

Paramètres :

• --fcidump : chemin vers le fichier FCIDump contenant l'hamiltonien
• --adetfile : chemin vers le fichier de déterminants alpha (chaînes CI issues de SqDRIFT)
• --method 0 : méthode de diagonalisation de Davidson
• --iteration 100 : nombre maximal d'itérations de Davidson
• --block 10 : taille de bloc pour l'algorithme de Davidson
• --tolerance 1e-8 : tolérance de convergence sur l'énergie

 Elapsed time for helper construction 0.001512 (sec)
 Elapsed time for init 1e-06 (sec)
 Davidson iteration 0.0 (tol=0.1857493937033264): -107.493531425221
 Davidson iteration 0.1 (tol=0.0455699976495391): -107.5098150493379 -105.93368397386
 Davidson iteration 0.2 (tol=0.009291923190816917): -107.5108356028703 -106.1245977311231 -105.6823366883346
 Davidson iteration 0.3 (tol=0.001488674426419356): -107.5108711084406 -106.3219246016818 -105.7376912083761 -105.4392974334491
 Davidson iteration 0.4 (tol=0.0003188697654737439): -107.5108722237947 -106.5935152393201 -106.1108017566682 -105.4932621297472
 Davidson iteration 0.5 (tol=3.080273116687244e-05): -107.5108722587463 -106.7245930484005 -106.2313070085243 -105.6423237790271
 Davidson iteration 0.6 (tol=2.981303155550973e-06): -107.5108722590151 -106.8008116251279 -106.330411722324 -105.69489630447
 Davidson iteration 0.7 (tol=4.565075388900609e-07): -107.5108722590186 -106.8090569201597 -106.3345525819528 -106.0578325555555
 Davidson iteration 0.8 (tol=4.870689101320012e-08): -107.5108722590186 -106.8098939189322 -106.3973756614146 -106.1617813840531
 Davidson iteration 0.9 (tol=6.224330513865847e-09): -107.5108722590187 -106.8159869231416 -106.4315891597535 -106.2222021964013
 Elapsed time for davidson 0.031276 (sec)
 Elapsed time for diagonalization 0.031286 (sec)
 Elapsed time for mult 0.002511 (sec)
 Energy = -107.5108722590187
 Elapsed time for measurement 4.6e-05 (sec)
 Sample-based diagonalization: Energy = -107.5108722590187
 Sample-based diagonalization: density = [1.999989340109303,1.99999662423008,1.999199121034164,1.981735389885314,1.999825046204946,1.999847567448985,1.995347419300126,0.01136512984757364,0.009844780833945199,0.002849581105572152
 Sample-based diagonalization: carryover bitstrings = [], size = 0

La ligne Energy = -107.5108722590187 est l'estimation de l'état fondamental sur le sous-espace échantillonné.

En quoi cela diffère de l'exemple Python

Les deux exemples diffèrent à plusieurs égards à la fois, et pas seulement par leur traitement des bitstrings de mauvaise symétrie, donc ils ne sont pas censés produire des énergies identiques :

	Onglet C++	Onglet Python
Bitstrings de mauvaise symétrie	Écartés par postselect_bitstrings sur le poids de Hamming	Réparés par la récupération itérative de configurations
Structure de boucle	Un seul passage : échantillonner → post-sélectionner → diagonaliser une fois	Boucle externe : diagonaliser, lire les occupations orbitalaires, récupérer les configurations, rediagonaliser
Nombre de circuits	100 circuits (num_circuits = 100)	400 circuits (200 par temps d'évolution)
Temps d'évolution	Un seul $\tau$ (time_step = 1)	Deux temps, times = [1.0, 10.0]
Termes diagonaux	Conservés dans l'opérateur échantillonné	Retirés par filter_diagonal_terms, ce qui modifie la distribution d'échantillonnage et les circuits
Ordre des termes	Groupés uniquement	canonical_order appliqué avant le regroupement pour la reproductibilité
Solveur de diagonalisation	SBD (diag, binaire externe MPI/OpenMP)	Solveur qiskit-addon-sqd en processus
Backend	Fixé par son nom (backend_name = "ibm_fez" ; l'exécution de cette page a utilisé ibm_pittsburgh)	Backend opérationnel least_busy avec assez de qubits
Transpilation	Chaque circuit séparément, niveau d'optimisation 2	Tous les circuits ensemble, niveau d'optimisation 3
Soumission des jobs	100 jobs, un circuit chacun	Un job Sampler contenant les 400 circuits
Exécution sur simulateur	Aucune ; matériel uniquement	Exemple de simulateur sans bruit avant l'exécution sur matériel

Les deux versions utilisent le même hamiltonien, regroupent ses termes par structure électronique, pondèrent chaque groupe par le coefficient absolu moyen de ses termes, partent de l'état de Hartree-Fock et échantillonnent 10 groupes et 100 shots par circuit.

Chaque différence affecte indépendamment le sous-espace échantillonné, et donc l'énergie. L'ensemble plus grand et le second temps d'évolution, plus long, élargissent tous deux les configurations explorées par l'exemple Python ; échantillonner à plus d'un temps donne à SqDRIFT sa couverture de type Krylov de l'espace de basse énergie. Filtrer les termes diagonaux change les opérateurs que l'échantillonneur peut tirer et leurs poids, donc les deux exécutions échantillonnent des distributions différentes.

La post-sélection est un filtre strict, et un shot avec un mauvais nombre de particules est tout simplement abandonné. Dans l'exécution ci-dessus, cela a écarté la plupart des données : 1137 shots sur 10000 ont passé la vérification du poids de Hamming (7,7) (environ 11 %), se réduisant à seulement 49 $\alpha$-déterminants distincts. La récupération de configurations utilise au contraire les occupations orbitalaires moyennes d'une diagonalisation précédente pour inverser des bits et réparer ces shots, recyclant ce que la post-sélection écarterait. Cette rétroaction la rend intrinsèquement itérative : recover_configurations dans qiskit-addon-sqd-hpc requiert un argument avg_occupancies qui n'existe qu'après l'exécution d'une diagonalisation.

Cette exécution C++ donne un sous-espace de 49 $\alpha$-déterminants et une énergie de -107.5109 Ha. L'espace complet des $\alpha$-déterminants pour sept électrons dans 10 orbitales ne compte que $\binom{10}{7} = 120$ chaînes, donc fournir les 120 au même binaire diag donne l'énergie FCI exacte -107.6482 Ha. L'exécution se situe donc environ 0.1373 Ha au-dessus de la valeur exacte tout en couvrant 49 des 120 $\alpha$-déterminants. Comme SQD est variationnel dans le sous-espace sélectionné, un sous-espace plus petit ne peut qu'augmenter l'énergie, donc le sous-espace incomplet est l'explication la plus directe de l'écart.

Cette vérification par FCI exacte n'est possible que parce que l'exemple est petit ; elle ne fait pas partie du workflow SqDRIFT, et pour les tailles de système où SQD est réellement nécessaire, aucune référence de ce type n'existe. Elle sert ici uniquement à quantifier la part de l'espace récupérée par cette exécution.

Une comparaison avec l'exemple Python n'est pas contrôlée, donc une différence ne doit pas être attribuée à la seule absence de récupération de configurations ; note que le nombre de circuits, le temps d'évolution unique, le filtrage des termes diagonaux et le bruit matériel diffèrent tous à la fois et façonnent tous le sous-espace.

Pour combler l'écart, voici quelques pistes à essayer après avoir terminé ce tutoriel : augmenter shots, num_circuits ou ops_per_circuit ; échantillonner à plus d'un time_step ; et ajouter une boucle de récupération au-dessus de Qiskit::addon::sqd::recover_configurations.

Prochaines étapes

Recommandations

• Diagonalisation quantique de Krylov basée sur des échantillons d'un modèle de réseau fermionique — un tutoriel connexe utilisant des circuits d'évolution temporelle au lieu d'un ansatz variationnel.

• Diagonalisation quantique basée sur des échantillons 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 littérature sur laquelle repose ce tutoriel. (Note que certaines des optimisations discutées dans cet article sont actuellement en cours de développement, et que ce tutoriel est susceptible de changer à l'avenir selon l'évolution des bibliothèques utilisées.)