{
"cells": [
{
"cell_type": "markdown",
"id": "090b6884",
"metadata": {},
"source": [
"---\n",
"title: \"Implementazione SQD\"\n",
"description: \"La diagonalizzazione quantistica basata su campioni (SQD) viene implementata nel contesto della risoluzione dello stato fondamentale di una molecola di azoto.\"\n",
"---\n",
"\n",
"{/* cspell:ignore fontdict fontsize milli hcore ccsd Motta */}\n",
"\n",
"\n",
"\n",
"# SQD per la stima energetica di un hamiltoniano chimico\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ee7de1bc",
"metadata": {},
"source": [
"In questa lezione, applicheremo l'SQD per stimare l'energia di stato fondamentale di una molecola.\n",
"\n",
"In particolare, verranno affrontati i seguenti argomenti utilizzando l'approccio del pattern $4$ Qiskit:\n",
"\n",
"1. Passo 1: mappare il problema in circuiti e operatori quantistici\n",
" * Impostazione dell'hamiltoniana molecolare per $N_2$.\n",
" * Spiega il cluster unitario locale Jastrow (LUCJ) ispirato alla chimica e compatibile con l'hardware [\\[1\\]](#references)\n",
"2. Fase 2: Ottimizzazione per l'hardware di destinazione\n",
" * Ottimizzare il numero di gate e il layout dell'ansatz per l'esecuzione hardware\n",
"3. Fase 3: Esecuzione sull'hardware di destinazione\n",
" * Eseguire il circuito ottimizzato su una QPU reale per generare campioni del sottospazio.\n",
"4. Fase 4: Postelaborazione dei risultati\n",
" * Introdurre il ciclo di recupero della configurazione autoconsistente [\\[2\\]](#references)\n",
" * 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.\n",
" * Creare in modo probabilistico lotti di sottocampioni dalle stringhe di bit recuperate.\n",
" * Proiettare e diagonalizzare l'hamiltoniana molecolare su ogni sottospazio campionato.\n",
" * Salva l'energia minima del ground state trovata in tutti i batch e aggiorna l'occupazione avg degli orbitali.\n",
"\n",
"Nel corso della lezione verranno utilizzati diversi pacchetti software.\n",
"\n",
"* `PySCF` per definire la molecola e impostare l'hamiltoniana.\n",
"* `ffsim` per costruire l'ansatz LUCJ.\n",
"* `Qiskit` per il transpiling dell'ansatz per l'esecuzione hardware.\n",
"* `Qiskit IBM Runtime` per eseguire il circuito su una QPU e raccogliere campioni.\n",
"* `Qiskit addon SQD` recupero della configurazione e stima dell'energia dello stato di terra utilizzando la proiezione del sottospazio e la diagonalizzazione della matrice.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "719a9c0e-8c00-4fab-ba23-f2c8a7ebb573",
"metadata": {},
"source": [
"\n",
"\n",
"## 1. Mappare il problema su circuiti quantistici e operatori\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d02f97af",
"metadata": {},
"source": [
"\n",
"\n",
"### Hamiltoniano molecolare\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "5dac27ce-56b1-4f7f-ab83-ac153c004e82",
"metadata": {},
"source": [
"L'hamiltoniana molecolare assume la forma generica:\n",
"\n",
"$$\n",
"\\hat{H} = \\sum_{ \\substack{pr\\\\\\sigma} } h_{pr} \\, \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{r\\sigma}\n",
"+\n",
"\\sum_{ \\substack{prqs\\\\\\sigma\\tau} }\n",
"\\frac{(pr|qs)}{2} \\,\n",
"\\hat{a}^\\dagger_{p\\sigma}\n",
"\\hat{a}^\\dagger_{q\\tau}\n",
"\\hat{a}_{s\\tau}\n",
"\\hat{a}_{r\\sigma}\n",
"$$\n",
"\n",
"$\\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`.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cc987c08-7261-4c4c-a06b-609d7003efe9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"converged SCF energy = -108.835236570774\n",
"CASCI E = -109.046671778080 E(CI) = -32.8155692383188 S^2 = 0.0000000\n"
]
}
],
"source": [
"import warnings\n",
"import pyscf\n",
"import pyscf.cc\n",
"import pyscf.mcscf\n",
"\n",
"warnings.filterwarnings(\"ignore\")\n",
"\n",
"# Specify molecule properties\n",
"open_shell = False\n",
"spin_sq = 0\n",
"\n",
"# Build N2 molecule\n",
"mol = pyscf.gto.Mole()\n",
"mol.build(\n",
" atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]], # Two N atoms 1 angstrom apart\n",
" basis=\"6-31g\",\n",
" symmetry=\"Dooh\",\n",
")\n",
"\n",
"# Define active space\n",
"n_frozen = 2\n",
"active_space = range(n_frozen, mol.nao_nr())\n",
"\n",
"# Get molecular integrals\n",
"scf = pyscf.scf.RHF(mol).run()\n",
"num_orbitals = len(active_space)\n",
"n_electrons = int(sum(scf.mo_occ[active_space]))\n",
"num_elec_a = (n_electrons + mol.spin) // 2\n",
"num_elec_b = (n_electrons - mol.spin) // 2\n",
"cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))\n",
"mo = cas.sort_mo(active_space, base=0)\n",
"hcore, nuclear_repulsion_energy = cas.get_h1cas(mo) # hcore: one-body integrals\n",
"eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals) # eri: two-body integrals\n",
"\n",
"# Compute exact energy for comparison\n",
"exact_energy = cas.run().e_tot"
]
},
{
"cell_type": "markdown",
"id": "3aa4e0f0",
"metadata": {},
"source": [
"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$.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "f6a89c89-85e1-4c0d-abeb-8ab8b154a7ba",
"metadata": {},
"source": [
"\n",
"\n",
"### 1.1 Circuito quantistico per la generazione di campioni: l'approccio LUCJ\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "a464c865-1528-45c2-8ede-325583b15976",
"metadata": {},
"source": [
"In questa lezione, utilizzeremo l'ipotesi di cluster unitario locale accoppiato di Jastrow (LUCJ) \\ \\[[1\\]](#references) 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.\n",
"\n",
"L'ansatz di UCJ ha la forma seguente (per un prodotto di $L$ strati o ripetizioni dell'operatore UCJ)\n",
"\n",
"$$\n",
"|\\psi\\rangle = \\prod_{\\mu=1}^{L}{(e^{K^{\\mu}} \\times {e^{iJ^{\\mu}}} \\times {e^{-K^{\\mu}}})} |\\Phi_{0}\\rangle\n",
"$$\n",
"\n",
"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:\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "0f0bb614",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "4a32d88b",
"metadata": {},
"source": [
"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)}}$ ).\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "963e8386-39b6-40b7-9740-fffcf1573fe6",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "271b1825-bc78-4957-bbf8-21a714e45ed3",
"metadata": {},
"source": [
"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$ ).\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "2e8ac1d2-8f04-4591-921b-8ba0174e4ad0",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ebc9b48c-e77d-4af5-b8cf-460daeeadcdb",
"metadata": {},
"source": [
"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)}}$ ).\n",
"\n",
"Tutti i blocchi di $e^{iJ^{(\\mu)}}$ sono costituiti da porte numeriche $U_{nn}(\\phi)$ [\\[1\\]](#references). 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.\n",
"\n",
"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.\n",
"\n",
"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).\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d3e24a20-1c86-4aea-8300-13bd2ead4e8f",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ef3f1d36-e8ab-46f0-bddb-fa8fb884e531",
"metadata": {},
"source": [
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9afe0036-318c-43b9-9b2c-39a808bda82c",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "6219aa8e-8a63-4bfc-b52b-64507292471d",
"metadata": {},
"source": [
"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.\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "24f54a55-4a09-4508-997a-96306835c7e6",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "165d0a19-0366-4af8-8ea0-20df351f49f5",
"metadata": {},
"source": [
"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.\n",
"\n",
"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.\n",
"\n",
"* **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}$.\n",
"* **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.\n",
"* **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.\n",
"* **Lineare** : Sono collegati solo un orbitale $\\alpha$ e uno $\\beta$, il che significa che il blocco $J_{\\alpha \\beta}$ avrà un solo gate.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8ec1e433-bc00-42bc-bdbb-288c56c32f9d",
"metadata": {},
"source": [
"Schemi di \n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "3e070615-2b4a-4866-8c6c-d2d038447f45",
"metadata": {},
"source": [
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "04367dac",
"metadata": {},
"source": [
"\n",
"\n",
"### 1.2 Inizializzazione dell'approccio LUCJ\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "f5eae55f",
"metadata": {},
"source": [
"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.\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "1835e74e-3354-425f-8596-c574f03e7a6e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"E(CCSD) = -109.0398256929733 E_corr = -0.20458912219883\n"
]
}
],
"source": [
"# Get CCSD t2 amplitudes for initializing the ansatz\n",
"ccsd = pyscf.cc.CCSD(\n",
" scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]\n",
")\n",
"ccsd.run()\n",
"\n",
"t1 = ccsd.t1\n",
"t2 = ccsd.t2"
]
},
{
"cell_type": "markdown",
"id": "9f0842fa",
"metadata": {},
"source": [
"\n",
"\n",
"### 1.3 Costruzione dell'approccio LUCJ utilizzando `ffsim`\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "e5f7f7e8-30e3-40e8-9b85-bf057b35b766",
"metadata": {},
"source": [
"Utilizzeremo il pacchetto [ffsim](https://github.com/qiskit-community/ffsim/tree/main) 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](https://qiskit-community.github.io/ffsim/api/ffsim.html#ffsim.UCJOpSpinBalanced).\n",
"\n",
"Poiché l'hardware di IBM ha una topologia heavy-hex, adotteremo lo schema *a zig-zag* utilizzato in [\\[1\\]](#references) 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).\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "6799421a-4425-404a-9994-47215957054d",
"metadata": {},
"outputs": [],
"source": [
"import ffsim\n",
"from qiskit import QuantumCircuit, QuantumRegister\n",
"\n",
"n_reps = 2\n",
"alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]\n",
"alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]\n",
"\n",
"ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
" t2=t2,\n",
" t1=t1,\n",
" n_reps=n_reps,\n",
" interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),\n",
")\n",
"\n",
"nelec = (num_elec_a, num_elec_b)\n",
"\n",
"# create an empty quantum circuit\n",
"qubits = QuantumRegister(2 * num_orbitals, name=\"q\")\n",
"circuit = QuantumCircuit(qubits)\n",
"\n",
"# prepare Hartree-Fock state as the reference state and append it to the quantum circuit\n",
"circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)\n",
"\n",
"# apply the UCJ operator to the reference state\n",
"circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n",
"circuit.measure_all()\n",
"# circuit.decompose().draw(\"mpl\", scale=0.5, fold=-1)"
]
},
{
"cell_type": "markdown",
"id": "3b2de316",
"metadata": {},
"source": [
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "7cb99cd9",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "e9ad460d",
"metadata": {},
"source": [
"\n",
"\n",
"## 2. Ottimizzazione per l'hardware di destinazione\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "cecf2994",
"metadata": {},
"source": [
"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.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4c8dbbef-378e-4ae3-9290-bde58b72024d",
"metadata": {},
"outputs": [],
"source": [
"from qiskit_ibm_runtime import QiskitRuntimeService\n",
"\n",
"service = QiskitRuntimeService()\n",
"# Use the least-busy backend or specify a quantum computer using the syntax commented out below.\n",
"backend = service.least_busy(operational=True, simulator=False)\n",
"# backend = service.backend(\"ibm_brisbane\")"
]
},
{
"cell_type": "markdown",
"id": "1b8ba388-6904-4725-ad97-2a31c0adaa07",
"metadata": {},
"source": [
"Per ottimizzare il surrogato e renderlo compatibile con l'hardware, si consiglia di eseguire le seguenti operazioni.\n",
"\n",
"* 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.\n",
"* Generare uno staged pass manager utilizzando la funzione [`generate_preset_pass_manager`](/docs/api/qiskit/qiskit.transpiler.generate_preset_pass_manager) di Qiskit con la scelta di `backend` e `initial_layout`.\n",
"* 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.\n",
"* Eseguire il pass manager sul circuito.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "ecd2675c-a302-43fb-80f4-d658d56360d5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Gate counts (w/o pre-init passes): OrderedDict({'rz': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})\n",
"Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})\n"
]
}
],
"source": [
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
"\n",
"spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]\n",
"spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]\n",
"\n",
"initial_layout = spin_a_layout + spin_b_layout\n",
"\n",
"pass_manager = generate_preset_pass_manager(\n",
" optimization_level=3, backend=backend, initial_layout=initial_layout\n",
")\n",
"\n",
"# without PRE_INIT passes\n",
"isa_circuit = pass_manager.run(circuit)\n",
"print(f\"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}\")\n",
"\n",
"# with PRE_INIT passes\n",
"# We will use the circuit generated by this pass manager for hardware execution\n",
"pass_manager.pre_init = ffsim.qiskit.PRE_INIT\n",
"isa_circuit = pass_manager.run(circuit)\n",
"print(f\"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}\")"
]
},
{
"cell_type": "markdown",
"id": "72a780f9",
"metadata": {},
"source": [
"\n",
"\n",
"## 3. Eseguire sull'hardware di destinazione\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "4cd9164c-697a-4aa6-b71f-86720e3d5b66",
"metadata": {},
"source": [
"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-runtime`[Modalità di esecuzione del processo](/docs/guides/execution-modes) eseguiremo il nostro circuito.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "9f542448-5767-4e12-85c8-dd8584545dbb",
"metadata": {},
"outputs": [],
"source": [
"from qiskit_ibm_runtime import SamplerV2 as Sampler\n",
"\n",
"sampler = Sampler(mode=backend)\n",
"sampler.options.dynamical_decoupling.enable = True\n",
"\n",
"job = sampler.run([isa_circuit], shots=10_000) # Takes approximately 5sec of QPU time"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8c5eb3e9-2a5a-423b-ac18-7bba7f69f18f",
"metadata": {},
"outputs": [],
"source": [
"# Run cell after IQX job completion\n",
"primitive_result = job.result()\n",
"pub_result = primitive_result[0]\n",
"counts = pub_result.data.meas.get_counts()"
]
},
{
"cell_type": "markdown",
"id": "e3ebf77a-8efd-43d6-bd19-e26423921c6f",
"metadata": {},
"source": [
"\n",
"\n",
"## 4. Risultati post-elaborazione\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "1e7f0ecd",
"metadata": {},
"source": [
"La parte di post-elaborazione del flusso di lavoro SQD può essere riassunta con il seguente diagramma.\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8c51a527",
"metadata": {},
"source": [
"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:\n",
"\n",
"$$\n",
"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\n",
"$$\n",
"\n",
"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.\n",
"\n",
"$$\n",
"H_{\\mathcal{S}^{(k)}} \\vert \\psi^{(k)} \\rangle = E^{(k)} \\vert \\psi^{(k)} \\rangle\n",
"$$\n",
"\n",
"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.\n",
"\n",
"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$.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "39997b5d-8bd3-4bc8-9ade-11fa5cfb34f9",
"metadata": {},
"source": [
"\n",
"\n",
"### 4.1 Ripristino della configurazione: panoramica\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "b05881ac-80ce-47a6-ac28-c1237f85bc0a",
"metadata": {},
"source": [
"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.\n",
"\n",
"La probabilità di capovolgimento può essere $|x[i] - avg\\_occupancy[i]|$ per `i`-esimo orbitale di spin. In [\\[2\\]](#references), gli autori hanno utilizzato una probabilità ponderata di ribaltamento utilizzando la funzione ReLU modificata.\n",
"\n",
"$$\n",
"\\begin{align}\n",
" w(y) = \\begin{cases}\n",
"\n",
" \\delta \\frac{y}{h} & \\text{if } y \\leq h\\\\ \\nonumber\n",
"\n",
" \\delta + (1 - \\delta) \\frac{y - h}{1 - h} & \\text{if } y > h\n",
"\n",
"\\end{cases}\n",
"\\end{align}\n",
"$$\n",
"\n",
"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).\n",
"\n",
"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.\n",
"\n",
"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).\n",
"\n",
"Supponiamo che nella prima iterazione vengano eseguiti due lotti e che gli stati a terra stimati da questi siano:\n",
"\n",
"$$\n",
"\\begin{align}\\nonumber\n",
" \\text{Batch0: } \\vert \\psi \\rangle &= 0.8 \\times \\vert 1001 \\rangle + 0.6 \\times \\vert 0110 \\rangle \\\\ \\nonumber\n",
" \\text{Batch1: } \\vert \\psi \\rangle &= \\frac{1}{\\sqrt{3}} \\left( \\vert 1001 \\rangle + \\vert 0101 \\rangle + \\vert 0110 \\rangle \\right) \\nonumber\n",
"\\end{align}\n",
"$$\n",
"\n",
"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.\n",
"\n",
"Occupazione ( Batch0 ):\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :--------------: | :------: | :------: | :------: | :------: |\n",
"| 1001 | 0.64 | 0.0 | 0.0 | 0.64 |\n",
"| 0110 | 0.0 | 0.36 | 0.36 | 0.0 |\n",
"| **n** *(Batch0)* | **0.64** | **0.36** | **0.36** | **0.64** |\n",
"\n",
"Occupazione ( Batch1 )\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :--------------: | :------: | :------: | :------: | :------: |\n",
"| 1001 | 0.33 | 0.00 | 0.00 | 0.33 |\n",
"| 0101 | 0.0 | 0.33 | 0.00 | 0.33 |\n",
"| 0110 | 0.0 | 0.33 | 0.33 | 0.00 |\n",
"| **n** *(Batch1)* | **0.33** | **0.66** | **0.33** | **0.66** |\n",
"\n",
"Occupazione (media dei lotti)\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :--------------: | :------: | :------: | :------: | :------: |\n",
"| **n** *(Batch0)* | 0.64 | 0.36 | 0.36 | 0.64 |\n",
"| **n** *(Batch1)* | 0.33 | 0.66 | 0.33 | 0.66 |\n",
"| **n** *(media)* | **0.49** | **0.51** | **0.35** | **0.65** |\n",
"\n",
"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.\n",
"\n",
"Probabilità di capovolgimento ( $x = \\vert 1000 \\rangle$, $\\delta = 0.01$, $h = N/M = 2/4 = 0.50$ )\n",
"\n",
"| | Q3 | Q2 | Q1 | Q0 |\n",
"| :------------------------------------------: | :-: | :--: | :---: | :--: |\n",
"| p(flip) ( $\\vert x[i] - \\text{n}[i] \\vert$ ) | 0 | 0.51 | 0.35 | 0.65 |\n",
"| w(p(flip)) | 0 | 0.03 | 0.007 | 0.31 |\n",
"\n",
"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$.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "c1940235",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "5c614290",
"metadata": {},
"source": [
"Il processo completo di ripristino della configurazione autoconsistente può essere riassunto come segue:\n",
"\n",
"**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.\n",
"\n",
"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.\n",
"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$.\n",
"3. A partire dagli autostati approssimati si costruisce la prima ipotesi per $n$.\n",
"\n",
"**Iterazioni successive:**\n",
"\n",
"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.\n",
"2. $\\widetilde{\\chi}_{R}$ viene campionata per creare lotti $\\mathcal{S}^{(1)}, \\cdots, \\mathcal{S}^{(K)}$.\n",
"3. Il solutore di autostati viene eseguito con nuovi lotti e genera nuove stime degli stati fondamentali $|\\psi^{(1)}\\rangle, \\cdots, |\\psi^{(K)}\\rangle$.\n",
"4. Dagli autostati approssimati si costruiscono ipotesi raffinate per $n$.\n",
"5. Se il criterio di arresto non è soddisfatto, si torna al passo `2.1`.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "5eea548a",
"metadata": {},
"source": [
"\n",
"\n",
"### 4.2 Stima dello stato fondamentale\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "31dfe00e",
"metadata": {},
"source": [
"Per prima cosa, trasformeremo i conteggi in una matrice di bitstring e in un array di probabilità per la post-elaborazione.\n",
"\n",
"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.\n",
"\n",
"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).\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "71550274",
"metadata": {},
"outputs": [],
"source": [
"from qiskit_addon_sqd.counts import counts_to_arrays\n",
"\n",
"# Convert counts into bitstring and probability arrays\n",
"bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)"
]
},
{
"cell_type": "markdown",
"id": "acbdbdc5",
"metadata": {},
"source": [
"Ci sono alcune opzioni controllate dall'utente che sono importanti per questa tecnica:\n",
"\n",
"* `iterations`: Numero di iterazioni di recupero della configurazione autoconsistente\n",
"* `n_batches`: Numero di batch di configurazioni utilizzati dalle diverse chiamate al risolutore di autovalori\n",
"* `samples_per_batch`: Numero di configurazioni uniche da includere in ogni lotto\n",
"* `max_davidson_cycles`: Numero massimo di cicli Davidson eseguiti da ciascun eigensolver\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a2d22e0b-8f51-42a0-858f-ad0297cb0bae",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting configuration recovery iteration 0\n",
" Batch 0 subspace dimension: 21609\n",
" Batch 1 subspace dimension: 21609\n",
" Batch 2 subspace dimension: 21609\n",
" Batch 3 subspace dimension: 21609\n",
" Batch 4 subspace dimension: 21609\n",
"Starting configuration recovery iteration 1\n",
" Batch 0 subspace dimension: 609961\n",
" Batch 1 subspace dimension: 616225\n",
" Batch 2 subspace dimension: 627264\n",
" Batch 3 subspace dimension: 633616\n",
" Batch 4 subspace dimension: 624100\n",
"Starting configuration recovery iteration 2\n",
" Batch 0 subspace dimension: 564001\n",
" Batch 1 subspace dimension: 605284\n",
" Batch 2 subspace dimension: 582169\n",
" Batch 3 subspace dimension: 559504\n",
" Batch 4 subspace dimension: 591361\n",
"Starting configuration recovery iteration 3\n",
" Batch 0 subspace dimension: 550564\n",
" Batch 1 subspace dimension: 549081\n",
" Batch 2 subspace dimension: 531441\n",
" Batch 3 subspace dimension: 527076\n",
" Batch 4 subspace dimension: 531441\n",
"Starting configuration recovery iteration 4\n",
" Batch 0 subspace dimension: 544644\n",
" Batch 1 subspace dimension: 580644\n",
" Batch 2 subspace dimension: 527076\n",
" Batch 3 subspace dimension: 531441\n",
" Batch 4 subspace dimension: 537289\n"
]
}
],
"source": [
"import numpy as np\n",
"from qiskit_addon_sqd.configuration_recovery import recover_configurations\n",
"from qiskit_addon_sqd.fermion import (\n",
" bitstring_matrix_to_ci_strs,\n",
" solve_fermion,\n",
")\n",
"from qiskit_addon_sqd.subsampling import postselect_and_subsample\n",
"\n",
"rng = np.random.default_rng(24)\n",
"# SQD options\n",
"iterations = 5\n",
"\n",
"# Eigenstate solver options\n",
"n_batches = 5\n",
"samples_per_batch = 500\n",
"max_davidson_cycles = 300\n",
"\n",
"# Self-consistent configuration recovery loop\n",
"e_hist = np.zeros((iterations, n_batches)) # energy history\n",
"s_hist = np.zeros((iterations, n_batches)) # spin history\n",
"occupancy_hist = []\n",
"avg_occupancy = None\n",
"for i in range(iterations):\n",
" print(f\"Starting configuration recovery iteration {i}\")\n",
" # On the first iteration, we have no orbital occupancy information from the\n",
" # solver, so we begin with the full set of noisy configurations.\n",
" if avg_occupancy is None:\n",
" bs_mat_tmp = bitstring_matrix_full\n",
" probs_arr_tmp = probs_arr_full\n",
"\n",
" # If we have average orbital occupancy information, we use it to refine\n",
" # the full set of noisy configurations.\n",
" else:\n",
" bs_mat_tmp, probs_arr_tmp = recover_configurations(\n",
" bitstring_matrix_full,\n",
" probs_arr_full,\n",
" avg_occupancy,\n",
" num_elec_a,\n",
" num_elec_b,\n",
" rand_seed=rng,\n",
" )\n",
"\n",
" # Create batches of subsamples. We postselect here to remove configurations\n",
" # with incorrect hamming weight during iteration 0, since no config recovery was performed.\n",
" batches = postselect_and_subsample(\n",
" bs_mat_tmp,\n",
" probs_arr_tmp,\n",
" hamming_right=num_elec_a,\n",
" hamming_left=num_elec_b,\n",
" samples_per_batch=samples_per_batch,\n",
" num_batches=n_batches,\n",
" rand_seed=rng,\n",
" )\n",
"\n",
" # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.\n",
" e_tmp = np.zeros(n_batches)\n",
" s_tmp = np.zeros(n_batches)\n",
" occs_tmp = []\n",
" coeffs = []\n",
" for j in range(n_batches):\n",
" strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])\n",
" print(f\" Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}\")\n",
" energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(\n",
" batches[j],\n",
" hcore,\n",
" eri,\n",
" open_shell=open_shell,\n",
" spin_sq=spin_sq,\n",
" max_cycle=max_davidson_cycles,\n",
" )\n",
" energy_sci += nuclear_repulsion_energy\n",
" e_tmp[j] = energy_sci\n",
" s_tmp[j] = spin\n",
" occs_tmp.append(avg_occs)\n",
" coeffs.append(coeffs_sci)\n",
"\n",
" # Combine batch results\n",
" avg_occupancy = tuple(np.mean(occs_tmp, axis=0))\n",
"\n",
" # Track optimization history\n",
" e_hist[i, :] = e_tmp\n",
" s_hist[i, :] = s_tmp\n",
" occupancy_hist.append(avg_occupancy)"
]
},
{
"cell_type": "markdown",
"id": "e4f6ec2c-032d-4e18-8ea4-cf867cba5054",
"metadata": {},
"source": [
"\n",
"\n",
"### 4.3 Discussione dei risultati\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "592eabb7-18f3-4622-8710-bfc43ad6cdec",
"metadata": {},
"source": [
"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.\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "61959e69-a182-4636-abcb-a32349fc9076",
"metadata": {},
"outputs": [],
"source": [
"# Data for energies plot\n",
"x1 = range(iterations)\n",
"min_e = [np.min(e) for e in e_hist]\n",
"e_diff = [abs(e - exact_energy) for e in min_e]\n",
"yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]\n",
"\n",
"# Chemical accuracy (+/- 1 milli-Hartree)\n",
"chem_accuracy = 0.001\n",
"\n",
"# Data for avg spatial orbital occupancy\n",
"y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]\n",
"x2 = range(len(y2))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "8cd90034-6ef3-41bd-a847-c115cade82f7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact energy: -109.04667 Ha\n",
"SQD energy: -109.02234 Ha\n",
"Absolute error: 0.02434 Ha\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n",
"\n",
"# Plot energies\n",
"axs[0].plot(x1, e_diff, label=\"energy error\", marker=\"o\")\n",
"axs[0].set_xticks(x1)\n",
"axs[0].set_xticklabels(x1)\n",
"axs[0].set_yticks(yt1)\n",
"axs[0].set_yticklabels(yt1)\n",
"axs[0].set_yscale(\"log\")\n",
"axs[0].set_ylim(1e-6)\n",
"axs[0].axhline(\n",
" y=chem_accuracy, color=\"#BF5700\", linestyle=\"--\", label=\"chemical accuracy\"\n",
")\n",
"axs[0].set_title(\"Approximated Ground State Energy Error vs SQD Iterations\")\n",
"axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n",
"axs[0].set_ylabel(\"Energy Error (Ha)\", fontdict={\"fontsize\": 12})\n",
"axs[0].legend()\n",
"\n",
"# Plot orbital occupancy\n",
"axs[1].bar(x2, y2, width=0.8)\n",
"axs[1].set_xticks(x2)\n",
"axs[1].set_xticklabels(x2)\n",
"axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n",
"axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n",
"axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n",
"\n",
"print(f\"Exact energy: {exact_energy:.5f} Ha\")\n",
"print(f\"SQD energy: {min_e[-1]:.5f} Ha\")\n",
"print(f\"Absolute error: {e_diff[-1]:.5f} Ha\")\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "e1530472",
"metadata": {},
"source": [
"\n",
"\n",
"#### Esercizio per il lettore\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d2b2241f",
"metadata": {},
"source": [
"\n",
"\n",
"## Riferimenti\n",
"\n",
"\\[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](https://pubs.rsc.org/en/content/articlehtml/2023/sc/d3sc02516k).\n",
"\n",
"\\[2] J. Robledo-Moreno e altri, \"Chimica oltre le soluzioni esatte su un supercomputer incentrato sui quanti\" (2024). [arXiv:quant-ph/2405.05068](https://arxiv.org/abs/2405.05068).\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"id": "a1b8767d",
"source": "© IBM Corp., 2017-2026"
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}