SQD pro odhad energie chemického Hamiltoniánu

V této lekci použijeme SQD k odhadu energie základního stavu molekuly.

Konkrétně probereme následující témata pomocí $4$krokového přístupu Qiskit pattern:

1. Krok 1: Mapování problému na kvantové obvody a operátory
   • Nastavení molekulárního Hamiltoniánu pro $N_2$.
   • Vysvětlení chemicky inspirovaného a hardwarově vhodného lokálního unitárního klastrového Jastrowa (LUCJ) [1]
2. Krok 2: Optimalizace pro cílový hardware
   • Optimalizace počtu hradel a rozložení ansatzu pro hardwarové spuštění
3. Krok 3: Spuštění na cílovém hardwaru
   • Spuštění optimalizovaného obvodu na skutečném QPU pro generování vzorků z podprostoru.
4. Krok 4: Postprocessing výsledků
   • Představení samokonzistentní smyčky pro obnovu konfigurací [2]
     • Postprocessing úplné sady vzorků bitstringů s využitím předchozí znalosti počtu částic a průměrné orbitální obsazenosti vypočítané v poslední iteraci.
     • Pravděpodobnostní tvorba dávek podvzorků z obnovených bitstringů.
     • Projekce a diagonalizace molekulárního Hamiltoniánu nad každým navzorkovaným podprostorem.
     • Uložení minimální nalezené energie základního stavu napříč všemi dávkami a aktualizace průměrné orbitální obsazenosti.

V průběhu lekce využijeme několik softwarových balíčků.

• PySCF pro definici molekuly a nastavení Hamiltoniánu.
• balíček ffsim pro konstrukci LUCJ ansatzu.
• Qiskit pro transpilaci ansatzu pro hardwarové spuštění.
• Qiskit IBM Runtime pro spuštění obvodu na QPU a sběr vzorků.
• Qiskit addon SQD pro obnovu konfigurací a odhad energie základního stavu pomocí projekce na podprostor a maticové diagonalizace.

1. Mapování problému na kvantové obvody a operátory

Molekulární hamiltonián

Molekulární Hamiltonián má obecný tvar:

$$\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}$ jsou fermionové kreační/anihilační operátory spojené s $p$-tým prvkem bázové sady a spinem $\sigma$. $h_{pr}$ a $(pr|qs)$ jsou jedno- a dvoučásticové elektronové integrály. Pomocí pySCF definujeme molekulu a spočítáme jedno- a dvoučásticové integrály Hamiltoniánu pro bázovou sadu 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

V této lekci použijeme Jordan-Wignerovu (JW) transformaci k mapování fermionové vlnové funkce na qubitovou vlnovou funkci, aby ji bylo možné připravit pomocí kvantového obvodu. JW transformace mapuje Fockův prostor fermionů v M prostorových orbitalech na Hilbertův prostor 2M qubitů, to znamená, že prostorový orbital se rozdělí na dva spinorbitaly, jeden spojený s elektronem se spinem nahoru ($\alpha$) a druhý se spinem dolů ($\beta$). Spinorbital může být obsazený nebo neobsazený. Obvykle, když mluvíme o počtu orbitalů, máme na mysli počet prostorových orbitalů. Počet spinorbitalů bude dvojnásobný. V kvantových obvodech budeme každý spinorbital reprezentovat jedním qubitem. Sada qubitů bude tedy reprezentovat orbitaly spin-up neboli $\alpha$-orbitaly a další sada bude reprezentovat orbitaly spin-down neboli $\beta$-orbitaly. Například molekula $N_2$ pro bázovou sadu 6-31g má $16$ prostorových orbitalů (to jest $16$ $\alpha$ + $16$ $\beta$ = $32$ spinorbitalů). Budeme tedy potřebovat $32$qubitový kvantový obvod (později můžeme potřebovat další pomocné qubity, jak bude popsáno). Qubity se měří ve výpočetní bázi, čímž se generují bitstringy, které reprezentují elektronické konfigurace neboli (Slaterovy) determinanty. V této lekci budeme termíny bitstringy, konfigurace a determinanty používat zaměnitelně. Bitstringy nám říkají obsazenost elektronů ve spinorbitalech: $1$ na pozici bitu znamená, že odpovídající spinorbital je obsazený, zatímco $0$ znamená, že spinorbital je prázdný. Protože problémy elektronické struktury zachovávají částice, musí být obsazen pouze pevný počet spinorbitalů. Molekula $N_2$ má $5$ elektronů se spinem nahoru ($\alpha$) a $5$ elektronů se spinem dolů ($\beta$). Každý bitstring reprezentující $\alpha$ a $\beta$ orbitaly tedy musí mít pro molekulu $N_2$ pět $1\text{}$ v každém.

1.1 Kvantový obvod pro generování vzorků: ansatz LUCJ

V této lekci použijeme ansatz lokálního unitárního coupled cluster Jastrow (LUCJ) \[1\] pro přípravu kvantového stavu a následné vzorkování. Nejprve vysvětlíme jednotlivé stavební bloky plného ansatzu UCJ a aproximace provedené v jeho lokální verzi. Následně pomocí balíčku ffsim sestavíme ansatz LUCJ a optimalizujeme jej pomocí transpileru Qiskit pro hardwarové spuštění.

Ansatz UCJ má následující tvar (pro součin $L$ vrstev nebo opakování operátoru UCJ).

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

kde $\vert \Phi_{0} \rangle$ je referenční stav, typicky brán jako Hartree-Fockův (HF) stav. Protože je Hartree-Fockův stav definován tak, že má obsazené orbitaly s nejnižšími čísly, příprava HF stavu bude zahrnovat aplikaci X gates pro nastavení qubits odpovídajících obsazeným orbitalům na jedničku. Například blok přípravy HF stavu pro 4 prostorové orbitaly a 2 spiny nahoru a 2 spiny dolů může vypadat následovně: Jedno opakování operátoru UCJ ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ se skládá z diagonální Coulombovy evoluce ($e^{iJ^{(\mu)}}$) obklopené orbitálními rotacemi ($e^{K^{(\mu)}}$ a $e^{-K^{(\mu)}}$). Bloky orbitální rotace pracují na jednom druhu spinu ($\alpha$ (spin nahoru)/$\beta$ (spin dolů)). Pro každý druh elektronu se orbitální rotace skládá z vrstvy jednoqubitových hradel $R_{z}$ následovaných posloupností dvouqubitových Givensových rotačních hradel ($XX + YY$).

Dvouqubitová hradla působí na sousední spin-orbitaly (nejbližší sousední qubity), a tedy jsou implementovatelná na IBM® QPU bez nutnosti hradel SWAP. $e^{iJ^{(\mu)}}$, také známý jako diagonální Coulombův operátor, se skládá ze tří bloků. Dva z nich pracují na stejných spinových sektorech ($e^{iJ_{\alpha \alpha}^{(\mu)}}$ a $e^{iJ_{\beta \beta}^{(\mu)}}$) a jeden pracuje mezi dvěma spinovými sektory ($e^{iJ_{\alpha \beta}^{(\mu)}}$).

Všechny bloky v $e^{iJ^{(\mu)}}$ se skládají z hradel number-number $U_{nn}(\phi)$ [1]. Hradlo $U_{nn}(\phi)$ lze dále rozložit na hradlo $R_{ZZ}(\frac{\phi}{2})$ následované dvěma jednoqubitovými $Rz(-\frac{\phi}{2})$ hradly působícími na dva samostatné qubity.

Stejnospinové složky ($J_{\alpha \alpha}$ a $J_{\beta \beta}$) mají hradla $U_{nn}$ mezi všemi možnými dvojicemi qubitů. Protože však supravodivé QPU mají omezenou konektivitu, je nutné qubity prohazovat, aby bylo možné realizovat hradla mezi nesousedními qubity.

Zvažte například následující blok $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (nebo $e^{iJ_{\beta \beta}^{(\mu)}}$) pro $N = 4$ prostorové orbitaly. Pro lineární konektivitu qubitů nejsou poslední tři hradla přímo implementovatelná, protože pracují mezi nesousedními qubity (například Q0 a Q2 nejsou přímo propojeny). Proto potřebujeme hradla SWAP, abychom je udělali sousedními (následující obrázek ukazuje příklad s $3$ hradly SWAP). Dále $J_{\alpha \beta}$ implementuje hradla mezi stejně indexovanými orbitaly z různých spinových sektorů (například mezi $0\alpha$ a $0\beta$). Podobně, pokud qubity nejsou fyzicky sousední na QPU, budou tato hradla také vyžadovat hradla SWAP. Z výše uvedené diskuse vyplývá, že ansatz UCJ čelí určitým překážkám při HW spuštění, protože potřebuje hradla SWAP kvůli interakcím mezi nesousedními qubity. Lokální varianta ansatzu UCJ, LUCJ, řeší tuto výzvu odstraněním některých $U_{nn}$ z diagonálního Coulombova operátoru.

V blocích stejných druhů elektronů, $J_{\alpha \alpha}$ a $J_{\beta \beta}$, ponecháváme pouze hradla $U_{nn}$ kompatibilní s konektivitou nejbližších sousedů a ve verzi LUCJ odstraňujeme hradla mezi nesousedními qubity. Následující obrázek ukazuje blok LUCJ po odstranění nesousedních hradel. Dále verze LUCJ bloku $J_{\alpha \beta}$, který pracuje mezi různými druhy elektronů, může mít různý tvar v závislosti na topologii zařízení.

I zde se verze LUCJ zbavuje nekompatibilních hradel. Obrázek níže ukazuje varianty bloku $J_{\alpha \beta}$ pro různé topologie qubitů včetně mřížkové, hexagonální, heavy-hex a lineární.

• Čtvercová: můžeme mít hradla $U_{nn}$ mezi všemi $\alpha$ a $\beta$ orbitaly bez jakýchkoli hradel SWAP, a proto nepotřebujeme odstranit žádná hradla $U_{nn}$.
• Heavy-hex: Interakce $\alpha$-$\beta$ jsou ponechány mezi každým $4$. indexovaným (například 0., 4. a 8.) spin-orbitalem a jsou zprostředkovány ancillou, to znamená, že potřebujeme ancilla qubity mezi lineárními řetězy reprezentujícími $\alpha$ a $\beta$ orbitaly. Toto uspořádání potřebuje omezený počet hradel SWAP.
• Hexagonální: Každý druhý orbital, například 0., 2. a 4. indexovaný orbital, se stává nejbližším sousedem, když jsou $\alpha$ a $\beta$ rozloženy do dvou sousedních lineárních řetězů.
• Lineární: Propojený je pouze jeden $\alpha$ a jeden $\beta$ orbital, což znamená, že blok $J_{\alpha \beta}$ bude mít pouze jedno hradlo. Odstraněním hradel z ansatzu UCJ pro sestavení verze LUCJ se sice stává HW kompatibilnější, ale ansatz ztrácí určitou expresivitu. Proto při použití ansatzu LUCJ může být potřeba více opakování ($L$) modifikovaného operátoru UCJ.

1.2 Inicializace ansatzu LUCJ

LUCJ je parametrizovaný ansatz a před hardwarovým spuštěním musíme parametry inicializovat. Jedním ze způsobů inicializace ansatzu je použití amplitud t1 a t2 z klasické metody coupled cluster singles and doubles (CCSD), kde amplitudy t1 jsou koeficienty operátorů jednoduchých excitací a amplitudy t2 jsou pro operátory dvojitých excitací.

Všimni si, že i když inicializace ansatzu LUCJ amplitudami t1 a t2 generuje slušné výsledky, parametry ansatzu mohou vyžadovat další optimalizaci.

# 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 Konstrukce ansatzu LUCJ pomocí ffsim

K vytvoření a inicializaci ansatzu s amplitudami t1 a t2 vypočítanými výše použijeme balíček ffsim. Protože naše molekula má Hartree-Fockův stav s uzavřenou slupkou, použijeme spinově vyváženou variantu ansatzu UCJ, UCJOpSpinBalanced.

Protože hardware IBM má topologii heavy-hex, přijmeme cik-cak vzor použitý v [1] a vysvětlený výše pro interakce qubits. V tomto vzoru jsou orbitaly (qubits) se stejným spinem propojeny lineární topologií (červené a modré kruhy). Kvůli topologii heavy-hex mají orbitaly pro různé spiny spojení mezi každým 4. orbitalem, to jest 0., 4., 8. a tak dále (fialové kruhy).

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)

Ansatz LUCJ s opakovanými vrstvami lze optimalizovat sloučením některých sousedních bloků. Zvažme případ pro n_reps=2. Dva bloky orbitální rotace uprostřed lze sloučit do jednoho bloku orbitální rotace. Balíček ffsim má pass manager s názvem ffsim.qiskit.PRE_INIT k optimalizaci obvodu sloučením takových sousedních bloků.

2. Optimalizace pro cílový hardware

Nejprve načteme backend dle našeho výběru. Pro tento backend optimalizujeme náš obvod a poté na stejném backendu spustíme optimalizovaný obvod, abychom vygenerovali vzorky pro podprostor.

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")

Dále doporučujeme následující kroky k optimalizaci ansatzu a jeho přizpůsobení hardwaru.

• Vyber fyzické qubity (initial_layout) z cílového hardwaru, které odpovídají cik-cak vzoru (dva lineární řetězce s pomocným qubitem mezi nimi) popsanému výše. Rozmístění qubitů v tomto vzoru vede k efektivnímu hardwarově kompatibilnímu obvodu s menším počtem hradel.
• Vygeneruj fázovaný pass manager pomocí funkce generate_preset_pass_manager z Qiskitu s tvým výběrem backend a initial_layout.
• Nastav fázi pre_init tvého fázovaného pass manageru na ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT obsahuje průchody transpileru Qiskit, které rozkládají hradla na orbitální rotace a poté tyto orbitální rotace slučují, což vede k menšímu počtu hradel ve výsledném obvodu.
• Spusť pass manager na svém obvodu.

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. Spuštění na cílovém hardware

Po optimalizaci obvodu pro spuštění na hardwaru jsme připraveni jej spustit na cílovém hardwaru a sesbírat vzorky pro odhad energie základního stavu. Protože máme pouze jeden obvod, použijeme režim Job prostředí Qiskit Runtime a náš obvod spustíme.

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. Zpracování výsledků

Část SQD workflow věnovanou zpracování výsledků lze shrnout následujícím diagramem.

Vzorkování LUCJ ansatzu ve výpočetní bázi generuje zásobu zašuměných konfigurací $\tilde{\mathcal{\chi}}$, které jsou použity v rutině zpracování výsledků. Ta zahrnuje metodu zvanou (detaily probereme později) obnova konfigurací k pravděpodobnostnímu opravení konfigurací s nesprávným počtem elektronů. Konfigurace pouze se správným počtem elektronů $\tilde{\mathcal{\chi}}_{R}$ jsou poté podvzorkovány a rozděleny do více dávek na základě frekvence výskytu každé unikátní konfigurace. Každá dávka vzorků definuje podprostor ($\mathcal{S^{(k)}}$). Dále je molekulární Hamiltonián $H$ promítnut na podprostory:

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

Každý promítnutý Hamiltonián $H_{\mathcal{S}^{(k)}}$ je poté předán do Eigensolveru, kde je diagonalizován pro výpočet vlastních čísel a vlastních vektorů k rekonstrukci vlastního stavu. V této lekci promítáme a diagonalizujeme Hamiltonián pomocí balíčku qiskit-addon-sqd, který pro diagonalizaci využívá Davidsonovu metodu z PySCF.

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

Poté z dávek sesbíráme nejnižší vlastní číslo (energii) a také spočítáme průměrnou orbitální obsazenost $\text{n}$. Informace o průměrné obsazenosti se používá v kroku obnovy konfigurací k pravděpodobnostnímu opravování zašuměných konfigurací.

Dále podrobně vysvětlíme samokonzistentní smyčku obnovy konfigurací a ukážeme konkrétní příklady kódu k implementaci výše zmíněných kroků pro odhad energie základního stavu Hamiltoniánu $N_2$.

4.1 Obnovení konfigurací: přehled

Každý bit v bitovém řetězci (Slaterově determinantu) reprezentuje jeden spinový orbital. Pravá polovina bitového řetězce reprezentuje orbitaly se spinem nahoru a levá polovina reprezentuje orbitaly se spinem dolů. Hodnota 1 znamená, že orbital je obsazený elektronem, a 0 znamená, že je prázdný. Správný počet částic (jak elektronů se spinem nahoru, tak se spinem dolů) známe předem. Předpokládejme, že máme determinant $x$ s $N_x$ elektrony (tedy s $N_x$ jedničkami v bitovém řetězci). Správný počet částic je $N$. Pokud $N_x \neq N$, víme, že bitový řetězec je porušený šumem. Rutina samokonzistentního obnovení konfigurací se pokouší bitový řetězec opravit tím, že pravděpodobnostně překlopí $|N_x - N|$ bitů s využitím informací o průměrné obsazenosti orbitalů. Průměrná obsazenost orbitalu ($n$) nám říká, jak pravděpodobně bude orbital obsazený elektronem. Pokud $N_x < N$, máme méně elektronů a je potřeba překlopit některé $0$ na $1$, a naopak.

Pravděpodobnost překlopení může být $|x[i] - avg\_occupancy[i]|$ pro i-tý spinový orbital. V článku [2] autoři použili váženou pravděpodobnost překlopení pomocí modifikované funkce ReLU.

$$\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}$$

Zde $h$ určuje polohu „rohu“ funkce ReLU a parametr $\delta$ určuje hodnotu funkce ReLU v tomto rohu. Pro $\delta = 0$ se $w$ stává skutečnou funkcí ReLU a pro $\delta >0$ se stává modifikovanou ReLU. V článku autoři použili $\delta = 0.01$ a $h =$ počet alfa (nebo beta) částic/počet alfa (nebo beta) spinových orbitalů $= N/M$ (plnicí faktor).

Průměrná obsazenost orbitalů ($n$) není předem známa. První iterace odhadu základního stavu začíná s konfiguracemi, které mají pouze správné počty částic v obou spinových složkách. Po první iteraci máme odhad základního stavu a pomocí tohoto odhadu můžeme sestrojit první odhad $n$. Tento odhad $n$ se použije k obnovení konfigurací, spuštění další iterace odhadu základního stavu a samokonzistentnímu zpřesnění odhadu $n$. Proces se opakuje, dokud není splněno kritérium ukončení.

Uvažujme následující příklad pro $N = 2$ a $x = |1000\rangle$ ($N_x = 1$). Potřebujeme překlopit jednu z nul na jedničku, abychom opravili počet částic, a na výběr máme 1100, 1010 a 1001. Na základě pravděpodobnosti překlopení bude jedna z možností vybrána jako obnovená konfigurace (tedy bitový řetězec se správným počtem částic).

Předpokládejme, že v první iteraci spustíme dvě dávky a odhadnuté základní stavy z nich jsou:

$$\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}$$

Pomocí stavů výpočetní báze a jejich amplitud můžeme spočítat pravděpodobnost obsazenosti elektrony (krátce obsazenost) na spinový orbital (qubit) (pozor, pravděpodobnost = |amplituda|$^2$). Níže tabelujeme obsazenosti po qubitech pro každý bitový řetězec vyskytující se v odhadnutém základním stavu a spočítáme celkovou obsazenost orbitalů pro danou dávku. Podle konvence řazení v Qiskitu reprezentuje nejpravější bit qubit-0 (Q0) a nejlevější bit Q3.

Obsazenost (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

Obsazenost (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

Obsazenost (průměr přes dávky)

	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 (average)	0.49	0.51	0.35	0.65

Pomocí průměrné obsazenosti orbitalů spočítané výše můžeme najít pravděpodobnosti překlopení pro různé orbitaly v konfiguraci $x = \vert 1000 \rangle$. Protože orbital reprezentovaný Q3 je již obsazený a není potřeba jej překlápět, nastavíme jeho p(flip) na $0$. Pro zbývající orbitaly, které jsou neobsazené, je pravděpodobnost překlopení $\vert x[i] - \text{n}[i] \vert$. Kromě p(flip) spočítáme také váhovou pravděpodobnost spojenou s překlopením pomocí výše popsané modifikované funkce ReLU.

Pravděpodobnost překlopení ($x = \vert 1000 \rangle$, $\delta = 0.01$, $h = N/M = 2/4 = 0.50$)

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

Nakonec pomocí vážených pravděpodobností výše můžeme překlopit jeden z neobsazených orbitalů Q2, Q1 a Q0. Na základě hodnot výše bude nejpravděpodobněji překlopen Q0 a možná obnovená konfigurace může být $\vert \text{1001} \rangle$. Celý proces samokonzistentního obnovení konfigurací lze shrnout takto:

První iterace: Předpokládejme, že bitové řetězce (konfigurace nebo Slaterovy determinanty) generované kvantovým počítačem tvoří množinu $\widetilde{\chi}$, která obsahuje jak konfigurace se správným ($\widetilde{\chi}_{correct}$), tak s nesprávným ($\widetilde{\chi}_{incorrect}$) počtem částic v každém spinovém sektoru.

1. Z konfigurací z ($\widetilde{\chi}_{correct}$) se náhodně vzorkují dávky $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ vektorů pro projekci do podprostoru. Počet dávek a vzorků v každé dávce jsou uživatelem definované parametry. Čím větší je počet vzorků v každé dávce, tím větší je dimenze podprostoru a tím výpočetně náročnější je diagonalizace. Na druhou stranu příliš malý počet vzorků může opomenout nosné vektory základního stavu a vést k nesprávnému odhadu.
2. Spusť řešič vlastních stavů (tedy projekci do podprostoru a diagonalizaci) na dávkách a získej přibližné vlastní stavy $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. Z přibližných vlastních stavů sestroj první odhad $n$.

Následné iterace:

1. Pomocí $n$ oprav konfigurace se špatným počtem částic v $\widetilde{\chi}_{incorrect}$. Označme je $\widetilde{\chi}_{correct\_new}$. Pak $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ tvoří novou množinu konfigurací se správnými počty částic.
2. $\widetilde{\chi}_{R}$ se vzorkuje k vytvoření dávek $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. Řešič vlastních stavů běží s novými dávkami a generuje nové odhady základních stavů $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. Z přibližných vlastních stavů sestroj zpřesněný odhad $n$.
5. Pokud není splněno kritérium ukončení, vrať se ke kroku 2.1.

4.2 Odhad základního stavu

Nejprve převedeme počty měření na matici bitstringů a pole pravděpodobností pro následné zpracování.

Každý řádek v matici představuje jeden jedinečný bitstring. Jelikož jsou qubity v Qiskitu indexovány zprava bitstringu, sloupec 0 odpovídá qubitu N-1 a sloupec N-1 odpovídá qubitu 0, kde N je počet qubitů.

Alfa orbitaly jsou reprezentovány v rozsahu indexů sloupců (N, N/2] (pravá polovina) a beta orbitaly jsou reprezentovány v rozsahu sloupců (N/2, 0] (levá polovina).

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)

Pro tuto techniku je důležitých několik uživatelsky řízených voleb:

• iterations: Počet iterací samokonzistentní obnovy konfigurací
• n_batches: Počet dávek konfigurací, které jsou používány jednotlivými voláními řešiče vlastních stavů
• samples_per_batch: Počet jedinečných konfigurací zahrnutých v každé dávce
• max_davidson_cycles: Maximální počet Davidsonových cyklů, které každý řešič vlastních stavů provede

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_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 Diskuze výsledků

První graf ukazuje, že po několika iteracích odhadujeme energii základního stavu s přesností ~24 mH (chemická přesnost se obvykle uvažuje jako 1 kcal/mol $\approx$ 1,6 mH). Druhý graf ukazuje průměrnou obsazenost každého prostorového orbitalu po poslední iteraci. Vidíme, že jak elektrony se spinem nahoru, tak elektrony se spinem dolů v našich řešeních obsazují prvních pět orbitalů s vysokou pravděpodobností.

Ačkoli je odhadnutá energie základního stavu slušná, není v mezích chemické přesnosti ($\pm \approx 1,6$ mH). Tento rozdíl lze přičíst malé dimenzi podprostoru, kterou jsme výše použili pro projekci a diagonalizaci. Protože jsme použili samples_per_batch=500, je podprostor generován maximálně $500$ vektory, přičemž chybí vektory z nosiče základního stavu. Zvětšení parametru samples_per_batch by mělo přesnost zlepšit na úkor větších nároků na klasické výpočetní zdroje a dobu běhu.

# 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

Cvičení pro čtenáře

Postupně zvětšuj parametr samples_per_batch (například od $1000$ do $10000$ s krokem $1000$; podle toho, co dovolí paměť tvého počítače) a porovnej odhadnuté energie základního stavu.

Reference

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