Algoritmo SqDRIFT per la stima dello stato fondamentale

Python

Scarica il notebook Python

Scarica il notebook Python per eseguire tu stesso questo tutorial. Il notebook scarica l'input FCIDump se non è già presente.

Stima di utilizzo: 180 secondi su un processore Heron r3 (NOTA: questa è solo una stima. I tempi di esecuzione effettivi potrebbero variare.)

Risultati di apprendimento

• Scopri come creare circuiti a profondità minore rispetto alla Trotterizzazione

• Percorri un workflow end-to-end per la stima dello stato fondamentale usando qDRIFT e SQD

• Scopri come usare qiskit-fermions insieme ad altri addon Qiskit per implementare tale workflow

Prerequisiti

• Leggi la panoramica Diagonalizzazione quantistica basata su campioni (SQD)

• Leggi la lezione Diagonalizzazione quantistica di Krylov basata su campioni (SKQD)

Contesto

SqDRIFT è una variante di SKQD che sostituisce la necessità di scegliere un ansatz da cui campionare bitstring con un insieme di circuiti di evoluzione temporale costruiti direttamente dalla Hamiltoniana target. Questo si ottiene sottocampionando operatori di evoluzione temporale più piccoli dalla Hamiltoniana in base ai suoi coefficienti, il che è noto come metodo di Trotterizzazione qDRIFT.

Entrambe le implementazioni usano Qiskit Fermions per rappresentare la Hamiltoniana molecolare e costruire i circuiti per l'algoritmo qDRIFT. L'implementazione Python usa pass del transpiler fermionici, mentre l'implementazione C++ campiona esplicitamente i gruppi e li mappa su operatori di qubit.

Sia la Hamiltoniana della forma:

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

dove, senza perdita di generalità, richiediamo $c_i > 0$ e che l'autovalore più grande di $h_i$ sia uguale, in valore assoluto, a $1$. Qualsiasi prefattore con segno o complesso viene assorbito in $h_i$, quindi i coefficienti $c_i$ sono pesi strettamente positivi mentre gli $h_i$ portano la direzione di ogni termine. Qui $N$ è il numero di termini (o, dopo il raggruppamento, il numero di gruppi) nella Hamiltoniana; è una proprietà della Hamiltoniana ed è distinto dal numero di operatori campionati in un singolo circuito, indicato con $n$ di seguito.

L'algoritmo qDRIFT realizza quindi, per il tempo target $t$, un operatore $V_k$, dove $k$ va da $1 \cdots K$ e indica il $k$-esimo circuito SqDRIFT, definito come:

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

Qui $n$ è il numero di operatori campionati per circuito e $K$ è il numero di circuiti nell'insieme. Il prodotto si estende sulle $n$ estrazioni, non su tutti gli $N$ termini della Hamiltoniana, e poiché i termini sono estratti con reinserimento, lo stesso $h_i$ può comparire più di una volta in un singolo $V_k$.

La quantità:

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

è la norma $L_1$ dei coefficienti, quindi ciascuno degli $n$ passi evolve per la stessa durata $\lambda t / n$ indipendentemente da quale termine è stato estratto. L'uniformità dell'angolo del passo è la caratteristica distintiva di qDRIFT: un coefficiente influenza il risultato tramite quanto spesso il suo termine viene estratto, non tramite quanto quel termine viene ruotato. Gli indici sono campionati dalla distribuzione:

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

quindi la serie $(k_1, \ldots, k_n)$ è una sequenza casuale di indici di termini estratti da questa distribuzione. Poiché i $c_i$ sono positivi e sommano a $\lambda$, questa è una distribuzione di probabilità normalizzata, e il valore atteso del canale risultante sulle estrazioni casuali approssima l'evoluzione sotto $H$, con un errore che diminuisce al crescere di $n$. Nota che l'errore di approssimazione dipende da $\lambda$ piuttosto che dal numero di termini $N$.

(L'articolo SqDRIFT scrive il numero di termini come $\mathcal{N}$ e la lunghezza della sequenza come $N$; qui usiamo $N$ e $n$ per mantenere i due chiaramente distinti.)

Questo tutorial mostra come generare un insieme di tali circuiti randomizzati. Dopo aver creato questi circuiti, in modo simile a come si crea un sottospazio di Krylov per operatori diversi, campioniamo bitstring da più di tali operatori con parametri temporali diversi. Questo assicura una maggiore sovrapposizione tra i vettori dello stato fondamentale e le bitstring campionate.

Requisiti

Prima di iniziare questo tutorial, assicurati di aver installato

• Un ambiente virtuale Python (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (nota che il nome è al plurale)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Puoi installare tutti i pacchetti richiesti con:

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

Configurazione

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

Implementazione

Passo 1: Mappa gli input classici a un problema quantistico

I passaggi Python riportati di seguito usano il simulatore Aer. L'esempio hardware segue più avanti.

Lettura e preparazione del FCIDump

Per questo tutorial, caricheremo la Hamiltoniana della struttura elettronica per l'azoto (N2). Esistono anche altri modi per creare operatori fermionici. Consulta la documentazione su qiskit_fermions.operators.library.

Informazioni su questo FCIDump. Il file N2_sto_3g descrive una molecola di azoto ($N_2$) nella base minimale STO-3G a una separazione interatomica di 1.09 $\AA$, la lunghezza di legame di equilibrio sperimentale. Il suo header dichiara NORB=10, NELEC=14 e MS2=0: 10 orbitali spaziali (quindi 20 orbitali di spin, e 20 qubit sotto Jordan-Wigner), 14 elettroni in un singoletto di spin, quindi sette elettroni $\alpha$ e sette $\beta$. A tutti gli orbitali è assegnata l'etichetta di simmetria 1, ovvero non viene sfruttata alcuna simmetria di gruppo puntuale. Trattandosi di un dump STO-3G a spazio completo, nessun orbitale è congelato e lo spazio di correlazione è abbastanza piccolo da permettere di calcolare classicamente un'energia di riferimento FCI esatta per il confronto, come mostrato nel blocco di codice successivo.

Un file equivalente può essere rigenerato con 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")

Poiché gli integrali dipendono dagli orbitali SCF a convergenza, un file rigenerato potrebbe differire da quello fornito nella fase o nell'ordinamento degli orbitali; le energie totali non ne sono influenzate.

Ottenere il file. Trova il FCIDump in questo repository GitHub. Puoi eseguire il codice qui sotto per scaricarlo nella posizione attesa dal resto del tutorial.

Per prima cosa usiamo il cisolver fornito da pyscf per ottenere l'energia di riferimento. Questa è la vera energia dello stato fondamentale della molecola con cui stiamo lavorando. Per farlo dichiareremo prima norb e nelec, che sono rispettivamente il numero di orbitali e il numero di elettroni. Poi dichiariamo h1e e h2e, che sono rispettivamente gli integrali a uno e due elettroni. Tutti questi verranno usati in seguito anche per 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

Caricamento dell'Hamiltoniana

Con i dati necessari pronti, leggiamo l'Hamiltoniana dal file FCI in un formato compatibile con qiskit-fermions

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

Workflow fermionici con qiskit-fermions

Per prima cosa mappiamo l'Hamiltoniana in un modello di circuito fermionico usando qiskit-fermions, che fornisce pass del transpiler e gate specifici per i circuiti fermionici. Questi verranno usati successivamente prima dei tradizionali pass del transpiler di Qiskit per questo workflow.

Raggruppamento dei termini

Per garantire la riproducibilità dei risultati, usiamo prima canonical_order per ordinare i termini basandoci solo sulla loro struttura. L'ordine degli operatori nella lista canon è quindi fisso. Questo garantisce la riproducibilità degli operatori creati perché il pass QDriftTrotterization che useremo in seguito campiona indici casuali per creare gli operatori qDRIFT.

In questo passaggio, sfruttiamo le numerose simmetrie presenti nell'Hamiltoniana della struttura elettronica raggruppando i termini correlati con coefficienti identici. Sebbene ciò cambi la distribuzione dei coefficienti dell'operatore da cui il protocollo qDRIFT campiona, questo non influisce sulle sue garanzie di convergenza. Fondamentalmente, raggruppare termini correlati per simmetria produce una cancellazione favorevole dei termini di Pauli e una profondità complessiva del circuito inferiore quando si evolve nel tempo uno stato sotto la loro azione.

qiskit-fermions fornisce la funzione group_terms_by_electronic_structure che esegue questo raggruppamento per noi.

Nota che group_terms_by_electronic_structure presuppone termini in ordinamento normale.

Filtraggio dei termini diagonali

Rimuoviamo i termini diagonali dall'Hamiltoniana usata per generare i circuiti, in modo che gli $n$ slot di campionamento qDRIFT vengano spesi su termini che spostano popolazione tra configurazioni. Tali termini è meglio filtrarli dall'Hamiltoniana a questo punto, prima che il gate Evolution venga costruito nel passaggio successivo.

I termini in questione sono quelli diagonali nella base dei numeri di occupazione, ovvero i prodotti di operatori numero $a^\dagger_i a_i$. Tre tipi di termini rientrano in questa descrizione:

• il termine costante di offset energetico, un prodotto di zero operatori numero, la cui evoluzione temporale contribuisce solo a una fase globale;

• i singoli operatori numero $n_i$, la cui evoluzione temporale si riduce a rotazioni $Z$ su singolo qubit;

• i prodotti di ordine superiore come $n_i n_j$.

Da soli, nessuno di questi sposta popolazione tra configurazioni del numero di occupazione; agiscono solo sulle fasi delle configurazioni già presenti. Non sono tuttavia inerti: quelle fasi relative alimentano l'interferenza generata dai termini di eccitazione più avanti nel circuito, quindi filtrarli cambia l'evoluzione effettivamente generata e può modificare la distribuzione di campionamento. Questa è un'approssimazione deliberata nel passaggio di generazione del circuito, fatta per concentrare il campionamento sui termini di eccitazione, piuttosto che un passaggio che lascia inalterata la distribuzione campionata. A differenza del raggruppamento per simmetria sopra descritto, che lascia intatte le garanzie di convergenza qDRIFT, questo filtro cambia l'operatore che viene evoluto. I circuiti quindi non approssimano più l'evoluzione sotto l'Hamiltoniana completa, e i limiti di errore qDRIFT si applicano all'operatore filtrato piuttosto che a quello originale. Questo è accettabile qui perché i circuiti sono solo un'euristica di campionamento usata per proporre configurazioni: nessun termine viene perso dalla stima dell'energia in sé, poiché il filtro si applica solo all'Hamiltoniana usata per costruire i circuiti, mentre la diagonalizzazione classica successiva usa l'Hamiltoniana completa, termini diagonali inclusi. L'accuratezza di SQD dipende da quel passaggio classico, che rimane variazionale nel sottospazio campionato indipendentemente da come sono state proposte le configurazioni.

La funzione filter_diagonal_terms() rimuove tali termini da un operatore in place. Li identifica dalla loro struttura in ordinamento normale — il multiset dei modi di creazione che corrisponde al multiset dei modi di annichilazione — quindi è valida solo su un operatore già in ordinamento normale. Questo presupposto non viene verificato a runtime.

# 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

Ora che abbiamo raggruppato i termini nell'Hamiltoniana, decideremo i seguenti parametri per generare l'insieme di circuiti:

• Il numero di circuiti da generare: num_circuits
• La lunghezza di ciascun circuito in termini di gruppi di eccitazione: num_exc
• Il fattore per i diversi tempi di evoluzione: times

Creazione dei circuiti fermionici

Ora creeremo circuiti fermionici per ciascuno dei passi temporali. Ogni circuito consisterà in un singolo gate di evoluzione, con il tempo di evoluzione dichiarato in precedenza. L'operatore di evoluzione è l'Hamiltoniana. In seguito eseguiremo pass del transpiler su questi circuiti per creare circuiti qDRIFT.

Preparazione dell'ansatz

Prepariamo lo stato di Hartree-Fock usando la classe InitializeModes. Per l'azoto, il processo consiste semplicemente nell'applicare gate X ai primi num_elec_a qubit e poi ai num_elec_b qubit, entrambi uguali a sette per l'azoto. Questo stato rappresenta i sette elettroni $\alpha$ e i sette $\beta$ dell'azoto.

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

Passo 2: Ottimizzare il problema per l'esecuzione su hardware quantistico

Ora che abbiamo i nostri circuiti, useremo prima i pass disponibili in qiskit-fermions per eseguire ottimizzazioni a livello fermionico, seguite dalla transpilazione del nostro circuito per il backend scelto. Poiché questo è un esperimento su simulatore, lo faremo prima per AerSimulator.

Calcolo del peso per ciascun gruppo

In questo passaggio, eseguiamo il campionamento qDRIFT dei termini stocasticamente con probabilità proporzionali ai loro coefficienti nell'Hamiltoniana. Il pass del transpiler qDRIFT fa questo per noi. Possiamo ora creare circuiti meno profondi che possono essere eseguiti sull'hardware in modo più efficiente nonostante la connettività limitata dei qubit, anche quando l'Hamiltoniana contiene accoppiamenti a lungo raggio e termini di ordine superiore al quadratico. Dopo il raggruppamento dei termini, campiona gli operatori in base ai loro pesi. Per ciascun operatore $h_i$, il peso $W_{h_i}$ è definito come segue:

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

Poiché i termini sono stati raggruppati nel Passo 1, ogni $h_i$ qui è un intero gruppo: $c_i$ è il coefficiente assoluto medio dei termini del gruppo $i$, e ogni termine del gruppo viene evoluto con il proprio coefficiente ridotto al solo segno.

Ottimizzazioni fermioniche e hardware-native

La funzione generate_preset_jw_pass_manager() restituisce un MultiStagePassManager che prende un FermionicCircuit e produce un circuito finale ottimizzato che possiamo compilare per eseguirlo sul nostro hardware. Sostituiamo la sua fase di ottimizzazione predefinita con un FermionicPassManager contenente il nostro pass QDriftTrotterization:

• Il pass QDriftTrotterization usa internamente il calcolo dei pesi e il campionamento per generare i circuiti che useremo per il campionamento

• Il pass RelabelModes è un altro pass di ottimizzazione che può essere usato per permutare i modi fermionici per ottimizzare la connettività tra i qubit e ridurre la profondità dei gate; leggi di più nel riferimento API

Le fasi rimanenti del MultiStagePassManager vengono eseguite automaticamente e gestiscono il mapping completo da fermione a qubit:

• F2QLayout: Il preset pass manager applica il pass TrivialF2QLayout, che mappa banalmente $n$ bit fermionici su $n$ qubit.

• F2QSynth: Un pass di transpilazione per mappare istruzioni di circuito basate su fermioni in istruzioni basate su qubit.

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

Ora che abbiamo terminato le ottimizzazioni a livello fermionico, possiamo compilare i circuiti per l'esecuzione sul simulatore.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Passo 3: Eseguire usando le primitive di Qiskit

Ora che abbiamo i nostri circuiti, possiamo eseguirli usando le primitive di Qiskit su AerSimulator. Combineremo tutti i conteggi dei diversi circuiti. Li convertiamo in vettori booleani prima di elaborarli infine con 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

Passo 4: Post-elaborare e restituire il risultato nel formato classico desiderato

Usare le bitstring per SQD

Possiamo ora eseguire lo schema di diagonalizzazione sulle bitstring selezionate per trovare l'autovalore più basso che corrisponderà all'energia dello stato fondamentale della molecola. Creiamo una funzione di callback, dichiariamo le occupazioni iniziali e impostiamo i parametri prima di eseguire infine lo schema di diagonalizzazione. La funzione di callback viene usata per stampare l'iterazione corrente e la stima corrente dell'autovalore a ogni iterazione.

Infine, per ottenere la stima dello stato fondamentale, aggiungiamo nuclear_repulsion_energy all'energia risultante.

Nota: la dimensione del sottospazio non è fissa tra le iterazioni, anche sul simulatore senza rumore — ogni sottocampione estrae un insieme diverso di configurazioni, e il passaggio di recupero rimodella il pool tra un'iterazione e l'altra, quindi la dimensione riportata varia da un sottocampione all'altro. Il campionamento senza rumore non fissa di per sé la dimensione del sottospazio selezionato. L'esecuzione su hardware, tuttavia, tende a dare sottospazi sistematicamente più grandi, perché gli shot rumorosi rompono la simmetria del numero di particelle e il recupero delle configurazioni li trasforma in vettori di base aggiuntivi. Per questo motivo, introdurremo anche un altro passaggio per la potatura delle bitstring nella sezione hardware.

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

Esempio hardware

Questo esempio usa 20 qubit (10 orbitali spaziali). Questa scelta è una comodità per un tutorial che deve eseguirsi rapidamente, non un limite rigido del metodo.

Il costo del passaggio classico non è determinato direttamente dal numero di qubit. SQD diagonalizza l'Hamiltoniana proiettata sul sottospazio generato dalle configurazioni campionate, quindi ciò che guida il costo classico è la dimensione di quel sottospazio selezionato — governata qui da samples_per_batch, num_batches e da quante configurazioni distinte i circuiti effettivamente producono — insieme all'algebra lineare sparsa necessaria per applicare l'Hamiltoniana proiettata. Lo spazio CI completo cresce combinatoriamente con orbitali ed elettroni, ma il sottospazio selezionato è una porzione piccola e regolabile di esso, e ne controlliamo direttamente la dimensione. Di conseguenza, il numero di qubit e la difficoltà classica si possono regolare in modo relativamente indipendente: uno spazio orbitale più ampio campionato in un sottospazio modesto può essere più economico di un sistema più piccolo diagonalizzato su uno molto grande.

In pratica, quindi, la dimensione del sistema fattibile dipende dalla dimensione del sottospazio necessaria per l'accuratezza desiderata e dalla memoria e dai core disponibili per l'eigensolver. Spazi orbitali più grandi in genere richiedono un sottospazio più ampio per raggiungere l'accuratezza chimica, ed è ciò che alla fine motiva l'uso di risorse distribuite — vedi qiskit-addon-sqd-hpc per scalare questo passaggio. Piuttosto che assumere un limite fisso, l'approccio pratico è osservare la dimensione del sottospazio riportata e la convergenza dell'energia tra le iterazioni e aumentare la dimensione del sottospazio finché l'energia smette di migliorare o si esaurisce la memoria disponibile.

Nota: a causa dell'errore di campionamento dovuto al rumore dell'hardware, il sottospazio creato per la diagonalizzazione nell'esecuzione su hardware sarà più grande di quello ottenuto usando il simulatore. Anche se aumenta la dimensione del sottospazio da diagonalizzare, il workflow fornisce comunque una risposta accurata grazie alla robustezza di SQD rispetto al rumore.

Potatura delle stringhe spurie

Qui possiamo scegliere di eseguire un passaggio aggiuntivo. Quando abbiamo tutte le bitstring dalle esecuzioni del circuito, possiamo filtrare le bitstring non valide prima di eseguire SQD, oppure procedere senza potatura. Saltare la potatura è generalmente preferibile per le esecuzioni su hardware, perché lascia gli shot con simmetria rotta disponibili al recupero delle configurazioni, che può ripararli in configurazioni valide e quindi ampliare il sottospazio invece di scartare del tutto quegli shot.

Poiché l'azoto può avere solo sette elettroni $\alpha$ e sette $\beta$, qualsiasi bitstring che non abbia esattamente sette 1 nella prima e nella seconda metà dell'output può essere scartata. Definiamo una funzione che verifica se le bitstring sono valide, e in caso contrario le scarta. Una volta filtrate le bitstring spurie, il resto viene inviato allo schema di diagonalizzazione. Usa il flag PRUNE qui sotto per passare da un comportamento all'altro.

Tieni presente che la potatura è solo una delle diverse scelte che determinano il sottospazio finale, insieme al numero di circuiti, all'insieme dei tempi di evoluzione e al filtraggio dei termini diagonali. Confrontare un'esecuzione con potatura rispetto a una senza è informativo solo se tutto il resto rimane fisso; la scheda C++ ne parla più in dettaglio, poiché effettua una postselezione anziché un recupero e differisce anche in quegli altri parametri.

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++

Ottieni il codice C++

Il sorgente C++, il progetto CMake e le istruzioni di build si trovano nella directory docs/tutorials/assets/sqdrift/cpp del repository della documentazione di Qiskit. L'input FCIDump si trova nella directory vicina fcidump_files. Clona il repository, oppure scarica entrambe le directory, per compilare ed eseguire tu stesso questo tutorial.

Stima di utilizzo: 180 secondi su un processore Heron r3 (NOTA: questa è solo una stima. I tempi di esecuzione effettivi potrebbero variare.)

Risultati di apprendimento (C++)

• Scopri come creare circuiti a profondità minore rispetto alla Trotterizzazione

• Percorri un workflow end-to-end per la stima dello stato fondamentale usando qDRIFT e SQD

• Scopri come usare qiskit-fermions insieme ad altri addon Qiskit per implementare tale workflow

Prerequisiti (C++)

• Leggi la panoramica Diagonalizzazione quantistica basata su campioni (SQD)

• Leggi la lezione Diagonalizzazione quantistica di Krylov basata su campioni (SKQD)

Contesto (C++)

SqDRIFT è una variante di SKQD che sostituisce la necessità di scegliere un ansatz da cui campionare bitstring con un insieme di circuiti di evoluzione temporale costruiti direttamente dalla Hamiltoniana target. Questo si ottiene sottocampionando operatori di evoluzione temporale più piccoli dalla Hamiltoniana in base ai suoi coefficienti, il che è noto come metodo di Trotterizzazione qDRIFT.

Entrambe le implementazioni usano Qiskit Fermions per rappresentare la Hamiltoniana molecolare e costruire i circuiti per l'algoritmo qDRIFT. L'implementazione Python usa pass del transpiler fermionici, mentre l'implementazione C++ campiona esplicitamente i gruppi e li mappa su operatori di qubit.

Sia la Hamiltoniana della forma:

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

dove, senza perdita di generalità, richiediamo $c_i > 0$ e che l'autovalore più grande di $h_i$ sia uguale, in valore assoluto, a $1$. Qualsiasi prefattore con segno o complesso viene assorbito in $h_i$, quindi i coefficienti $c_i$ sono pesi strettamente positivi mentre gli $h_i$ portano la direzione di ogni termine. Qui $N$ è il numero di termini (o, dopo il raggruppamento, il numero di gruppi) nella Hamiltoniana; è una proprietà della Hamiltoniana ed è distinto dal numero di operatori campionati in un singolo circuito, indicato con $n$ di seguito.

L'algoritmo qDRIFT realizza quindi, per il tempo target $t$, un operatore $V_k$, dove $k$ va da $1 \cdots K$ e indica il $k$-esimo circuito SqDRIFT, definito come:

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

Qui $n$ è il numero di operatori campionati per circuito e $K$ è il numero di circuiti nell'insieme. Il prodotto si estende sulle $n$ estrazioni, non su tutti gli $N$ termini della Hamiltoniana, e poiché i termini sono estratti con reinserimento, lo stesso $h_i$ può comparire più di una volta in un singolo $V_k$.

La quantità:

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

è la norma $L_1$ dei coefficienti, quindi ciascuno degli $n$ passi evolve per la stessa durata $\lambda t / n$ indipendentemente da quale termine è stato estratto. L'uniformità dell'angolo del passo è la caratteristica distintiva di qDRIFT: un coefficiente influenza il risultato tramite quanto spesso il suo termine viene estratto, non tramite quanto quel termine viene ruotato. Gli indici sono campionati dalla distribuzione:

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

quindi la serie $(k_1, \ldots, k_n)$ è una sequenza casuale di indici di termini estratti da questa distribuzione. Poiché i $c_i$ sono positivi e sommano a $\lambda$, questa è una distribuzione di probabilità normalizzata, e il valore atteso del canale risultante sulle estrazioni casuali approssima l'evoluzione sotto $H$, con un errore che diminuisce al crescere di $n$. Nota che l'errore di approssimazione dipende da $\lambda$ piuttosto che dal numero di termini $N$.

(L'articolo SqDRIFT scrive il numero di termini come $\mathcal{N}$ e la lunghezza della sequenza come $N$; qui usiamo $N$ e $n$ per mantenere i due chiaramente distinti.)

Questo tutorial mostra come generare un insieme di tali circuiti randomizzati. Dopo aver creato questi circuiti, in modo simile a come si crea un sottospazio di Krylov per operatori diversi, campioniamo bitstring da più di tali operatori con parametri temporali diversi. Questo assicura una maggiore sovrapposizione tra i vettori dello stato fondamentale e le bitstring campionate.

Requisiti (C++)

Questa implementazione viene eseguita su hardware IBM Quantum® e usa un eseguibile SBD esterno per la diagonalizzazione classica. Non include l'esempio con il simulatore Python.

Segui le istruzioni di build C++ per installare il compilatore, CMake, Rust, Boost e le dipendenze di Qiskit. Installa Rust con rustup: le dipendenze di Qiskit fissano la propria versione di Rust, che rustup scarica automaticamente, ma un Rust installato da Homebrew o da un gestore di pacchetti Linux potrebbe essere troppo vecchio. Il passaggio finale con SBD richiede inoltre MPI, OpenMP e BLAS/LAPACK. Windows al momento non è supportato da questa implementazione.

Usa il sorgente C++ completo e il progetto CMake dal repository. Mantieni la directory adiacente fcidump_files quando scarichi l'esempio. Da un checkout del repository della documentazione, compila con i comandi seguenti. La prima build scarica e compila tutte le dipendenze e richiede diversi minuti.

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

La build scrive l'eseguibile SqDRIFT in docs/tutorials/assets/sqdrift/cpp (la radice del progetto), non in build/. Esegui ./SqDRIFT dalla radice del progetto in modo che possa trovare ../fcidump_files/N2_sto_3g. Tutti gli altri comandi in questa scheda vengono eseguiti anch'essi dalla radice del progetto, salvo diversa indicazione. Se ./SqDRIFT fallisce con un errore di libreria condivisa mancante, segui le istruzioni specifiche della piattaforma sul percorso delle librerie nella guida alla build; su macOS, l'eseguibile di solito funziona senza di esse.

Configurazione (C++)

Gli snippet seguenti sono estratti del programma SqDRIFT.cpp completo. Compila ed esegui quel file come un unico programma.

Prima di eseguire SqDRIFT, assicurati che il tuo token API di IBM Quantum sia configurato per IBM Cloud Platform.

Se hai Python con qiskit-ibm-runtime installato, il modo più semplice per creare un file di configurazione formattato correttamente è:

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

Questo scrive ~/.qiskit/qiskit-ibm.json (Linux/macOS) o %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) con la struttura corretta. Se preferisci creare il file manualmente:

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

Ottieni la tua chiave API da IBM Quantum Platform.

Nota: Un file di configurazione mancante o malformato causa all'avvio un panic Rust unwrap() poco utile anziché un errore chiaro. Ricontrolla il nome della chiave (default-ibm-quantum-platform) e il valore del canale (ibm_quantum_platform) se incontri un problema.

Implementazione (C++)

Passo 1: Mappa gli input classici a un problema quantistico (C++)

Carica la Hamiltoniana molecolare

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)

Carica la Hamiltoniana della molecola $N_2$ da un file FCIDump. Il numero di qubit è pari al doppio del numero di orbitali (spin-up e spin-down).

Crea la Hamiltoniana e ordinala in forma normale

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

Converte i dati FCIDump in un operatore fermionico e applica l'ordinamento normale per semplificare la struttura dell'operatore.

Raggruppa i termini per struttura elettronica

// 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

Raggruppa i termini correlati della Hamiltoniana per struttura elettronica prima di costruire i loro circuiti di evoluzione.

Calcola i pesi di campionamento, poi normalizza e mappa sugli operatori di qubit

L'ordine conta qui. La normalizzazione riscala ogni coefficiente a magnitudine unitaria sul posto, quindi i pesi di campionamento devono essere calcolati per primi. Se la normalizzazione viene eseguita per prima, ogni $|c_j|$ risulta esattamente 1.0 e il peso di ciascun gruppo si riduce al semplice conteggio dei suoi termini — la distribuzione di campionamento sarebbe quindi determinata dalla dimensione del gruppo anziché dai coefficienti fisici.

Il peso di ciascun gruppo è il coefficiente assoluto medio dei suoi termini. Il pass QDriftTrotterization della scheda Python usa lo stesso peso per una Hamiltoniana raggruppata. Usare invece la somma farebbe crescere $\lambda$ con la dimensione del gruppo, il che a sua volta distorce il tempo di evoluzione $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

Il peso $\lambda_i = \frac{1}{|G_i|}\sum_{j \in G_i} |c_{ij}|$ è il coefficiente assoluto medio sui termini del gruppo $G_i$, letto dai coefficienti originali. La normalizzazione riscala poi ogni termine fermionico a magnitudine unitaria prima della trasformazione di Jordan-Wigner, preservando nell'operatore evoluto solo la sua fase/segno. Questo corrisponde all'evoluzione a termini normalizzati di qDRIFT: le magnitudini dei coefficienti determinano le probabilità di campionamento tramite $\lambda$, mentre l'evoluzione del circuito usa il termine normalizzato, così che i coefficienti grandi non vengano contati due volte.

Campiona gli insiemi di operatori

// 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...

Estrae gli insiemi di operatori per 100 circuiti casuali, ciascuno contenente 10 gruppi campionati dalla distribuzione pesata. Il passaggio successivo trasforma questi insiemi nell'ensemble di circuiti qDRIFT.

Passo 2: Costruisci i circuiti qDRIFT (C++)

Nella versione C++, questo passaggio costruisce i circuiti a partire dagli insiemi di operatori campionati. La transpilazione per il backend di destinazione avviene nel Passo 3, immediatamente prima dell'invio di ciascun circuito.

    // 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

Il programma completo stampa poi le statistiche dei circuiti e gli indici degli operatori del circuito 0, che qui sono omessi.

Per ogni operatore campionato, crea un circuito di evoluzione Suzuki-Trotter con tempo $t_k = \frac{\lambda \tau}{n}$, dove $n$ è il numero di operatori campionati per circuito e $\tau$ è il passo temporale. Poiché ogni gruppo campionato è stato normalizzato termine per termine prima della mappatura, questa evoluzione applica durante la rotazione solo la fase/segno di ciascun termine; le magnitudini contribuiscono soltanto attraverso la distribuzione di campionamento. Lo stato iniziale di Hartree-Fock viene preparato applicando gate X ai qubit 0-6 e 10-16 (sette elettroni in ciascun settore di spin). Le istruzioni di ogni circuito di evoluzione vengono aggiunte manualmente al circuito principale tramite l'API C++.

Passo 3: Eseguire usando le primitive di 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;

Si connette a IBM Quantum, seleziona il backend per nome, transpila ogni circuito per quel backend e invia un job per circuito: 100 job Sampler separati da 100 shot ciascuno su un registro di 20 qubit. Poi interroga ciascun job finché non termina e raccoglie le bitstring da result[0]. Ogni job usa tempo QPU e attende separatamente in coda, quindi aspettati che questo passaggio richieda molto più tempo di un singolo job in batch.

I risultati hardware riportati in questa pagina provengono da un'esecuzione su ibm_pittsburgh (100 circuiti × 100 shot), con backend_name modificato rispetto al valore predefinito ibm_fez. I campioni salvati sono stati post-elaborati con l'ordinamento dei bit dei determinanti corretto mostrato nel Passo 4. Il rumore hardware varia da un'esecuzione all'altra, quindi i tuoi conteggi e l'energia finale saranno diversi.

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

Cambiare il backend: Modifica backend_name in SqDRIFT.cpp per scegliere un dispositivo diverso; service.backends() elenca ciò che il tuo account può raggiungere.

Passo 4: Post-elaborare e restituire il risultato nel formato classico desiderato (C++)

Postseleziona le bitstring e scrivi le stringhe 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!

Postseleziona le bitstring con il corretto peso di Hamming (n_alpha elettroni spin-up, n_beta elettroni spin-down, derivati dal campo NELEC di FCIDump), le converte in stringhe di Configuration Interaction (CI) usando la simmetrizzazione di spin e le scrive in alphadets_from_sqd.txt per la successiva Selected Basis Diagonalization.

Libera le risorse

// 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;

Libera correttamente tutta la memoria e le risorse allocate. Il std::vector<QuantumCircuit> C++ viene ripulito automaticamente tramite RAII. Libera manualmente gli operatori di qubit, gli operatori fermionici e i dati FCIDump.

Ora che abbiamo creato la base su cui proiettare la nostra Hamiltoniana, possiamo procedere con il processo di diagonalizzazione per ottenere la stima dello stato fondamentale. Esegui i comandi seguenti dalla radice del progetto, docs/tutorials/assets/sqdrift/cpp.

Compila il binario SBD diag

CMake scarica SBD solo come sorgente in deps/sbd/ (la build principale salta intenzionalmente la sua compilazione). L'app SBD usa un semplice Makefile, e SBD stesso richiede MPI, OpenMP e BLAS/LAPACK — quindi un semplice make non è sufficiente finché questi non sono installati e il file Configuration dell'app non corrisponde alla tua toolchain.

Il Makefile legge i flag del compilatore e del linker dal file Configuration nella stessa directory. La versione distribuita upstream è calibrata su uno specifico layout macOS/Homebrew e non verrà compilata così com'è su una macchina standard — in particolare -I/opt/homebrew/opt/llvm/include non contiene omp.h (si trova in una directory di risorse clang specifica per versione), e clang++ di Apple rifiuta -fopenmp del tutto:

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

La procedura seguente è verificata funzionante su macOS (Apple Silicon, Darwin 25.6, Open MPI 5.x + Homebrew LLVM).

macOS (Apple Silicon) — testato

# 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

Poi sostituisci deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization/Configuration con:

# 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 ../../../..

Una build riuscita stampa le due righe di compilazione/link e lascia un binario diag accanto al Makefile. Un ld: warning: ignoring duplicate libraries: '-lomp' è innocuo.

macOS Intel: Il prefisso di Homebrew è /usr/local anziché /opt/homebrew, quindi sostituiscilo in entrambi i percorsi precedenti.

Linux (Ubuntu/Debian) — non testato

I prerequisiti seguenti coprono MPI, OpenMP e BLAS/LAPACK. GCC supporta -fopenmp in modo nativo, quindi non serve alcuna sostituzione del compilatore. Questo non è stato verificato su 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

Esegui SqDRIFT per generare le stringhe CI

Se hai già eseguito ./SqDRIFT attraverso i Passi 1–4, alphadets_from_sqd.txt esiste già nella radice del progetto e puoi saltare questo passaggio. Altrimenti, eseguilo dalla radice del progetto:

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

Esegui la diagonalizzazione

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 ../../../..

Parametri:

• --fcidump: percorso del file FCIDump contenente la Hamiltoniana
• --adetfile: percorso del file dei determinanti alpha (stringhe CI da SqDRIFT)
• --method 0: metodo di diagonalizzazione di Davidson
• --iteration 100: numero massimo di iterazioni di Davidson
• --block 10: dimensione del blocco per l'algoritmo di Davidson
• --tolerance 1e-8: tolleranza di convergenza per l'energia

 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 riga Energy = -107.5108722590187 è la stima dello stato fondamentale sul sottospazio campionato.

Come differisce dall'esempio Python

I due esempi differiscono sotto diversi aspetti contemporaneamente, non solo nel modo in cui trattano le bitstring con simmetria errata, quindi non ci si aspetta che producano energie identiche:

	Scheda C++	Scheda Python
Bitstring con simmetria errata	Scartate da postselect_bitstrings in base al peso di Hamming	Riparate dal recupero iterativo della configurazione
Struttura del ciclo	Singolo passaggio: campiona → postseleziona → diagonalizza una volta	Ciclo esterno: diagonalizza, leggi le occupazioni orbitali, recupera le configurazioni, ridiagonalizza
Numero di circuiti	100 circuiti (num_circuits = 100)	400 circuiti (200 per tempo di evoluzione)
Tempi di evoluzione	Un singolo $\tau$ (time_step = 1)	Due tempi, times = [1.0, 10.0]
Termini diagonali	Mantenuti nell'operatore campionato	Rimossi da filter_diagonal_terms, cambiando la distribuzione di campionamento e i circuiti
Ordinamento dei termini	Solo raggruppati	canonical_order applicato prima del raggruppamento per la riproducibilità
Solver di diagonalizzazione	SBD (diag, binario MPI/OpenMP esterno)	Solver in-process qiskit-addon-sqd
Backend	Fissato per nome (backend_name = "ibm_fez"; l'esecuzione in questa pagina ha usato ibm_pittsburgh)	Backend operativo least_busy con un numero sufficiente di qubit
Transpilazione	Ogni circuito separatamente, livello di ottimizzazione 2	Tutti i circuiti insieme, livello di ottimizzazione 3
Invio dei job	100 job, un circuito ciascuno	Un job Sampler contenente tutti i 400 circuiti
Esecuzione su simulatore	Nessuna; solo hardware	Esempio con simulatore senza rumore prima dell'esecuzione hardware

Entrambe le versioni usano la stessa Hamiltoniana, ne raggruppano i termini per struttura elettronica, pesano ogni gruppo con il coefficiente assoluto medio dei suoi termini, partono dallo stato di Hartree-Fock e campionano 10 gruppi e 100 shot per circuito.

Ciascuna differenza influisce indipendentemente sul sottospazio campionato, e quindi sull'energia. L'ensemble più grande e il secondo tempo di evoluzione, più lungo, ampliano entrambi le configurazioni esplorate dall'esempio Python; campionare a più di un tempo conferisce a SqDRIFT la sua copertura di tipo Krylov dello spazio a bassa energia. Il filtraggio dei termini diagonali cambia quali operatori il campionatore può estrarre e i loro pesi, quindi le due esecuzioni campionano distribuzioni diverse.

La postselezione è un filtro rigido, e uno shot con il numero di particelle sbagliato viene scartato del tutto. Nell'esecuzione sopra ciò ha scartato la maggior parte dei dati: 1137 shot su 10000 hanno superato il controllo del peso di Hamming (7,7) (circa l'11%), condensandosi in soli 49 $\alpha$-determinanti distinti. Il recupero della configurazione, invece, usa le occupazioni orbitali medie di una diagonalizzazione precedente per invertire i bit e riparare tali shot, riciclando ciò che la postselezione scarterebbe. Questo feedback lo rende intrinsecamente iterativo: recover_configurations in qiskit-addon-sqd-hpc richiede un argomento avg_occupancies che esiste solo dopo l'esecuzione di una diagonalizzazione.

Questa esecuzione C++ produce un sottospazio di 49 $\alpha$-determinanti e un'energia di -107.5109 Ha. L'intero spazio degli $\alpha$-determinanti per sette elettroni in 10 orbitali è di soli $\binom{10}{7} = 120$ stringhe, quindi fornendo tutte le 120 allo stesso binario diag si ottiene l'energia FCI esatta -107.6482 Ha. L'esecuzione si colloca quindi circa 0.1373 Ha sopra il valore esatto, pur coprendo 49 dei 120 $\alpha$-determinanti. Poiché SQD è variazionale nel sottospazio selezionato, un sottospazio più piccolo può solo alzare l'energia, quindi il sottospazio incompleto è la lettura più diretta del divario.

Questo controllo FCI esatto è possibile solo perché l'esempio è piccolo; non fa parte del flusso di lavoro SqDRIFT, e per le dimensioni di sistema in cui SQD è effettivamente necessario non esiste alcun riferimento di questo tipo. Qui serve solo a quantificare quanta parte dello spazio ha recuperato questa esecuzione.

Un confronto con l'esempio Python non è controllato, quindi una differenza non dovrebbe essere attribuita alla sola assenza del recupero della configurazione; nota che il numero di circuiti, il singolo tempo di evoluzione, il filtraggio dei termini diagonali e il rumore hardware differiscono tutti contemporaneamente e tutti influenzano il sottospazio.

Per colmare il divario, ecco alcuni possibili passi successivi da provare dopo aver completato questo tutorial: aumentare shots, num_circuits o ops_per_circuit; campionare a più di un time_step; e aggiungere un ciclo di recupero sopra Qiskit::addon::sqd::recover_configurations.

Prossimi passi

Raccomandazioni

• Diagonalizzazione quantistica di Krylov basata su campioni di un modello reticolare fermionico — un tutorial correlato che usa circuiti di evoluzione temporale anziché un ansatz variazionale.

• Diagonalizzazione quantistica basata su campioni di una Hamiltoniana chimica — un tutorial su come costruire un circuito local unitary cluster Jastrow (LUCJ) per la simulazione di chimica quantistica.

• Il paper SqDRIFT — la letteratura su cui si basa questo tutorial. (Nota che alcune delle ottimizzazioni discusse in questo paper sono attualmente in corso di sviluppo, e questo tutorial è soggetto a modifiche in futuro in base all'evoluzione delle librerie utilizzate.)