SQD per la stima dell'energia di un Hamiltoniano di chimica

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

In particolare, discuteremo i seguenti argomenti seguendo l'approccio a $4$ passi del pattern Qiskit:

1. Passo 1: Mappare il problema su circuiti e operatori quantistici
   • Configurare l'Hamiltoniano molecolare per $N_2$.
   • Spiegare il locale unitary cluster Jastrow (LUCJ) [1], ispirato alla chimica e compatibile con l'hardware
2. Passo 2: Ottimizzare per l'hardware di destinazione
   • Ottimizzare il conteggio dei gate e il layout dell'ansatz per l'esecuzione su hardware
3. Passo 3: Eseguire sull'hardware di destinazione
   • Eseguire il circuito ottimizzato su una vera QPU per generare campioni del sottospazio
4. Passo 4: Post-elaborare i risultati
   • Introdurre il ciclo di recupero della configurazione auto-consistente [2]
     • Post-elaborare l'insieme completo di campioni di bitstring, usando la conoscenza a priori del numero di particelle e l'occupazione orbitale media calcolata all'iterazione più recente.
     • Creare probabilisticamente batch di sotto-campioni dalle bitstring recuperate.
     • Proiettare e diagonalizzare l'Hamiltoniano molecolare su ogni sottospazio campionato.
     • Salvare l'energia minima dello stato fondamentale trovata tra tutti i batch e aggiornare l'occupazione orbitale media.

Nel corso della lezione useremo diversi pacchetti software.

• PySCF per definire la molecola e configurare l'Hamiltoniano.
• Il pacchetto ffsim per costruire l'ansatz LUCJ.
• Qiskit per il transpiling dell'ansatz per l'esecuzione su hardware.
• Qiskit IBM Runtime per eseguire il circuito su una QPU e raccogliere i campioni.
• Qiskit addon SQD per il recupero della configurazione e la stima dell'energia dello stato fondamentale tramite proiezione sul sottospazio e diagonalizzazione matriciale.

1. Mappare il problema su circuiti e operatori quantistici

Hamiltoniano molecolare

Un Hamiltoniano molecolare ha la forma generica:

$$\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \sum_{ \substack{prqs\\\sigma\tau} } \frac{(pr|qs)}{2} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma}$$

$\hat{a}^\dagger_{p\sigma}$/$\hat{a}_{p\sigma}$ sono gli operatori fermonici di creazione/annichilazione associati al $p$-esimo elemento del set di base e allo spin $\sigma$. $h_{pr}$ e $(pr|qs)$ sono gli integrali elettronici a uno e due corpi. Usando pySCF, definiremo la molecola e calcoleremo gli integrali a uno e due corpi dell'Hamiltoniano per il set di base 6-31g.

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

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

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

In questa lezione useremo la trasformazione di Jordan-Wigner (JW) per mappare una funzione d'onda fermionica in una funzione d'onda di qubit, così 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: in pratica, un orbitale spaziale viene suddiviso in due spin-orbitali, uno associato a un elettrone con spin su ($\alpha$) e l'altro con spin giù ($\beta$). Uno spin-orbitale può essere occupato o non occupato. Di norma, quando si parla di numero di orbitali si intende il numero di orbitali spaziali; il numero di spin-orbitali sarà il doppio. Nei circuiti quantistici, ogni spin-orbitale è rappresentato da un qubit. Quindi un insieme di qubit rappresenterà gli orbitali spin-up o $\alpha$, e un altro insieme gli orbitali spin-down o $\beta$. Ad esempio, la molecola $N_2$ per il set di base 6-31g ha $16$ orbitali spaziali (cioè $16$ $\alpha$ + $16$ $\beta$ = $32$ spin-orbitali). Sono quindi necessari $32$ qubit (più eventuali qubit ancilla, come discusso in seguito). I qubit vengono misurati nella base computazionale per generare bitstring, che rappresentano configurazioni elettroniche o determinanti (di Slater). In questa lezione useremo i termini bitstring, configurazioni e determinanti in modo intercambiabile. Le bitstring indicano l'occupazione degli elettroni negli spin-orbitali: un 1 in una posizione significa che il corrispondente spin-orbitale è occupato, mentre uno 0 significa che è vuoto. Poiché i problemi di struttura elettronica conservano il numero di particelle, solo un numero fisso di spin-orbitali deve essere occupato. La molecola $N_2$ ha $5$ elettroni con spin-up ($\alpha$) e $5$ con spin-down ($\beta$). Pertanto, ogni bitstring che rappresenta gli orbitali $\alpha$ e $\beta$ deve avere cinque 1 per ciascun tipo per la molecola $N_2$.

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

In questa lezione useremo l'ansatz local unitary coupled cluster Jastrow (LUCJ) \[1\] per la preparazione dello stato quantistico e il successivo campionamento. Prima spiegheremo i diversi blocchi costitutivi dell'ansatz UCJ completo e le approssimazioni introdotte nella versione locale. Poi, usando il pacchetto ffsim, costruiremo l'ansatz LUCJ e lo ottimizzeremo con il transpiler di Qiskit per l'esecuzione su hardware.

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

$$|\psi\rangle = \prod_{\mu=1}^{L}{(e^{K^{\mu}} \times {e^{iJ^{\mu}}} \times {e^{-K^{\mu}}})} |\Phi_{0}\rangle$$

dove $\vert \Phi_{0} \rangle$ è uno stato di riferimento, tipicamente lo stato Hartree-Fock (HF). Poiché lo stato Hartree-Fock è definito come lo stato con gli orbitali a numerazione più bassa occupati, la sua preparazione richiede l'applicazione di gate X sui qubit corrispondenti agli orbitali occupati, per portarli a uno. Ad esempio, il blocco di preparazione dello stato HF per 4 orbitali spaziali con 2 elettroni spin-up e 2 spin-down può essere simile al seguente: Una singola ripetizione dell'operatore UCJ ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ consiste in un'evoluzione di Coulomb diagonale ($e^{iJ^{(\mu)}}$) affiancata da rotazioni orbitali ($e^{K^{(\mu)}}$ e $e^{-K^{(\mu)}}$). I blocchi di rotazione orbitale agiscono su una singola specie di spin ($\alpha$ (spin-up)/$\beta$ (spin-down)). Per ogni specie elettronica, la rotazione orbitale è composta da uno strato di gate $R_{z}$ a singolo qubit seguito da una sequenza di gate di rotazione di Givens a 2 qubit (gate $XX + YY$).

I gate a 2 qubit agiscono su spin-orbitali adiacenti (qubit vicini più prossimi) e sono quindi implementabili su QPU IBM® senza la necessità di gate SWAP. $e^{iJ^{(\mu)}}$, noto anche come operatore di Coulomb diagonale, è composto da tre blocchi. Due di essi agiscono sui settori dello stesso spin ($e^{iJ_{\alpha \alpha}^{(\mu)}}$ e $e^{iJ_{\beta \beta}^{(\mu)}}$), e uno agisce tra i due settori di spin ($e^{iJ_{\alpha \beta}^{(\mu)}}$).

Tutti i blocchi in $e^{iJ^{(\mu)}}$ sono composti da gate numero-numero $U_{nn}(\phi)$ [1]. Un gate $U_{nn}(\phi)$ può essere ulteriormente scomposto in un gate $R_{ZZ}(\frac{\phi}{2})$ seguito da due gate $Rz(-\frac{\phi}{2})$ a singolo qubit che agiscono su due qubit separati.

Le componenti dello stesso spin ($J_{\alpha \alpha}$ e $J_{\beta \beta}$) hanno gate $U_{nn}$ tra tutte le coppie possibili di qubit. Tuttavia, poiché le QPU superconduttive hanno una connettività limitata, i qubit devono essere scambiati (swapped) per realizzare gate tra qubit non adiacenti.

Ad esempio, considera il seguente blocco $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (o $e^{iJ_{\beta \beta}^{(\mu)}}$) per $N = 4$ orbitali spaziali. Per una connettività lineare dei qubit, gli ultimi tre gate non sono direttamente implementabili poiché agiscono tra qubit non adiacenti (ad esempio, Q0 e Q2 non sono direttamente connessi). È quindi necessario usare gate SWAP per renderli adiacenti (la figura seguente mostra un esempio con $3$ gate SWAP). Il blocco $J_{\alpha \beta}$ implementa gate tra orbitali con lo stesso indice ma di specie di spin diverse (ad esempio, tra $0\alpha$ e $0\beta$). Analogamente, se i qubit non sono fisicamente adiacenti su una QPU, anche questi gate richiedono SWAP. Dalla discussione precedente emerge che l'ansatz UCJ presenta alcune difficoltà per l'esecuzione su hardware, poiché richiede gate SWAP dovuti alle interazioni tra qubit non adiacenti. La variante locale dell'ansatz UCJ, LUCJ, affronta questa sfida rimuovendo alcuni gate $U_{nn}$ dall'operatore di Coulomb diagonale.

Nei blocchi della stessa specie elettronica, $J_{\alpha \alpha}$ e $J_{\beta \beta}$), nella versione LUCJ si mantengono solo i gate $U_{nn}$ compatibili con la connettività ai primi vicini e si rimuovono i gate tra qubit non adiacenti. La figura seguente mostra il blocco LUCJ dopo la rimozione dei gate non adiacenti. La versione LUCJ del blocco $J_{\alpha \beta}$, che agisce tra specie elettroniche diverse, può assumere forme diverse in base alla topologia del dispositivo.

Anche qui, la versione LUCJ elimina i gate non compatibili. La figura seguente mostra le varianti del blocco $J_{\alpha \beta}$ per diverse topologie di qubit: griglia, esagonale, heavy-hex e lineare.

• Quadrata: è possibile avere gate $U_{nn}$ tra tutti gli orbitali $\alpha$ e $\beta$ senza SWAP, quindi non è necessario rimuovere alcun gate $U_{nn}$.
• Heavy-hex: Le interazioni $\alpha$-$\beta$ sono mantenute ogni $4$ orbitali indicizzati (ad esempio il 0°, il 4° e l'8°) e sono mediate da qubit ancilla, ovvero servono qubit ancilla tra le catene lineari che rappresentano gli orbitali $\alpha$ e $\beta$. Questa configurazione richiede un numero limitato di SWAP.
• Esagonale: Ogni altro orbitale, come il 0°, il 2° e il 4° indicizzati, diventa vicino più prossimo quando $\alpha$ e $\beta$ sono disposti in due catene lineari adiacenti.
• Lineare: Solo un orbitale $\alpha$ e uno $\beta$ sono connessi, il che significa che il blocco $J_{\alpha \beta}$ avrà un solo gate. Sebbene la rimozione di gate dall'ansatz UCJ per costruire la versione LUCJ la renda più compatibile con l'hardware, l'ansatz perde parte dell'espressività. Pertanto, quando si usa l'ansatz LUCJ potrebbero essere necessarie più ripetizioni ($L$) dell'operatore UCJ modificato.

1.2 Inizializzazione dell'ansatz LUCJ

Il LUCJ è un ansatz parametrizzato e occorre inizializzare i parametri prima dell'esecuzione su hardware. Un modo per inizializzare l'ansatz consiste nell'usare le ampiezze t1 e t2 del metodo classico coupled cluster con singole e doppie eccitazioni (CCSD), dove le ampiezze t1 sono i coefficienti degli operatori di singola eccitazione e le ampiezze t2 sono quelli degli operatori di doppia eccitazione.

Tieni presente che, sebbene l'inizializzazione dell'ansatz LUCJ con le ampiezze t1 e t2 produca risultati discreti, 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

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 Costruire l'ansatz LUCJ con ffsim

Useremo 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 shell chiusa, useremo la variante bilanciata in spin dell'ansatz UCJ, UCJOpSpinBalanced.

Poiché l'hardware IBM ha una topologia heavy-hex, adotteremo il pattern zig-zag usato in [1] e descritto sopra per le interazioni tra qubit. In questo pattern, gli orbitali (qubit) con lo stesso spin sono connessi con una topologia a linea (cerchi rossi e blu). A causa della topologia heavy-hex, gli orbitali di spin diversi hanno connessioni ogni 4° orbitale, ovvero il 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'ansatz LUCJ con strati ripetuti può essere ottimizzato unendo alcuni blocchi adiacenti. Considera il caso n_reps=2: i due blocchi di rotazione orbitale al centro possono essere fusi in un unico blocco. Il pacchetto ffsim dispone di un pass manager chiamato ffsim.qiskit.PRE_INIT per ottimizzare il circuito unendo questi blocchi adiacenti.

2. Ottimizzare per l'hardware di destinazione

Per prima cosa, recuperiamo il backend di nostra scelta. Ottimizzeremo il circuito per il backend scelto 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")

Di seguito raccomandiamo i seguenti passi per ottimizzare l'ansatz e renderlo compatibile con l'hardware.

• Seleziona i qubit fisici (initial_layout) dall'hardware di destinazione che rispettano il pattern zig-zag (due catene lineari con qubit ancilla tra di esse) descritto sopra. Disporre i qubit in questo pattern porta a un circuito efficiente e compatibile con l'hardware, con meno gate.
• Genera un pass manager a stadi usando la funzione generate_preset_pass_manager di Qiskit con il backend e l'initial_layout scelti.
• Imposta lo stadio pre_init del tuo pass manager a stadi su ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT include pass del transpiler Qiskit che decompongono i gate in rotazioni orbitali e poi le uniscono, riducendo il numero di gate nel circuito finale.
• Esegui il pass manager sul tuo 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()}")

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 su hardware, siamo pronti a eseguirlo sull'hardware di destinazione e raccogliere campioni per la stima dell'energia dello stato fondamentale. Poiché abbiamo un solo circuito, useremo la modalità di esecuzione Job di Qiskit Runtime ed 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. Post-elaborare i risultati

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

Il campionamento dell'ansatz LUCJ nella base computazionale genera un pool di configurazioni rumorose $\tilde{\mathcal{\chi}}$, usate nella routine di post-elaborazione. Questa include un metodo chiamato recupero della configurazione (descritto in dettaglio più avanti) per correggere probabilisticamente le configurazioni con un numero errato di elettroni. Solo le configurazioni con il numero corretto di elettroni $\tilde{\mathcal{\chi}}_{R}$ vengono poi sottocampionate e distribuite in più batch in base alla frequenza di comparsa di ogni configurazione unica. Ogni batch di campioni definisce un sottospazio ($\mathcal{S^{(k)}}$). L'Hamiltoniano molecolare $H$ viene quindi proiettato sui sottospazi:

$$H_{\mathcal{S}^{(k)}} = P_{\mathcal{S}^{(k)}} H _{\mathcal{S}^{(k)}} \text{ with } P_{\mathcal{S}^{(k)}} = \sum_{x \in \mathcal{S}^{(k)}} \vert x \rangle \langle x \vert$$

Ogni Hamiltoniano proiettato $H_{\mathcal{S}^{(k)}}$ viene poi fornito a un eigensolver, che lo diagonalizza per calcolare autovalori e autovettori e ricostruire un autostato. In questa lezione, proiettiamo e diagonalizziamo l'Hamiltoniano usando il pacchetto qiskit-addon-sqd, che utilizza il metodo di Davidson di PySCF per la diagonalizzazione.

$$H_{\mathcal{S}^{(k)}} \vert \psi^{(k)} \rangle = E^{(k)} \vert \psi^{(k)} \rangle$$

Raccogliamo poi l'autovalore (energia) minimo dai batch e calcoliamo l'occupazione orbitale media $\text{n}$. L'informazione sull'occupazione media viene usata nel passo di recupero della configurazione per correggere probabilisticamente le configurazioni rumorose.

Successivamente, spieghiamo in dettaglio il ciclo di recupero della configurazione auto-consistente e mostriamo esempi di codice concreti per implementare i passi descritti sopra e stimare l'energia dello stato fondamentale dell'Hamiltoniano $N_2$.

4.1 Recupero della configurazione: panoramica

Ogni bit in una bitstring (determinante di Slater) rappresenta uno spin-orbitale. La metà destra di una bitstring rappresenta gli orbitali spin-up, e la metà sinistra gli orbitali spin-down. Un 1 significa che l'orbitale è occupato da un elettrone, e uno 0 che l'orbitale è vuoto. Conosciamo a priori il numero corretto di particelle (sia elettroni con spin-up che con spin-down). Supponiamo di avere un determinante $x$ con $N_x$ elettroni (ovvero $N_x$ valori 1 nella bitstring). Il numero corretto di particelle è $N$. Se $N_x \neq N$, sappiamo che la bitstring è stata corrotta dal rumore. La routine di configurazione auto-consistente tenta di correggere la bitstring ribaltando probabilisticamente $|N_x - N|$ bit, sfruttando le informazioni sull'occupazione orbitale media. L'occupazione orbitale media ($n$) ci dice quanto è probabile che un orbitale sia occupato da un elettrone. Se $N_x < N$, abbiamo meno elettroni del necessario e dobbiamo ribaltare alcuni 0 in 1, e viceversa.

La probabilità di ribaltamento può essere $|x[i] - avg\_occupancy[i]|$ per il $i$-esimo spin-orbitale. In [2], gli autori hanno usato una probabilità di ribaltamento pesata tramite una funzione ReLU modificata.

$$\begin{align} w(y) = \begin{cases} \delta \frac{y}{h} & \text{if } y \leq h\\ \nonumber \delta + (1 - \delta) \frac{y - h}{1 - h} & \text{if } y > h \end{cases} \end{align}$$

Qui $h$ definisce la posizione del "gomito" della funzione ReLU, e il parametro $\delta$ il valore della funzione ReLU al gomito. Per $\delta = 0$, $w$ diventa la vera funzione ReLU, e per $\delta > 0$ diventa la ReLU modificata. Nell'articolo, gli autori hanno usato $\delta = 0.01$ e $h =$ numero di particelle alfa (o beta) / numero di spin-orbitali 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 parte dalle configurazioni con solo il numero corretto di particelle in entrambe le specie di spin. Dopo la prima iterazione, abbiamo una stima dello stato fondamentale e, usando questa stima, possiamo costruire la prima approssimazione di $n$. Questa approssimazione di $n$ viene usata per recuperare le configurazioni, eseguire la successiva iterazione della stima dello stato fondamentale e raffinare auto-consistentemente la stima di $n$. Il processo si ripete finché non viene soddisfatto un criterio di arresto.

Considera il seguente esempio per $N = 2$ e $x = |1000\rangle$ ($N_x = 1$). Dobbiamo ribaltare uno degli 0 in 1 per correggere il numero di particelle, e le scelte sono 1100, 1010 e 1001. In base alla probabilità di ribaltamento, una delle scelte sarà selezionata come configurazione recuperata (ovvero la bitstring con il numero corretto di particelle).

Supponiamo che nella prima iterazione eseguiamo due batch e che gli stati fondamentali stimati da essi siano:

$$\begin{align}\nonumber \text{Batch0: } \vert \psi \rangle &= 0.8 \times \vert 1001 \rangle + 0.6 \times \vert 0110 \rangle \\ \nonumber \text{Batch1: } \vert \psi \rangle &= \frac{1}{\sqrt{3}} \left( \vert 1001 \rangle + \vert 0101 \rangle + \vert 0110 \rangle \right) \nonumber \end{align}$$

Usando gli stati della base computazionale e le loro ampiezze, possiamo calcolare la probabilità di occupazione degli elettroni (in breve, occupazioni) per spin-orbitale (qubit) (nota che probabilità = |ampiezza|$^2$). Nella tabella seguente vengono riportate le occupazioni per qubit di ogni bitstring che appare nello stato fondamentale stimato, e viene calcolata l'occupazione orbitale totale per un batch. Nota che, secondo la convenzione di ordinamento di Qiskit, il bit più a destra rappresenta il qubit-0 (Q0), e il bit più a sinistra rappresenta Q3.

Occupazione (Batch0):

	Q3	Q2	Q1	Q0
1001	0.64	0.0	0.0	0.64
0110	0.0	0.36	0.36	0.0
n (Batch0)	0.64	0.36	0.36	0.64

Occupazione (Batch1)

	Q3	Q2	Q1	Q0
1001	0.33	0.00	0.00	0.33
0101	0.0	0.33	0.00	0.33
0110	0.0	0.33	0.33	0.00
n (Batch1)	0.33	0.66	0.33	0.66

Occupazione (media tra i batch)

	Q3	Q2	Q1	Q0
n (Batch0)	0.64	0.36	0.36	0.64
n (Batch1)	0.33	0.66	0.33	0.66
n (media)	0.49	0.51	0.35	0.65

Usando l'occupazione orbitale media calcolata sopra, possiamo trovare le probabilità di ribaltamento per i diversi orbitali nella configurazione $x = \vert 1000 \rangle$. Poiché l'orbitale rappresentato da Q3 è già occupato e non deve essere ribaltato, impostiamo p(ribaltamento) = $0$. Per i rimanenti orbitali non occupati, la probabilità di ribaltamento è $\vert x[i] - \text{n}[i] \vert$ per ciascuno. Insieme a p(ribaltamento), calcoliamo anche il peso di probabilità associato al ribaltamento usando la funzione ReLU modificata descritta sopra.

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

	Q3	Q2	Q1	Q0
p(ribaltamento) ($\vert x[i] - \text{n}[i] \vert$)	0	0.51	0.35	0.65
w(p(ribaltamento))	0	0.03	0.007	0.31

Infine, usando le probabilità pesate sopra, possiamo ribaltare uno degli orbitali non occupati Q2, Q1 e Q0. In base ai valori calcolati, Q0 verrà ribaltato con la probabilità maggiore, e una possibile configurazione recuperata può essere $\vert \text{1001} \rangle$. Il processo completo di recupero della configurazione auto-consistente 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 include sia configurazioni con il numero corretto ($\widetilde{\chi}_{correct}$) che errato ($\widetilde{\chi}_{incorrect}$) di particelle in ciascun settore di spin.

1. Le configurazioni da ($\widetilde{\chi}_{correct}$) vengono campionate casualmente per creare batch $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ di vettori per la proiezione sul sottospazio. Il numero di batch e di campioni per ogni batch sono parametri definiti dall'utente. Più grande è il numero di campioni per batch, maggiore è la dimensione del sottospazio e più onerosa diventa la diagonalizzazione dal punto di vista computazionale. D'altro canto, un numero troppo piccolo di campioni potrebbe non includere i vettori di supporto dello stato fondamentale, portando a una stima errata.
2. Esegui l'eigensolver (ovvero la proiezione sul sottospazio e la diagonalizzazione) sui batch e ottieni gli autostati approssimati $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. Dagli autostati approssimati, costruisci la prima stima di $n$.

Iterazioni successive:

1. Usando $n$, correggi le configurazioni con il numero errato di particelle in $\widetilde{\chi}_{incorrect}$. Chiamiamo queste $\widetilde{\chi}_{correct\_new}$. Allora $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ forma il nuovo insieme di configurazioni con il numero corretto di particelle.
2. $\widetilde{\chi}_{R}$ viene campionato per creare i batch $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. L'eigensolver viene eseguito con i nuovi batch e genera nuove stime degli stati fondamentali $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. Dagli autostati approssimati, costruisci una stima raffinata di $n$.
5. Se il criterio di arresto non è soddisfatto, torna al passo 2.1.

4.2 Stima dello stato fondamentale

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

Ogni riga della matrice rappresenta una bitstring unica. Poiché i qubit sono indicizzati da destra in una bitstring 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), e gli orbitali beta nell'intervallo (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 controllabili dall'utente che sono importanti per questa tecnica:

• iterations: Numero di iterazioni del recupero della configurazione auto-consistente
• n_batches: Numero di batch di configurazioni usati dalle diverse chiamate all'eigensolver
• samples_per_batch: Numero di configurazioni uniche da includere in ogni batch
• max_davidson_cycles: Numero massimo di cicli di Davidson eseguiti da ogni 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 post-select 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_davidson=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)

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 poche iterazioni stimiamo l'energia dello stato fondamentale con un errore di circa 24 mH (la precisione chimica è tipicamente accettata come 1 kcal/mol $\approx$ 1.6 mH). Il secondo grafico mostra l'occupazione media di ogni orbitale spaziale dopo l'ultima iterazione. Si può vedere 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 stimata dello stato fondamentale sia discreta, non è entro il limite di precisione chimica ($\pm \approx 1.6$ mH). Questo divario può essere attribuito alla piccola dimensione del sottospazio usata sopra per la proiezione e la diagonalizzazione. Poiché abbiamo usato samples_per_batch=500, il sottospazio è spannato da al massimo $500$ vettori, che non coprono tutti i vettori di supporto dello stato fondamentale. Aumentare il parametro samples_per_batch dovrebbe migliorare la precisione a scapito di maggiori risorse di calcolo classico e 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()

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

Esercizio per il lettore

Aumenta progressivamente il parametro samples_per_batch (ad esempio, da $1000$ a $10000$ con un passo di $1000$; compatibilmente con la memoria del tuo computer) e confronta le energie stimate dello stato fondamentale.

Riferimenti

[1] M. Motta et al., "Bridging physical intuition and hardware efficiency for correlated electronic states: the local unitary cluster Jastrow ansatz for electronic structure" (2023). Chem. Sci., 2023, 14, 11213.

[2] J. Robledo-Moreno et al., "Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer" (2024). arXiv:quant-ph/2405.05068.