Algoritmo SqDRIFT per la stima dello stato fondamentale

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

Cerchi la versione C++?

Questo tutorial usa Python. Per l'implementazione C++, incluso il codice sorgente e le istruzioni di compilazione, consulta il tutorial SqDRIFT in C++.

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

Questo tutorial è presentato come notebook Python a scopo didattico.

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.

Questo tutorial utilizza Qiskit Fermions per creare i circuiti fermionici più naturali per l'algoritmo qDRIFT, seguito dall'uso di layout fermionico e passaggi di sintesi prima di inserire i circuiti nella pipeline tradizionale di Qiskit per l'esecuzione hardware.

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

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

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

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

# Qiskit Aer
from qiskit_aer import AerSimulator

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

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

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

Esempio con simulatore

Passo 1: Mappa gli input classici a un problema quantistico

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 nella cella successiva.

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 convergenti, 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 la cella 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 = "assets/sqdrift/fcidump_files/N2_sto_3g"

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

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

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 compilazione 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$$

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 compilazione 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 possono variare in modo in qualche misura 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 abbia più o meno di 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; il companion in C++ ne parla più in dettaglio, poiché effettua una postselezione anziché un recupero e differisce anche in quegli altri parametri.

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

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

print(len(canon.groups))

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

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

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

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

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

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

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

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

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

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

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

shots = 100

sampler = Sampler(mode=backend)

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

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

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

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

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

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

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

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

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

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

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

result_history = []

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

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

computed_energy = result.energy + nuclear_repulsion_energy

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

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

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

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

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

Prossimi passi

Raccomandazioni

Se hai trovato interessante questo lavoro, potrebbero interessarti anche i seguenti materiali:

• Diagonalizzazione quantistica di Krylov basata su campionamento di un modello reticolare fermionico - un tutorial correlato che usa circuiti di evoluzione temporale invece di un ansatz variazionale.
• Diagonalizzazione quantistica basata su campionamento di un'Hamiltoniana chimica - un tutorial su come costruire un circuito local unitary cluster Jastrow (LUCJ) per la simulazione della 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 fase di sviluppo, e questo tutorial è soggetto a cambiamenti in futuro in base all'evoluzione delle librerie usate.)