SqDRIFT algoritmo per la stima dello stato fondamentale

Stima del tempo di esecuzione: 180 secondi su un processore Heron r3 (NOTA: si tratta solo di una stima. (La durata potrebbe variare.)

Risultati di apprendimento

• Scopri come realizzare circuiti con profondità inferiore rispetto alla tecnica di Trotter
• Esplora un flusso di lavoro completo per la stima dello stato fondamentale utilizzando l' qDRIFT e e l'SQD
• Scopri come utilizzarlo in qiskit-fermions combinazione con altri componenti aggiuntivi di Qiskit per implementare un flusso di lavoro di questo tipo

Prerequisiti

• Leggi la panoramica sulla diagonalizzazione quantistica basata su campioni (SQD)
• Leggi la lezione sulla diagonalizzazione quantistica di Krylov basata su campioni (SKQD)

Sfondo

SqDRIFT è una variante dello SKQD che sostituisce la necessità di scegliere un ansatz da cui campionare le stringhe di bit con un insieme di circuiti di evoluzione temporale costruiti direttamente a partire dall'Hamiltoniano di destinazione. Ciò si ottiene mediante il sottocampionamento di operatori di evoluzione temporale più piccoli a partire dall’Hamiltoniano, sulla base dei suoi coefficienti; tale metodo è noto come metodo di Trotterizzazione “ qDRIFT ”.

Entrambe le implementazioni utilizzano Qiskit Fermions per rappresentare l'Hamiltoniano molecolare e costruire circuiti per l'algoritmo " qDRIFT ". L'implementazione in Python utilizza passaggi di transpiler fermionici, mentre l'implementazione in C++ campiona esplicitamente i gruppi e li mappa agli operatori sui qubit.

Sia l'Hamiltoniano della forma:

dove, senza perdita di generalità, si pone $c_i > 0$ e si richiede che il massimo autovalore di $h_i$ sia uguale, in valore assoluto, a $1$. Qualsiasi prefattore con segno o complesso viene assorbito in $h_i$, per cui i coefficienti $c_i$ sono pesi strettamente positivi mentre i $h_i$ indicano la direzione di ciascun termine. Qui $N$ indica il numero di termini (o, dopo il raggruppamento, il numero di gruppi) nell’Hamiltoniano; si tratta di una proprietà dell’Hamiltoniano ed è distinta dal numero di operatori campionati in un singolo circuito, indicato di seguito c $n$.

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

Qui $n$ indica il numero di operatori campionati per circuito e $K$ indica il numero di circuiti nell'insieme. Il prodotto viene calcolato sulle estrazioni $n$, non su tutti i termini hamiltoniani $N$; inoltre, poiché le estrazioni avvengono con sostituzione, lo stesso $h_i$ può comparire più di una volta in un singolo $V_k$.

La quantità:

è la norma $L_1$ dei coefficienti, quindi ciascuno dei passi $n$ si evolve per la stessa durata $\lambda t / n$ indipendentemente dal termine estratto. L'uniformità dell'angolo di passo è la caratteristica distintiva dell' qDRIFT: : un coefficiente influenza il risultato in base alla frequenza con cui il suo termine viene estratto, non in base all'ampiezza della rotazione di tale termine. Gli indici sono estratti dalla distribuzione:

quindi la serie $(k_1, \ldots, k_n)$ è una sequenza casuale di indici di termini estratti da questa distribuzione. Poiché i valori di $c_i$ sono positivi e la loro somma è pari a $\lambda$, si tratta di una distribuzione di probabilità normalizzata, e l’aspettativa del canale risultante rispetto alle estrazioni casuali approssima l’evoluzione in base a $H$, con un errore che diminuisce all’aumentare di $n$. Si noti che l'errore di approssimazione dipende da $\lambda$ piuttosto che dal numero di termini $N$.

(L'articolo pubblicato su * SqDRIFT * indica il numero di termini come $\mathcal{N}$ e la lunghezza della sequenza come $N$; qui utilizziamo $N$ e $n$ per distinguere chiaramente le due cose.)

Questo tutorial illustra come generare un insieme di tali circuiti randomizzati. Dopo aver creato questi circuiti, analogamente a come si crea un sottospazio di Krylov per diversi operatori, campioniamo stringhe di bit provenienti da più operatori di questo tipo con parametri temporali diversi. Ciò garantisce una maggiore sovrapposizione tra i vettori dello stato fondamentale e le stringhe di bit 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 (Si noti che il nome è al plurale)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

È possibile installare tutti i pacchetti necessari con:

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

Configura

# 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

Fase 1: Mappare gli input classici su un problema quantistico

I passaggi riportati di seguito relativi a “ Python ” utilizzano il simulatore Aer. Di seguito è riportato l'esempio relativo all'hardware.

Lettura e preparazione dell'FCIDump

Per questo tutorial, caricheremo l'Hamiltoniano della struttura elettronica dell'azoto ( N2 ). Esistono anche altri modi per creare operatori fermionici. Consultare la documentazione all'indirizzo qiskit_fermions.operators.library.

Informazioni su questo FCIDump. Il file N2_sto_3g descrive una molecola di azoto ( $N_2$ ) nella base minima STO-3G con una distanza interatomica di 1.09 $\AA$, ovvero la lunghezza di legame di equilibrio sperimentale. L'intestazione dichiara NORB=10, NELEC=14, e MS2=0: 10 orbitali spaziali (quindi 20 orbitali di spin e 20 qubit secondo la rappresentazione di Jordan-Wigner), 14 elettroni in uno stato di singoletto di spin, ovvero sette elettroni $\alpha$ e sette $\beta$. A tutti gli orbitali viene assegnata l'etichetta di simmetria 1, ovvero non viene sfruttata alcuna simmetria di gruppo puntuale. Trattandosi di un dump “ STO-3G ” a spazio pieno, nessun orbitale è congelato e lo spazio di correlazione è sufficientemente piccolo da consentire il calcolo classico di un’energia di riferimento FCI esatta a scopo di confronto, come illustrato nel blocco di codice seguente.

È possibile rigenerare un file equivalente con il comando 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 convergenti, un file rigenerato potrebbe differire da quello fornito in termini di fase o ordine degli orbitali; le energie totali non ne risentono.

Come procurarsi il file. Trova FCIDump in questo repository GitHub. Puoi eseguire il codice riportato di seguito per scaricarlo nella posizione indicata nel resto del tutorial.

Per prima cosa utilizziamo la funzione fornita cisolver da pyscf per ottenere l'energia di riferimento. Questa è la vera energia dello stato fondamentale della molecola su cui stiamo lavorando. A tal fine, dichiareremo innanzitutto norb e nelec, che rappresentano rispettivamente il numero di orbitali e il numero di elettroni. Quindi definiamo e h1e h2e, che sono rispettivamente gli integrali a un elettrone e a due elettroni. Tutti questi elementi verranno utilizzati in seguito anche per l'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'Hamiltoniano

Una volta preparati i dati necessari, leggiamo l’Hamiltoniano 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

Flussi di lavoro fermionici con qiskit-fermions

Per prima cosa mapperemo l'Hamiltoniano su un modello di circuito fermionico utilizzando qiskit-fermions, che fornisce passaggi di transpiler e porte specifiche per i circuiti fermionici. Questi verranno successivamente utilizzati prima che il transpiler tradizionale di Qiskit esegua le fasi previste da questo flusso di lavoro.

Raggruppamento per termini

Per garantire la riproducibilità dei risultati, utilizziamo innanzitutto canonical_order per ordinare i termini basandoci esclusivamente sulla loro struttura. L'ordine degli operatori nell'elenco canon è quindi fisso. Ciò garantisce la riproducibilità degli operatori creati, poiché il pass QDriftTrotterization che utilizzeremo negli esempi successivi estrae indici casuali per generare gli operatori qDRIFT.

In questa fase, sfruttiamo le numerose simmetrie presenti nell'Hamiltoniano della struttura elettronica raggruppando i termini correlati con coefficienti identici. Sebbene ciò modifichi la distribuzione dei coefficienti degli operatori da cui il protocollo “ qDRIFT ” preleva i campioni, ciò non influisce sulle sue garanzie di convergenza. È fondamentale sottolineare che il raggruppamento dei termini correlati per simmetria comporta una cancellazione vantaggiosa 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 al posto nostro.

Si noti che la formula group_terms_by_electronic_structure presuppone che i termini siano in ordine normale.

Filtraggio dei termini diagonali

Eliminiamo i termini diagonali dall’Hamiltoniano utilizzato per generare i circuiti, in modo che $n$ qDRIFT gli slot di campionamento siano dedicati ai termini che spostano la popolazione tra le configurazioni. È preferibile eliminare tali termini dall’Hamiltoniano in questa fase, prima che il gate Evolution venga costruito nella fase successiva.

I termini in questione sono quelli diagonali nella base "occupazione-numero", ovvero i prodotti degli operatori numerici $a^\dagger_i a_i$. Rientrano in questa descrizione tre tipi di termini:

• l' offset energetico costante, un prodotto di operatori a numero zero, la cui evoluzione temporale contribuisce solo con una fase globale;
• gli operatori a numero singolo $n_i$, la cui evoluzione temporale si riduce a rotazioni di un singolo qubit $Z$;
• i prodotti di ordine superiore, come ad esempio $n_i n_j$.

Di per sé, nessuna di queste mosse fa passare la popolazione da una configurazione occupazionale all’altra; esse agiscono solo sulle fasi delle configurazioni già presenti. Tuttavia, non sono inerti: quelle fasi relative contribuiscono all’interferenza generata dai termini di eccitazione più avanti nel circuito, quindi filtrarle modifica l’evoluzione effettivamente generata e può alterare la distribuzione di campionamento. Si tratta di un'approssimazione intenzionale nella fase di generazione del circuito, volta a concentrare il campionamento sui termini di eccitazione, piuttosto che di una fase che lasci inalterata la distribuzione campionata. A differenza del raggruppamento di simmetria sopra descritto, che lascia intatte le garanzie di convergenza di tipo “ qDRIFT ”, questo filtro modifica l’operatore oggetto dell’evoluzione. I circuiti, quindi, non approssimano più l’evoluzione nell’ambito dell’Hamiltoniano completo, e i limiti di errore di tipo “ qDRIFT ” si applicano all’operatore filtrato anziché a quello originale. In questo caso ciò è accettabile perché i circuiti rappresentano solo un metodo euristico di campionamento utilizzato per proporre configurazioni: nessun termine viene perso nella stima dell'energia stessa, poiché il filtro si applica solo all'Hamiltoniano utilizzato per costruire i circuiti, mentre la diagonalizzazione classica successiva utilizza l'Hamiltoniano completo, termini diagonali inclusi. L'accuratezza dell'SQD dipende da quel passo classico, che rimane variazionale nel sottospazio campionato indipendentemente dal modo in cui sono state proposte le configurazioni.

La funzione filter_diagonal_terms() rimuove tali termini da un operatore direttamente sul posto. Li identifica in base alla loro struttura ordinata normalmente — il multinseme delle modalità di creazione che corrisponde al multinseme delle modalità di annichilimento — quindi è valida solo su un operatore che sia già ordinato normalmente. Questa ipotesi non viene verificata durante l'esecuzione.

# 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'Hamiltoniano, stabiliremo i seguenti parametri per generare l'insieme dei circuiti:

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

Creazione di circuiti fermionici

Ora creeremo circuiti fermionici per ciascuno degli intervalli di tempo. Ogni circuito sarà costituito da un unico gate di evoluzione, con il tempo di evoluzione che abbiamo indicato in precedenza. L'operatore di evoluzione è l'hamiltoniano. Successivamente eseguiamo i passaggi del transpiler su questi circuiti per creare circuiti di tip qDRIFT.

Preparazione dell'Ansatz

Prepariamo lo stato di Hartree-Fock utilizzando la classe InitializeModes . Per l'azoto, il processo consiste semplicemente nell'applicare X gate ai primi num_elec_a qubit e poi ai qubit num_elec_b , entrambi pari a sette nel caso dell'azoto. Questo stato rappresenta i sette elettroni $\alpha$ e i sette elettroni $\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)

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

Ora che disponiamo dei nostri circuiti, utilizzeremo innanzitutto i passaggi disponibili in qiskit-fermions per eseguire ottimizzazioni a livello fermionico, per poi traspilare il nostro circuito per il backend di nostra scelta. Trattandosi di un esperimento al simulatore, lo faremo innanzitutto per l’ AerSimulator.

Calcolo del peso per ciascun gruppo

In questa fase, eseguiamo il campionamento “ qDRIFT ” dei termini in modo stocastico, con probabilità proporzionali ai loro coefficienti nell’Hamiltoniano. Il passaggio del transpiler “ qDRIFT ” lo fa per noi. Ora siamo in grado di creare circuiti meno profondi che possono essere eseguiti sull'hardware in modo più efficiente nonostante la connettività limitata dei qubit, anche quando l'hamiltoniano contiene accoppiamenti a lungo raggio e termini superiori al quadrato. Dopo il raggruppamento per termini, effettua un campionamento degli operatori in base ai loro pesi. Per ciascun operatore $h_i$, il peso $W_{h_i}$ è definito come segue:

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

Ottimizzazioni fermioniche e native per l'hardware

La funzione restituisce generate_preset_jw_pass_manager() un oggetto MultiStagePassManager che accetta un parametro e FermionicCircuit produce un circuito finale ottimizzato che possiamo transpilare per eseguirlo sul nostro hardware. Sostituiamo la sua fase di ottimizzazione predefinita con una che contenga FermionicPassManager il nostro passaggio QDriftTrotterization :

• Il pass QDriftTrotterization utilizza internamente il calcolo del peso e il campionamento per generare i circuiti che useremo per il campionamento
• Il passaggio RelabelModes è un’ulteriore fase di ottimizzazione che può essere utilizzata per permutare le modalità fermioniche al fine di ottimizzare la connettività tra i qubit e ridurre la profondità dei gate; per ulteriori informazioni, consultare la documentazione di riferimento dell’API

Le fasi restanti vengono eseguite MultiStagePassManager automaticamente e gestiscono l'intera mappatura da fermioni a qubit:

• F2QLayout : Il gestore di passaggi predefinito applica il passaggio TrivialF2QLayout , che mappa in modo banale i bit fermionici $n$ i ai qubit $n$ i.
• F2QSynth : un passaggio di transpilazione per mappare le istruzioni dei circuiti basati sui fermioni a quelle basate sui 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 completato le ottimizzazioni a livello fermionico, possiamo traspilare i circuiti per eseguirli sul simulatore.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Fase 3: Eseguire il comando utilizzando Qiskit primitives

Ora che abbiamo i nostri circuiti, possiamo eseguirli utilizzando Qiskit primitives su AerSimulator. Sommaremo tutti i conteggi provenienti dai diversi circuiti. Li convertiamo in vettori booleani prima di procedere alla post-elaborazione finale 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

Fase 4: Elaborazione finale e restituzione del risultato nel formato classico desiderato

Utilizzo delle stringhe di bit per SQD

Ora possiamo applicare lo schema di diagonalizzazione alle stringhe di bit selezionate per individuare 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 utilizzata per visualizzare l'iterazione corrente e la stima dell'autovalore corrente ad ogni iterazione.

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

Nota : la dimensione del sottospazio non è fissa da un'iterazione all'altra, nemmeno sul simulatore privo di rumore: ogni sottocampione genera un insieme diverso di configurazioni e la fase di recupero rimodella il pool tra un'iterazione e l'altra, pertanto la dimensione riportata varia da un sottocampione all'altro. Il campionamento silenzioso non determina di per sé la dimensione del sottospazio selezionato. L'esecuzione su hardware, tuttavia, tende a produrre sottospazi sistematicamente più grandi, poiché gli scatti rumorosi rompono la simmetria del numero di particelle e il recupero della configurazione li trasforma in vettori di base aggiuntivi. Per questo motivo, nella sezione dedicata all'hardware introdurremo anche un'ulteriore fase per la potatura delle stringhe di bit.

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 di hardware

In questo esempio vengono utilizzati 20 qubit (10 orbitali spaziali). Questa scelta è dettata da ragioni di praticità per un tutorial che dovrebbe essere eseguito rapidamente, non rappresenta un limite inderogabile per il metodo.

Il costo della fase classica non dipende direttamente dal numero di qubit. SQD diagonalizza l’Hamiltoniano proiettato sul sottospazio generato dalle configurazioni campionate, quindi ciò che determina il costo classico è la dimensione di quel sottospazio selezionato — determinata in questo caso da samples_per_batch, num_batches, e dal numero di configurazioni distinte effettivamente prodotte dai circuiti — insieme all’algebra lineare sparsa necessaria per applicare l’Hamiltoniano proiettato. Lo spazio CI completo cresce in modo combinatorio con il numero di orbitali e di elettroni, ma il sottospazio selezionato ne costituisce una piccola porzione regolabile, di cui controlliamo direttamente le dimensioni. Di conseguenza, il numero di qubit e la difficoltà classica possono essere variati in modo relativamente indipendente: uno spazio orbitale più ampio, campionato in un sottospazio di dimensioni modeste, può risultare meno oneroso rispetto a un sistema più piccolo diagonalizzato su uno spazio molto ampio.

In pratica, quindi, la dimensione massima del sistema dipende dalla dimensione del sottospazio necessaria per ottenere la precisione desiderata e dalla memoria e dai core a disposizione del risolutore di autovalori. Spazi orbitali più ampi richiedono in genere un sottospazio più ampio per raggiungere l'accuratezza chimica, ed è proprio questo che alla fine giustifica l'uso di risorse distribuite — si veda qiskit-addon-sqd-hpc per la scalabilità di questa fase. Anziché ipotizzare un valore limite fisso, l'approccio pratico consiste nell'osservare la dimensione del sottospazio riportata e la convergenza dell'energia nel corso delle iterazioni, aumentando la dimensione del sottospazio fino a quando l'energia smette di migliorare o si esaurisce la memoria disponibile.

Nota: a causa dell'errore di campionamento dovuto al rumore presente nell'hardware, il sottospazio creato per la diagonalizzazione nell'esecuzione su hardware risulterà più ampio rispetto a quello ottenuto utilizzando il simulatore. Sebbene aumenti la dimensione del sottospazio che vogliamo diagonalizzare, il flusso di lavoro ci fornisce comunque una risposta accurata grazie alla robustezza dell'SQD nei confronti del rumore.

Eliminazione delle stringhe superflue

A questo punto possiamo scegliere di eseguire un passaggio aggiuntivo. Una volta ottenute tutte le stringhe di bit derivanti dall’esecuzione dei circuiti, possiamo scegliere se filtrare le stringhe di bit non valide prima di eseguire l’SQD oppure procedere senza effettuare alcuna eliminazione. In genere, è preferibile non effettuare la potatura nelle esecuzioni hardware, poiché in questo modo le soluzioni con simmetria violata rimangono disponibili per il recupero della configurazione, che può ripararle trasformandole in configurazioni valide e quindi ampliare il sottospazio, anziché scartarle del tutto.

Poiché l'azoto può avere solo sette $\alpha$ e sette $\beta$ elettroni, qualsiasi stringa di bit che contenga più o meno di sette 1s nella prima e nella seconda metà dell'output può essere scartata. Definiamo una funzione che verifica se le stringhe di bit sono valide e, in caso contrario, le scarta. Una volta eliminate le stringhe di bit spurie, quelle rimanenti vengono inviate allo schema di diagonalizzazione. Utilizza il flag PRUNE riportato di seguito per passare da un comportamento all'altro.

Tieni presente che la potatura è solo una delle diverse opzioni che determinano il sottospazio finale, insieme al numero di circuiti, all’insieme dei tempi di evoluzione e al filtraggio dei termini diagonali. Il confronto tra un'esecuzione potata e una non potata è significativo solo se tutti gli altri fattori rimangono costanti; la scheda dedicata al C++ approfondisce questo aspetto, poiché tale linguaggio effettua una post-selezione anziché un recupero e presenta differenze anche in relazione a 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

Stima del tempo di esecuzione: 180 secondi su un processore Heron r3 (NOTA: si tratta solo di una stima. (La durata potrebbe variare.)

Risultati di apprendimento (C++)

• Scopri come realizzare circuiti con profondità inferiore rispetto alla tecnica di Trotter
• Esplora un flusso di lavoro completo per la stima dello stato fondamentale utilizzando l' qDRIFT e e l'SQD
• Scopri come utilizzarlo in qiskit-fermions combinazione con altri componenti aggiuntivi di Qiskit per implementare un flusso di lavoro di questo tipo

Prerequisiti (C++)

• Leggi la panoramica sulla diagonalizzazione quantistica basata su campioni (SQD)
• Leggi la lezione sulla diagonalizzazione quantistica di Krylov basata su campioni (SKQD)

Contesto (C++)

SqDRIFT è una variante dello SKQD che sostituisce la necessità di scegliere un ansatz da cui campionare le stringhe di bit con un insieme di circuiti di evoluzione temporale costruiti direttamente a partire dall'Hamiltoniano di destinazione. Ciò si ottiene mediante il sottocampionamento di operatori di evoluzione temporale più piccoli a partire dall’Hamiltoniano, sulla base dei suoi coefficienti; tale metodo è noto come metodo di Trotterizzazione “ qDRIFT ”.

Entrambe le implementazioni utilizzano Qiskit Fermions per rappresentare l'Hamiltoniano molecolare e costruire circuiti per l'algoritmo " qDRIFT ". L'implementazione in Python utilizza passaggi di transpiler fermionici, mentre l'implementazione in C++ campiona esplicitamente i gruppi e li mappa agli operatori sui qubit.

Sia l'Hamiltoniano della forma:

dove, senza perdita di generalità, si pone $c_i > 0$ e si richiede che il massimo autovalore di $h_i$ sia uguale, in valore assoluto, a $1$. Qualsiasi prefattore con segno o complesso viene assorbito in $h_i$, per cui i coefficienti $c_i$ sono pesi strettamente positivi mentre i $h_i$ indicano la direzione di ciascun termine. Qui $N$ indica il numero di termini (o, dopo il raggruppamento, il numero di gruppi) nell’Hamiltoniano; si tratta di una proprietà dell’Hamiltoniano ed è distinta dal numero di operatori campionati in un singolo circuito, indicato di seguito c $n$.

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

Qui $n$ indica il numero di operatori campionati per circuito e $K$ indica il numero di circuiti nell'insieme. Il prodotto viene calcolato sulle estrazioni $n$, non su tutti i termini hamiltoniani $N$; inoltre, poiché le estrazioni avvengono con sostituzione, lo stesso $h_i$ può comparire più di una volta in un singolo $V_k$.

La quantità:

è la norma $L_1$ dei coefficienti, quindi ciascuno dei passi $n$ si evolve per la stessa durata $\lambda t / n$ indipendentemente dal termine estratto. L'uniformità dell'angolo di passo è la caratteristica distintiva dell' qDRIFT: : un coefficiente influenza il risultato in base alla frequenza con cui il suo termine viene estratto, non in base all'ampiezza della rotazione di tale termine. Gli indici sono estratti dalla distribuzione:

quindi la serie $(k_1, \ldots, k_n)$ è una sequenza casuale di indici di termini estratti da questa distribuzione. Poiché i valori di $c_i$ sono positivi e la loro somma è pari a $\lambda$, si tratta di una distribuzione di probabilità normalizzata, e l’aspettativa del canale risultante rispetto alle estrazioni casuali approssima l’evoluzione in base a $H$, con un errore che diminuisce all’aumentare di $n$. Si noti che l'errore di approssimazione dipende da $\lambda$ piuttosto che dal numero di termini $N$.

(L'articolo pubblicato su * SqDRIFT * indica il numero di termini come $\mathcal{N}$ e la lunghezza della sequenza come $N$; qui utilizziamo $N$ e $n$ per distinguere chiaramente le due cose.)

Questo tutorial illustra come generare un insieme di tali circuiti randomizzati. Dopo aver creato questi circuiti, analogamente a come si crea un sottospazio di Krylov per diversi operatori, campioniamo stringhe di bit provenienti da più operatori di questo tipo con parametri temporali diversi. Ciò garantisce una maggiore sovrapposizione tra i vettori dello stato fondamentale e le stringhe di bit campionate.

Requisiti (C++)

Questa implementazione gira su hardware IBM Quantum® e utilizza un eseguibile SBD esterno per la diagonalizzazione classica. Non include l'esempio del simulatore " Python ".

Segui le istruzioni di compilazione in C++ per installare il compilatore, CMake, Rust, Boost e le dipendenze di Qiskit. Installa Rust con rustup : le dipendenze di Qiskit specificano una versione propria di Rust, che rustup scarica automaticamente; tuttavia, una versione di Rust installata tramite Homebrew o un gestore di pacchetti Linux potrebbe essere troppo vecchia. La fase finale dell'SBD richiede inoltre MPI, l' OpenMP, e e BLAS/LAPACK. Windows non è attualmente supportato da questa implementazione.

Utilizza il codice sorgente completo in C++ e il progetto CMake presenti nel repository. Quando scarichi l'esempio, mantieni la directory fcidump_files adiacente. Dopo aver effettuato il checkout del repository della documentazione, eseguire la compilazione con i seguenti comandi. La prima compilazione scarica e compila tutte le dipendenze e richiede alcuni minuti.

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

La compilazione salva l'eseguibile SqDRIFT in docs/tutorials/assets/sqdrift/cpp (la directory principale del progetto), non in build/. Eseguire ./SqDRIFT dalla directory principale del progetto in modo che possa individuare ../fcidump_files/N2_sto_3g. Anche tutti gli altri comandi presenti in questa scheda vengono eseguiti dalla directory principale del progetto, salvo diversa indicazione. Se l'esecuzione fallisce ./SqDRIFT a causa di un errore relativo alla mancanza di una libreria condivisa, seguire le istruzioni relative al percorso delle librerie specifiche per la piattaforma riportate nella guida alla compilazione; su macOS, l'eseguibile di solito funziona senza di esse.

Configurazione (C++)

I brani riportati di seguito sono estratti dal programma completo di “ SqDRIFT.cpp ”. Compila ed esegui quel file come un unico programma.

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

Se hai installato qiskit-ibm-runtime Python, 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 ( Linux/macOS~/.qiskit/qiskit-ibm.json ) o %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) con la forma 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 non valido provoca un " unwrap() panic" di Rust poco chiaro all'avvio, anziché un errore chiaro. Se riscontri un problema, ricontrolla il nome della chiave (default-ibm-quantum-platform) e il valore del canale (ibm_quantum_platform).

Implementazione (C++)

Fase 1: Mappare gli input classici su un problema quantistico (C++)

Carica l'Hamiltoniano 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 l'hamiltoniano della molecola " $N_2$ " da un file FCIDump. Il numero di qubit è pari al doppio del numero di orbitali (spin-up e spin-down).

Creare e ordinare secondo l'ordine normale l'hamiltoniano

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

Classificazione dei gruppi in base alla 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

I gruppi hanno messo in relazione i termini hamiltoniani con la struttura elettronica prima di costruire i relativi circuiti di evoluzione.

Calcolare i pesi di campionamento, quindi normalizzarli e mapparli sugli operatori dei qubit

Qui l'ordine è importante. La normalizzazione ridimensiona ogni coefficiente portandolo a grandezza unitaria direttamente, pertanto è necessario calcolare prima i pesi di campionamento. Se la normalizzazione viene eseguita per prima, ogni $|c_j|$ viene interpretata esattamente come 1.0 e il peso di ciascun gruppo si riduce al semplice conteggio dei suoi termini — la distribuzione campionaria sarebbe quindi determinata dalla dimensione del gruppo piuttosto che dai coefficienti fisici.

Il peso di ciascun gruppo è dato dal coefficiente assoluto medio dei suoi termini. Il passaggio QDriftTrotterization nella scheda " Python " utilizza lo stesso peso per un hamiltoniano raggruppato. Se si utilizzasse invece la somma, l' $\lambda$ e varierebbe in base alla dimensione del gruppo, il che a sua volta altererebbe 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 calcolato sui termini del gruppo $G_i$, ricavato dai coefficienti originali. La normalizzazione ridimensiona quindi ciascun termine fermionico portandolo a grandezza unitaria prima della trasformazione di Jordan-Wigner, conservandone solo la fase e il segno nell'operatore evoluto. Ciò corrisponde a qDRIFT's, dove l’evoluzione dei termini normalizzati: le grandezze dei coefficienti determinano le probabilità di campionamento tramite $\lambda$, mentre l’evoluzione del circuito utilizza il termine normalizzato in modo che i coefficienti di grande valore non vengano conteggiati due volte.

Esempi di 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...

Genera gli insiemi di operatori per 100 circuiti casuali, ciascuno contenente 10 gruppi estratti dalla distribuzione ponderata. Il passo successivo trasforma questi insiemi nell’insieme dei circuiti di tipo “ qDRIFT ”.

Fase 2: Realizzare i circuiti " qDRIFT " (C++)

Nella versione in C++, questa fase costruisce i circuiti a partire dagli insiemi di operatori campionati. La transpilazione per il backend di destinazione avviene nella Fase 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 visualizza quindi le statistiche del circuito e gli indici dell'operatore del circuito 0, che qui sono stati omessi.

Per ciascun operatore campionato, crea un circuito di evoluzione di 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é ciascun gruppo campionato è stato normalizzato termine per termine prima della mappatura, questa evoluzione riguarda solo la fase e il segno di ciascun termine durante la rotazione; le grandezze contribuiscono esclusivamente attraverso la distribuzione campionaria. Lo stato iniziale di Hartree-Fock viene preparato applicando porte X ai qubit da 0 a 6 e da 10 a 16 (sette elettroni in ciascun settore di spin). Le istruzioni di ciascun circuito di evoluzione vengono aggiunte manualmente al circuito principale tramite l'API C++.

Passaggio 3: Esecuzione tramite " Qiskit primitives " (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 in base al nome, traspila ciascun circuito per quel backend e invia un lavoro per ogni circuito : 100 lavori Sampler distinti, ciascuno composto da 100 shot su un registro da 20 qubit. Quindi controlla lo stato di ogni processo fino al suo completamento e raccoglie le stringhe di bit da result[0]. Ogni processo utilizza tempo di elaborazione della QPU e rimane in coda separatamente, quindi è prevedibile che questa fase richieda molto più tempo rispetto a un singolo processo eseguito in batch.

I risultati relativi all'hardware riportati in questa pagina derivano da una serie di test ( ibm_pittsburgh 100 circuiti × 100 shot), in cui il valore è stato modificato backend_name rispetto a quello predefinito ibm_fez. I campioni salvati sono stati sottoposti a post-elaborazione utilizzando l'ordinamento corretto dei bit del determinante illustrato nella Fase 4. Il rumore dell'hardware varia da una corsa all'altra, quindi i conteggi e l'energia finale risulteranno 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

Modifica del backend: modifica backend_name in SqDRIFT.cpp per selezionare un dispositivo diverso; service.backends() elenca i dispositivi accessibili dal tuo account.

Fase 4: Elaborazione successiva e restituzione del risultato nel formato classico desiderato (C++)

Pubblica stringhe di bit selezionate e scrivi 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!

Seleziona stringhe di bit con il peso di Hamming corretto (n_alpha elettroni con spin verso l’alto, n_beta elettroni con spin verso il basso, derivati dal campo NELEC FCIDump), le converte in stringhe di interazione di configurazione (CI) utilizzando la simmetrizzazione di spin e le scrive in alphadets_from_sqd.txt per la successiva diagonalizzazione a basi selezionate.

Pulisci 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 codice C++ std::vector<QuantumCircuit> viene automaticamente liberato tramite RAII. Libera manualmente gli operatori dei qubit, gli operatori fermionici e i dati FCIDump.

Ora che abbiamo creato le basi per la proiezione del nostro hamiltoniano, possiamo procedere con il processo di diagonalizzazione per ottenere la stima dello stato fondamentale. Esegui i seguenti comandi dalla directory principale del progetto, docs/tutorials/assets/sqdrift/cpp.

Compilare il file binario diag SBD

CMake recupera SBD solo come codice sorgente in deps/sbd/ (la compilazione principale ne tralascia intenzionalmente la compilazione). L'app SBD utilizza un file di configurazione standard Makefile, e SBD stesso richiede MPI, OpenMP, e BLAS/LAPACK — pertanto il file di configurazione standard non è make sufficiente finché tali componenti non sono installati e il Configuration file di configurazione dell'app non è compatibile con la propria toolchain.

Il programma Makefile legge i parametri del compilatore e del linker dal file Configuration presente nella stessa directory. La versione fornita dall'upstream è ottimizzata per una specifica configurazione di macOS/Homebrew e non può essere compilata così com'è su un computer di serie — in particolare, non -I/opt/homebrew/opt/llvm/include contiene omp.h (che si trova in una directory delle risorse di Clang specifica per quella versione), e il sistema di Apple la clang++ rifiuta -fopenmp categoricamente:

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

La seguente ricetta è stata verificata e funziona 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

Quindi sostituire con deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization/Configuration :

# 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 compilazione riuscita stampa le due righe relative alla compilazione e al collegamento e crea un file binario diag accanto al file Makefile. Uno ld: warning: ignoring duplicate libraries: '-lomp' è innocuo.

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

Linux ( Ubuntu / Debian ) — non verificato

I seguenti prerequisiti riguardano MPI, l' OpenMP, e e BLAS/LAPACK. GCC lo supporta in modo nativo -fopenmp , quindi non è necessario alcun intervento sul 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

Eseguire il comando SqDRIFT per generare le stringhe CI

Se hai già completato i passaggi da ./SqDRIFT 1 a 4, il file esiste già alphadets_from_sqd.txt nella directory principale del progetto e puoi saltare questo passaggio. Altrimenti, eseguilo dalla directory principale del progetto:

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

Eseguire 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 l'hamiltoniano
• --adetfile: Percorso del file dei determinanti alfa (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 linea Energy = -107.5108722590187 rappresenta la stima dello stato fondamentale sul sottospazio campionato.

In che modo questo differisce dall'esempio " Python "

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

Entrambe le versioni utilizzano lo stesso hamiltoniano, raggruppano i suoi termini in base alla struttura elettronica, ponderano ciascun gruppo in base al coefficiente assoluto medio dei suoi termini, partono dallo stato di Hartree-Fock e campionano 10 gruppi e 100 shot per circuito.

Ogni differenza influisce in modo indipendente sul sottospazio campionato e, di conseguenza, sull'energia. Sia l'insieme più ampio che il secondo tempo di evoluzione, più lungo, ampliano le configurazioni esplorate dall'esempio di “ Python ”; il campionamento in più di un istante conferisce a “ SqDRIFT ” una copertura simile a quella di Krylov dello spazio a bassa energia. Il filtraggio dei termini diagonali modifica gli operatori che il campionatore può estrarre e i relativi pesi, pertanto le due esecuzioni campionano distribuzioni diverse.

La post-selezione è un filtro rigoroso, e un’immagine con un numero di particelle errato viene scartata immediatamente. Nell'esecuzione sopra descritta, che ha scartato la maggior parte dei dati, 1137 scatti su 10.000 hanno superato il controllo del peso di Hamming (7,7) (circa l'11%), raggruppandosi in soli 49 determinanti distinti di $\alpha$. Il recupero della configurazione utilizza invece le occupazioni orbitali medie ricavate da una precedente diagonalizzazione per invertire i bit e correggere tali stati, riutilizzando ciò che la post-selezione scarterebbe. Questo feedback lo rende intrinsecamente iterativo: recover_configurations in qiskit-addon-sqd-hpc richiede un avg_occupancies argomento che esiste solo dopo l'esecuzione di una diagonalizzazione.

Questa esecuzione in C++ produce un sottospazio di 49 determinanti di $\alpha$ e un’energia pari a Ha -107.5109 . Lo spazio completo dei determinanti di $\alpha$ per sette elettroni in 10 orbitali è costituito da sole stringhe “ $\binom{10}{7} = 120$ ”, quindi inserendo tutte e 120 nello stesso codice binario diag si ottiene l’energia FCI esatta Ha -107.6482 . La serie si discosta quindi di circa Ha 0.1373 dal valore esatto, coprendo 49 dei 120 determinanti di $\alpha$. Poiché l’SQD è variazionale nel sottospazio selezionato, un sottospazio più piccolo non può che aumentare l’energia; pertanto, il sottospazio incompleto rappresenta l’interpretazione più diretta del gap.

Questa verifica esatta dell'FCI è possibile solo perché l'esempio è di piccole dimensioni; non fa parte del flusso di lavoro SqDRIFT e, per sistemi di dimensioni tali da richiedere effettivamente l'uso dell'SQD, non esiste alcun riferimento di questo tipo. In questo caso serve solo a quantificare la porzione di spazio recuperata da questa esecuzione.

Il confronto con l'esempio " Python " non è controllato, pertanto una differenza non dovrebbe essere attribuita esclusivamente all'assenza del recupero della configurazione; si noti che il numero di circuiti, il tempo di evoluzione singolo, il filtraggio dei termini diagonali e il rumore hardware differiscono tutti contemporaneamente e contribuiscono tutti a modellare il sottospazio.

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