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

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

Tu cherches la version C++ ?

Ce tutoriel utilise Python. Pour l'implémentation en C++, y compris le code source et les instructions de construction, consulte le tutoriel SqDRIFT en C++.

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

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

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.

Ce tutoriel utilise Qiskit Fermions pour créer les circuits fermioniques plus naturels pour l'algorithme qDRIFT, suivi de l'utilisation de passes de disposition (layout) et de synthèse fermioniques avant de brancher les circuits dans le pipeline Qiskit traditionnel pour l'exécution matérielle.

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

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

# 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 : associer les entrées classiques à un problème quantique

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 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 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 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. 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 = "assets/sqdrift/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 = "assets/sqdrift/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 de l'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 l'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 dont échantillonne le protocole qDRIFT, 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$$

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 de terminer le post-traitement 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 comportant 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, 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 ; le compagnon 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 = "assets/sqdrift/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 t'intéresser :

• Diagonalisation quantique de Krylov basée sur des échantillons pour un modèle de réseau fermionique - un tutoriel connexe utilisant des circuits d'évolution temporelle plutôt qu'un ansatz variationnel.
• Diagonalisation quantique basée sur des échantillons d'un hamiltonien de chimie - un tutoriel sur la construction d'un circuit de cluster de Jastrow unitaire local (LUCJ) pour la simulation en chimie quantique.
• L'article SqDRIFT - la référence sur laquelle ce tutoriel est basé. (Note que certaines des optimisations abordées dans cet article sont actuellement en cours de développement, et ce tutoriel est susceptible d'évoluer à l'avenir en fonction de l'évolution des bibliothèques utilisées.)