SQD per la stima energetica di un hamiltoniano chimico

In questa lezione, applicheremo l'SQD per stimare l'energia di stato fondamentale di una molecola.

In particolare, verranno affrontati i seguenti argomenti utilizzando l'approccio del pattern $4$ Qiskit:

1. Passo 1: mappare il problema in circuiti e operatori quantistici
   • Impostazione dell'hamiltoniana molecolare per $N_2$.
   • Spiega il cluster unitario locale Jastrow (LUCJ) ispirato alla chimica e compatibile con l'hardware [1]
2. Fase 2: Ottimizzazione per l'hardware di destinazione
   • Ottimizzare il numero di gate e il layout dell'ansatz per l'esecuzione hardware
3. Fase 3: Esecuzione sull'hardware di destinazione
   • Eseguire il circuito ottimizzato su una QPU reale per generare campioni del sottospazio.
4. Fase 4: Postelaborazione dei risultati
   • Introdurre il ciclo di recupero della configurazione autoconsistente [2]
     • Postelaborazione dell'intero set di campioni di bitstring, utilizzando la conoscenza preliminare del numero di particelle e dell'occupazione orbitale media calcolata nell'iterazione più recente.
     • Creare in modo probabilistico lotti di sottocampioni dalle stringhe di bit recuperate.
     • Proiettare e diagonalizzare l'hamiltoniana molecolare su ogni sottospazio campionato.
     • Salva l'energia minima del ground state trovata in tutti i batch e aggiorna l'occupazione avg degli orbitali.

Nel corso della lezione verranno utilizzati diversi pacchetti software.

• PySCF per definire la molecola e impostare l'hamiltoniana.
• ffsim per costruire l'ansatz LUCJ.
• Qiskit per il transpiling dell'ansatz per l'esecuzione hardware.
• Qiskit IBM Runtime per eseguire il circuito su una QPU e raccogliere campioni.
• Qiskit addon SQD recupero della configurazione e stima dell'energia dello stato di terra utilizzando la proiezione del sottospazio e la diagonalizzazione della matrice.

1. Mappare il problema su circuiti quantistici e operatori

Hamiltoniano molecolare

L'hamiltoniana molecolare assume la forma generica:

$\hat{a}^\dagger_{p\sigma}$ / $\hat{a}_{p\sigma}$ sono gli operatori di creazione/annichilazione fermionici associati all'elemento $p$ dell'insieme base e allo spin $\sigma$. $h_{pr}$ e $(pr|qs)$ sono gli integrali elettronici a uno e due corpi. Utilizzando pySCF, definiremo la molecola e calcoleremo gli integrali a uno e due corpi dell'hamiltoniana per l'insieme di basi 6-31g.

import warnings
import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# Specify molecule properties
open_shell = False
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],  # Two N atoms 1 angstrom apart
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)  # hcore: one-body integrals
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)  # eri: two-body integrals

# Compute exact energy for comparison
exact_energy = cas.run().e_tot

Output:

converged SCF energy = -108.835236570774
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

In questa lezione utilizzeremo la trasformazione Jordan-Wigner (JW) per mappare una funzione d'onda fermionica in una funzione d'onda qubit, in modo da poterla preparare con un circuito quantistico. La trasformazione JW mappa lo spazio di Fock dei fermioni in M orbitali spaziali sullo spazio di Hilbert di 2M qubit, cioè un orbitale spaziale viene diviso in due orbitali di spin, uno associato a un elettrone con spin up ( $\alpha$ ) e un altro con spin down ( $\beta$ ). Un orbitale di spin può essere occupato o non occupato. Di solito, quando si parla di numero di orbitali, si fa riferimento al numero di orbitali spaziali. Il numero di orbitali di spin sarà doppio. Nei circuiti quantistici, rappresenteremo ogni orbitale di spin con un qubit. Così, un insieme di qubit rappresenterà orbitali di spin-up o $\alpha$, e un altro insieme rappresenterà orbitali di spin-down o $\beta$. Ad esempio, la molecola $N_2$ per l'insieme di basi 6-31g ha orbitali spaziali $16$ (cioè, $16$ $\alpha$ + $16$ $\beta$ = $32$ orbitali di spin). Avremo quindi bisogno di un circuito quantistico di $32$ -qubit (potremmo avere bisogno di qubit ancilla aggiuntivi, come discusso in seguito). I qubit sono misurati in base al calcolo per generare bitstring, che rappresentano configurazioni elettroniche o determinanti (Slater). In questa lezione utilizzeremo i termini bitstring, configurazioni e determinanti in modo intercambiabile. Le stringhe di bit ci dicono l'occupazione degli elettroni negli orbitali di spin: un $1$ in una posizione di bit significa che l'orbitale di spin corrispondente è occupato, mentre un $0$ significa che l'orbitale di spin è vuoto. Dato che i problemi di struttura elettronica sono a conservazione particellare, deve essere occupato solo un numero fisso di orbitali di spin. La molecola $N_2$ ha $5$ elettroni con spin-up ( $\alpha$ ) e $5$ con spin-down ( $\beta$ ). Pertanto, qualsiasi bitstring che rappresenti gli orbitali $\alpha$ e $\beta$ deve avere cinque $1\text{s}$ ciascuno per la molecola $N_2$.

1.1 Circuito quantistico per la generazione di campioni: l'approccio LUCJ

In questa lezione, utilizzeremo l'ipotesi di cluster unitario locale accoppiato di Jastrow (LUCJ) \ [1] per la preparazione dello stato quantistico e il successivo campionamento. In primo luogo, spiegheremo i diversi elementi costitutivi dell'ansatz UCJ completo e le approssimazioni fatte nella versione locale di esso. Successivamente, utilizzando il pacchetto ffsim, costruiremo l'ansatz LUCJ e lo ottimizzeremo utilizzando il transpiler Qiskit per l'esecuzione hardware.

L'ansatz di UCJ ha la forma seguente (per un prodotto di $L$ strati o ripetizioni dell'operatore UCJ)

dove, $\vert \Phi_{0} \rangle$ è uno stato di riferimento, tipicamente assunto come stato Hartree-Fock (HF). Poiché lo stato Hartree-Fock è definito come quello in cui gli orbitali più bassi sono occupati, la preparazione dello stato HF prevede l'applicazione di porte X per impostare su uno i qubit corrispondenti agli orbitali occupati. Ad esempio, il blocco di preparazione dello stato HF per 4 orbitali spaziali e 2 up- e 2 down-spin può essere il seguente:

Una singola ripetizione dell'operatore UCJ ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ consiste in un'evoluzione diagonale di Coulomb ( $e^{iJ^{(\mu)}}$ ) intramezzata da rotazioni orbitali ( $e^{K^{(\mu)}}$ e $e^{-K^{(\mu)}}$ ).

I blocchi di rotazione orbitale lavorano su una singola specie di spin ( $\alpha$ (up-spin)/ $\beta$ (down-spin)). Per ogni specie di elettrone, la rotazione degli orbitali consiste in uno strato di porte a singolo qubit $R_{z}$ seguito da una sequenza di porte di rotazione di Given a 2 qubit (porte $XX + YY$ ).

Le porte a 2 qubit agiscono su orbitali di spin adiacenti (qubit vicini) e sono quindi implementabili sulle QPU IBM® senza bisogno di porte SWAP.

L'operatore $e^{iJ^{(\mu)}}$, noto anche come operatore diagonale di Coulomb, è composto da tre blocchi. Due di essi lavorano sugli stessi settori di spin ( $e^{iJ_{\alpha \alpha}^{(\mu)}}$ e $e^{iJ_{\beta \beta}^{(\mu)}}$ ), mentre uno lavora tra due settori di spin ( $e^{iJ_{\alpha \beta}^{(\mu)}}$ ).

Tutti i blocchi di $e^{iJ^{(\mu)}}$ sono costituiti da porte numeriche $U_{nn}(\phi)$ [1]. Una porta $U_{nn}(\phi)$ può essere ulteriormente suddivisa in una porta $R_{ZZ}(\frac{\phi}{2})$ seguita da due porte $Rz(-\frac{\phi}{2})$ a singolo qubit che agiscono su due qubit separati.

I componenti same-spin ( $J_{\alpha \alpha}$ e $J_{\beta \beta}$ ) hanno $U_{nn}$ porte tra tutte le possibili coppie di qubit. Tuttavia, poiché le QPU superconduttrici hanno una connettività restrittiva, i qubit devono essere scambiati per realizzare porte tra qubit non adiacenti.

Ad esempio, si consideri il seguente blocco $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (o $e^{iJ_{\beta \beta}^{(\mu)}}$ ) per gli orbitali spaziali $N = 4$. Per una connettività lineare dei qubit, le ultime tre porte non sono direttamente implementabili, poiché funzionano tra qubit non adiacenti (ad esempio, Q0 e Q2 non sono direttamente connessi). Pertanto, abbiamo bisogno di porte SWAP per renderle adiacenti (la figura seguente mostra un esempio con $3$ porte SWAP).

Poi, $J_{\alpha \beta}$ implementa i gate tra gli stessi orbitali indicizzati di settori di spin diversi (ad esempio, tra $0\alpha$ e $0\beta$ ). Analogamente, se i qubit non sono fisicamente adiacenti su una QPU, anche questi gate richiederanno degli SWAP.

Dalla discussione precedente, l'ansatz UCJ incontra alcuni ostacoli per l'esecuzione HW, in quanto necessita di porte SWAP a causa delle interazioni tra qubit non adiacenti. La variante locale dell'ansatz UCJ, LUCJ, affronta questa sfida eliminando alcuni $U_{nn}$ dall'operatore di Coulomb diagonale.

Negli stessi blocchi di specie di elettroni, $J_{\alpha \alpha}$ e $J_{\beta \beta}$ ), manteniamo solo le porte di $U_{nn}$ compatibili con la connettività nearest-neighbor e rimuoviamo le porte tra qubit non adiacenti nella versione LUCJ. La figura seguente mostra il blocco LUCJ dopo la rimozione delle porte non adiacenti.

Successivamente, la versione LUCJ del blocco $J_{\alpha \beta}$ che lavora tra diverse specie di elettroni può assumere forme diverse in base alla topologia del dispositivo.

Anche in questo caso, la versione LUCJ elimina le porte non compatibili. La figura seguente mostra le varianti del blocco $J_{\alpha \beta}$ per diverse topologie di qubit, tra cui griglia, esagono, heavy-hex e lineare.

• Quadrato : possiamo avere porte $U_{nn}$ tra tutti gli orbitali $\alpha$ e $\beta$ senza alcuno SWAP, e quindi non è necessario rimuovere alcuna porta $U_{nn}$.
• Pesantezza dell'esagono : Le interazioni $\alpha$ - $\beta$ sono mantenute tra ogni $4$ -esimo orbitale di spin indicizzato (come lo 0°, il 4° e l'8°) e sono mediate da ancilla, cioè abbiamo bisogno di qubit ancilla tra le catene lineari che rappresentano gli orbitali $\alpha$ e $\beta$. Questa disposizione richiede un numero limitato di SWAP.
• Esagonale : Ogni altro orbitale, come gli orbitali indicizzati 0, 2 e 4, diventano vicini di casa quando $\alpha$ e $\beta$ sono disposti in due catene lineari adiacenti.
• Lineare : Sono collegati solo un orbitale $\alpha$ e uno $\beta$, il che significa che il blocco $J_{\alpha \beta}$ avrà un solo gate.

Schemi di

Se da un lato l'eliminazione dei gate dall'ansatz UCJ per costruire la versione LUCJ la rende più compatibile con l'HW, dall'altro l'ansatz perde un po' di espressività. Pertanto, potrebbero essere necessarie più ripetizioni ( $L$ ) dell'operatore UCJ modificato quando si utilizza l'ansatz LUCJ.

1.2 Inizializzazione dell'approccio LUCJ

Il LUCJ è un ansatz parametrizzato e dobbiamo inizializzare i parametri prima dell'esecuzione hardware. Un modo per inizializzare l'ansatz è quello di utilizzare le ampiezze t1 e t2 del metodo classico di cluster accoppiati singoli e doppi (CCSD), dove le ampiezze t1 sono il coefficiente degli operatori di eccitazione singola e le ampiezze t2 sono per gli operatori di eccitazione doppia.

Si noti che, mentre l'inizializzazione dell'ansatz LUCJ con le ampiezze t1 e t2 genera risultati decenti, i parametri dell'ansatz potrebbero richiedere un'ulteriore ottimizzazione.

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
)
ccsd.run()

t1 = ccsd.t1
t2 = ccsd.t2

Output:

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 Costruzione dell'approccio LUCJ utilizzando ffsim

Utilizzeremo il pacchetto ffsim per creare e inizializzare l'ansatz con le ampiezze t1 e t2 calcolate in precedenza. Poiché la nostra molecola ha uno stato Hartree-Fock a guscio chiuso, utilizzeremo la variante con equilibrio di spin dell'ansatz UCJ, UCJOpSpinBalanced.

Poiché l'hardware di IBM ha una topologia heavy-hex, adotteremo lo schema a zig-zag utilizzato in [1] e spiegato sopra per le interazioni tra qubit. In questo schema, gli orbitali (qubit) con lo stesso spin sono collegati con una topologia a linee (cerchi rossi e blu). A causa della topologia heavy-hex, gli orbitali di spin diversi hanno connessioni tra ogni quarto orbitale, cioè lo 0°, il 4°, l'8° e così via (cerchi viola).

import ffsim
from qiskit import QuantumCircuit, QuantumRegister

n_reps = 2
alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]
alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]

ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),
)

nelec = (num_elec_a, num_elec_b)

# create an empty quantum circuit
qubits = QuantumRegister(2 * num_orbitals, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()
# circuit.decompose().draw("mpl", scale=0.5, fold=-1)

L'ipotesi LUCJ con strati ripetuti può essere ottimizzata unendo alcuni blocchi adiacenti. Si consideri il caso di n_reps=2. I due blocchi di rotazione orbitale al centro possono essere uniti in un unico blocco di rotazione orbitale. Il pacchetto ffsim ha un gestore di passaggi chiamato ffsim.qiskit.PRE_INIT per ottimizzare il circuito unendo tali blocchi adiacenti.

2. Ottimizzazione per l'hardware di destinazione

Per prima cosa, recuperiamo un backend di nostra scelta. Ottimizzeremo il nostro circuito per il backend e poi eseguiremo il circuito ottimizzato sullo stesso backend per generare campioni per il sottospazio.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

Per ottimizzare il surrogato e renderlo compatibile con l'hardware, si consiglia di eseguire le seguenti operazioni.

• Selezionare i qubit fisici (initial_layout) dall'hardware di destinazione che aderisce allo schema a zig-zag (due catene lineari con in mezzo un qubit ancilla) descritto sopra. La disposizione dei qubit secondo questo schema porta a un circuito efficiente e compatibile con l'hardware, con meno porte.
• Generare uno staged pass manager utilizzando la funzione generate_preset_pass_manager di Qiskit con la scelta di backend e initial_layout.
• Impostate lo stage pre_init del vostro gestore di staged pass su ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT include i passaggi del transpiler Qiskit che scompongono i gate in rotazioni orbitali e poi uniscono le rotazioni orbitali, ottenendo un minor numero di gate nel circuito finale.
• Eseguire il pass manager sul circuito.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]
spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]

initial_layout = spin_a_layout + spin_b_layout

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend, initial_layout=initial_layout
)

# without PRE_INIT passes
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}")

# with PRE_INIT passes
# We will use the circuit generated by this pass manager for hardware execution
pass_manager.pre_init = ffsim.qiskit.PRE_INIT
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}")

Output:

Gate counts (w/o pre-init passes): OrderedDict({'rz': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})

3. Eseguire sull'hardware di destinazione

Dopo aver ottimizzato il circuito per l'esecuzione hardware, siamo pronti a eseguirlo sull'hardware di destinazione e a raccogliere campioni per la stima dell'energia dello stato fondamentale. Dato che abbiamo un solo circuito, useremo il ed qiskit-ibm-runtimeModalità di esecuzione del processo eseguiremo il nostro circuito.

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
sampler.options.dynamical_decoupling.enable = True

job = sampler.run([isa_circuit], shots=10_000)  # Takes approximately 5sec of QPU time

# Run cell after IQX job completion
primitive_result = job.result()
pub_result = primitive_result[0]
counts = pub_result.data.meas.get_counts()

4. Risultati post-elaborazione

La parte di post-elaborazione del flusso di lavoro SQD può essere riassunta con il seguente diagramma.

Il campionamento dell'ansatz LUCJ nella base di calcolo genera un pool di configurazioni rumorose $\tilde{\mathcal{\chi}}$, che vengono utilizzate nella routine di post-elaborazione. Si tratta di un metodo chiamato (i dettagli saranno discussi più avanti) recupero della configurazione per correggere probabilisticamente le configurazioni con numeri di elettroni errati. Le configurazioni solo con numeri di elettroni corretti $\tilde{\mathcal{\chi}}_{R}$ vengono quindi sottocampionate e distribuite in lotti multipli in base alla frequenza di comparsa di ciascuna configurazione unica. Ogni gruppo di campioni definisce un sottospazio ( $\mathcal{S^{(k)}}$ ). Quindi, l'hamiltoniana molecolare, $H$, viene proiettata sui sottospazi:

Ogni Hamiltoniana proiettata $H_{\mathcal{S}^{(k)}}$ viene quindi inserita in un Eigensolver, dove viene diagonalizzata per calcolare autovalori e autovettori e ricostruire un autostato. In questa lezione, proiettiamo e diagonalizziamo l'hamiltoniana usando il pacchetto qiskit-addon-sqd che utilizza il metodo di Davidson di PySCF per la diagonalizzazione.

Raccogliamo quindi l'autovalore più basso (energia) dai lotti e calcoliamo anche l'occupazione orbitale media, $\text{n}$. Le informazioni sull'occupazione media vengono utilizzate nella fase di recupero della configurazione per correggere probabilisticamente le configurazioni di disturbo.

Successivamente, spieghiamo in dettaglio il ciclo di recupero della configurazione autoconsistente e mostriamo esempi concreti di codice per implementare i passaggi sopra citati per stimare l'energia di stato fondamentale dell'hamiltoniana di $N_2$.

4.1 Ripristino della configurazione: panoramica

Ogni bit di una stringa di bit (determinante di Slater) rappresenta un orbitale di spin. La metà destra di una stringa di bit rappresenta gli orbitali di spin-up, mentre la metà sinistra rappresenta gli orbitali di spin-down. Un 1 significa che l'orbitale è occupato da un elettrone, mentre un 0 significa che l'orbitale è vuoto. Conosciamo a priori il numero corretto di particelle (sia l'elettrone up-spin che l'elettrone down-spin). Supponiamo di avere un determinante $x$ con $N_x$ elettroni (cioè, ci sono $N_x$ numeri di $1$ s nella stringa di bit). Il numero corretto di particelle è $N$. Se $N_x \neq N$, allora sappiamo che la stringa di bit è corrotta dal rumore. La routine di configurazione autoconsistente tenta di correggere la stringa di bit ribaltando in modo probabilistico i bit $|N_x - N|$ sfruttando le informazioni sull'occupazione orbitale media. L'occupazione media degli orbitali ( $n$ ) indica la probabilità che un orbitale sia occupato da un elettrone. Se $N_x < N$, abbiamo meno elettroni e dobbiamo trasformare alcuni $0$ s in $1$ s e viceversa.

La probabilità di capovolgimento può essere $|x[i] - avg\_occupancy[i]|$ per i-esimo orbitale di spin. In [2], gli autori hanno utilizzato una probabilità ponderata di ribaltamento utilizzando la funzione ReLU modificata.

Qui $h$ definisce la posizione dell'"angolo" della funzione ReLU e il parametro $\delta$ definisce il valore della funzione ReLU nell'angolo. Per $\delta = 0$, $w$ diventa la vera funzione ReLU, e per $\delta >0$, diventa ReLU* modificata*. Nel documento, gli autori hanno utilizzato $\delta = 0.01$ e $h =$ numero di particelle alfa (o beta)/numero di orbitali di spin alfa (o beta) $= N/M$ (fattore di riempimento).

L'occupazione orbitale media ( $n$ ) non è nota a priori. La prima iterazione della stima dello stato fondamentale inizia con configurazioni che presentano solo numeri di particelle corretti in entrambe le specie di spin. Dopo la prima iterazione, abbiamo una stima dello stato fondamentale e, utilizzando questa stima, possiamo costruire la prima ipotesi di $n$. Questa ipotesi di $n$ viene utilizzata per recuperare le configurazioni, eseguire l'iterazione successiva di stima dello stato fondamentale e perfezionare in modo autoconsistente l'ipotesi di $n$. Il processo si ripete finché non viene soddisfatto un criterio di arresto.

Si consideri il seguente esempio per $N = 2$ e $x = |1000\rangle$ ( $N_x = 1$ ). Dobbiamo invertire uno dei 0s in 1 per correggere il numero di particelle, e le scelte sono 1100, 1010, e 1001. In base alla probabilità di capovolgimento, una delle scelte sarà selezionata come configurazione recuperata (o la stringa di bit con il numero corretto di particelle).

Supponiamo che nella prima iterazione vengano eseguiti due lotti e che gli stati a terra stimati da questi siano:

Utilizzando gli stati base computazionali e le loro ampiezze, possiamo calcolare la probabilità di occupazione degli elettroni (in breve le occupazioni ) per spin-orbitale (qubit) (si noti che probabilità = |ampiezza| $^2$ ). Qui di seguito vengono riportate le occupazioni dei qubit per ogni stringa di bit che compare nello stato fondamentale stimato e viene calcolata l'occupazione totale degli orbitali per un gruppo. Si noti che, come da convenzione di ordinamento di Qiskit, il bit più a destra rappresenta qubit-0 ( Q0 ), e il bit più a sinistra rappresenta Q3.

Occupazione ( Batch0 ):

Occupazione ( Batch1 )

Occupazione (media dei lotti)

Utilizzando l'occupazione orbitale media calcolata sopra, possiamo trovare le probabilità di flip per i diversi orbitali nella configurazione $x = \vert 1000 \rangle$. Poiché l'orbitale rappresentato da Q3 è già occupato e non deve essere capovolto, impostiamo la sua p(flip) a $0$. Per i restanti orbitali, che non sono occupati, la probabilità di flip è $\vert x[i] - \text{n}[i] \vert$. Oltre a p(flip), calcoliamo anche il peso della probabilità associata al flip utilizzando la funzione ReLU modificata descritta in precedenza.

Probabilità di capovolgimento ( $x = \vert 1000 \rangle$, $\delta = 0.01$, $h = N/M = 2/4 = 0.50$ )

Infine, utilizzando le probabilità ponderate di cui sopra, possiamo capovolgere uno degli orbitali non occupati Q2, Q1 e Q0. In base ai valori sopra indicati, Q0 sarà probabilmente invertito e una possibile configurazione recuperata può essere $\vert \text{1001} \rangle$.

Il processo completo di ripristino della configurazione autoconsistente può essere riassunto come segue:

Prima iterazione: Supponiamo che le bitstring (configurazioni o determinanti di Slater) generate dal computer quantistico formino un insieme $\widetilde{\chi}$, che comprende sia le configurazioni con un numero corretto ( $\widetilde{\chi}_{correct}$ ) sia quelle con un numero errato ( $\widetilde{\chi}_{incorrect}$ ) di particelle in ogni settore di spin.

1. Le configurazioni di ( $\widetilde{\chi}_{correct}$ ) vengono campionate in modo casuale per creare lotti $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ di vettori per la proiezione del sottospazio. Il numero di lotti e di campioni in ciascun lotto sono parametri definiti dall'utente. Più grande è il numero di campioni in ogni lotto, più grande è la dimensione del sottospazio e più impegnativa dal punto di vista computazionale diventa la diagonalizzazione. D'altra parte, un numero troppo basso di campioni può perdere i vettori di supporto dello stato a terra e portare a una stima errata.
2. Eseguire il solutore di autostati (cioè proiezione sul sottospazio e diagonalizzazione) sui batch e ottenere autostati approssimati. $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. A partire dagli autostati approssimati si costruisce la prima ipotesi per $n$.

Iterazioni successive:

1. Utilizzando $n$ si correggono le configurazioni con numero di particelle errato in $\widetilde{\chi}_{incorrect}$. Supponiamo di chiamarle $\widetilde{\chi}_{correct\_new}$. Quindi, $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ forma il nuovo insieme di configurazioni con numeri di particelle corretti.
2. $\widetilde{\chi}_{R}$ viene campionata per creare lotti $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. Il solutore di autostati viene eseguito con nuovi lotti e genera nuove stime degli stati fondamentali $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. Dagli autostati approssimati si costruiscono ipotesi raffinate per $n$.
5. Se il criterio di arresto non è soddisfatto, si torna al passo 2.1.

4.2 Stima dello stato fondamentale

Per prima cosa, trasformeremo i conteggi in una matrice di bitstring e in un array di probabilità per la post-elaborazione.

Ogni riga della matrice rappresenta una stringa di bit unica. Poiché i qubit sono indicizzati a partire dalla destra di una stringa di bit in Qiskit, la colonna 0 rappresenta il qubit N-1, e la colonna N-1 rappresenta il qubit 0, dove N è il numero di qubit.

Gli orbitali alfa sono rappresentati nell'intervallo di indici di colonna (N, N/2] (metà destra), mentre gli orbitali beta sono rappresentati nell'intervallo di colonne (N/2, 0] (metà sinistra).

from qiskit_addon_sqd.counts import counts_to_arrays

# Convert counts into bitstring and probability arrays
bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)

Ci sono alcune opzioni controllate dall'utente che sono importanti per questa tecnica:

• iterations: Numero di iterazioni di recupero della configurazione autoconsistente
• n_batches: Numero di batch di configurazioni utilizzati dalle diverse chiamate al risolutore di autovalori
• samples_per_batch: Numero di configurazioni uniche da includere in ogni lotto
• max_davidson_cycles: Numero massimo di cicli Davidson eseguiti da ciascun eigensolver

import numpy as np
from qiskit_addon_sqd.configuration_recovery import recover_configurations
from qiskit_addon_sqd.fermion import (
    bitstring_matrix_to_ci_strs,
    solve_fermion,
)
from qiskit_addon_sqd.subsampling import postselect_and_subsample

rng = np.random.default_rng(24)
# SQD options
iterations = 5

# Eigenstate solver options
n_batches = 5
samples_per_batch = 500
max_davidson_cycles = 300

# Self-consistent configuration recovery loop
e_hist = np.zeros((iterations, n_batches))  # energy history
s_hist = np.zeros((iterations, n_batches))  # spin history
occupancy_hist = []
avg_occupancy = None
for i in range(iterations):
    print(f"Starting configuration recovery iteration {i}")
    # On the first iteration, we have no orbital occupancy information from the
    # solver, so we begin with the full set of noisy configurations.
    if avg_occupancy is None:
        bs_mat_tmp = bitstring_matrix_full
        probs_arr_tmp = probs_arr_full

    # If we have average orbital occupancy information, we use it to refine
    # the full set of noisy configurations.
    else:
        bs_mat_tmp, probs_arr_tmp = recover_configurations(
            bitstring_matrix_full,
            probs_arr_full,
            avg_occupancy,
            num_elec_a,
            num_elec_b,
            rand_seed=rng,
        )

    # Create batches of subsamples. We postselect here to remove configurations
    # with incorrect hamming weight during iteration 0, since no config recovery was performed.
    batches = postselect_and_subsample(
        bs_mat_tmp,
        probs_arr_tmp,
        hamming_right=num_elec_a,
        hamming_left=num_elec_b,
        samples_per_batch=samples_per_batch,
        num_batches=n_batches,
        rand_seed=rng,
    )

    # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.
    e_tmp = np.zeros(n_batches)
    s_tmp = np.zeros(n_batches)
    occs_tmp = []
    coeffs = []
    for j in range(n_batches):
        strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])
        print(f"  Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}")
        energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(
            batches[j],
            hcore,
            eri,
            open_shell=open_shell,
            spin_sq=spin_sq,
            max_cycle=max_davidson_cycles,
        )
        energy_sci += nuclear_repulsion_energy
        e_tmp[j] = energy_sci
        s_tmp[j] = spin
        occs_tmp.append(avg_occs)
        coeffs.append(coeffs_sci)

    # Combine batch results
    avg_occupancy = tuple(np.mean(occs_tmp, axis=0))

    # Track optimization history
    e_hist[i, :] = e_tmp
    s_hist[i, :] = s_tmp
    occupancy_hist.append(avg_occupancy)

Output:

Starting configuration recovery iteration 0
  Batch 0 subspace dimension: 21609
  Batch 1 subspace dimension: 21609
  Batch 2 subspace dimension: 21609
  Batch 3 subspace dimension: 21609
  Batch 4 subspace dimension: 21609
Starting configuration recovery iteration 1
  Batch 0 subspace dimension: 609961
  Batch 1 subspace dimension: 616225
  Batch 2 subspace dimension: 627264
  Batch 3 subspace dimension: 633616
  Batch 4 subspace dimension: 624100
Starting configuration recovery iteration 2
  Batch 0 subspace dimension: 564001
  Batch 1 subspace dimension: 605284
  Batch 2 subspace dimension: 582169
  Batch 3 subspace dimension: 559504
  Batch 4 subspace dimension: 591361
Starting configuration recovery iteration 3
  Batch 0 subspace dimension: 550564
  Batch 1 subspace dimension: 549081
  Batch 2 subspace dimension: 531441
  Batch 3 subspace dimension: 527076
  Batch 4 subspace dimension: 531441
Starting configuration recovery iteration 4
  Batch 0 subspace dimension: 544644
  Batch 1 subspace dimension: 580644
  Batch 2 subspace dimension: 527076
  Batch 3 subspace dimension: 531441
  Batch 4 subspace dimension: 537289

4.3 Discussione dei risultati

Il primo grafico mostra che dopo alcune iterazioni si stima l'energia dello stato fondamentale entro ~24 mH (l'accuratezza chimica è tipicamente accettata come 1 kcal/mol $\approx$ 1.6 mH ). Il secondo grafico mostra l'occupazione media di ciascun orbitale spaziale dopo l'iterazione finale. Possiamo notare che sia gli elettroni con spin-up che quelli con spin-down occupano i primi cinque orbitali con alta probabilità nelle nostre soluzioni.

Sebbene l'energia di stato fondamentale stimata sia discreta, non rientra nel limite di precisione chimica ( $\pm \approx 1.6$ mH ). Questo divario può essere attribuito alla piccola dimensione del sottospazio utilizzata in precedenza per la proiezione e la diagonalizzazione. Poiché abbiamo usato samples_per_batch=500, il sottospazio è spaziato da max $500$ vettori, che sono vettori mancanti dal supporto dello stato fondamentale. L'aumento del parametro samples_per_batch dovrebbe migliorare l'accuratezza a scapito di risorse di calcolo più classiche e del tempo di esecuzione.

# Data for energies plot
x1 = range(iterations)
min_e = [np.min(e) for e in e_hist]
e_diff = [abs(e - exact_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]
x2 = range(len(y2))

import matplotlib.pyplot as plt

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-6)
axs[0].axhline(
    y=chem_accuracy, color="#BF5700", linestyle="--", label="chemical accuracy"
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

print(f"Exact energy: {exact_energy:.5f} Ha")
print(f"SQD energy: {min_e[-1]:.5f} Ha")
print(f"Absolute error: {e_diff[-1]:.5f} Ha")
plt.tight_layout()
plt.show()

Output:

Exact energy: -109.04667 Ha
SQD energy: -109.02234 Ha
Absolute error: 0.02434 Ha

Esercizio per il lettore

Aumentare progressivamente il parametro samples_per_batch (ad esempio, da $1000$ a $10000$ con un passo di $1000$; consentito dalla memoria del computer) e confrontare le energie di stato fondamentale stimate.

Riferimenti

[1] M. Motta et al, "Un ponte tra l'intuizione fisica e l'efficienza hardware per gli stati elettronici correlati: il cluster unitario locale di Jastrow ansatz per la struttura elettronica" (2023). Chimica. Sci, 2023, 14, 11213.

[2] J. Robledo-Moreno e altri, "Chimica oltre le soluzioni esatte su un supercomputer incentrato sui quanti" (2024). arXiv:quant-ph/2405.05068.