SQD do estymacji energii hamiltonianu chemicznego

W tej lekcji zastosujemy SQD do oszacowania energii stanu podstawowego cząsteczki.

W szczególności omówimy następujące tematy, korzystając z $4$-etapowego podejścia wzorca Qiskit:

1. Krok 1: Odwzorowanie problemu na obwody i operatory kwantowe
   • Skonfigurowanie molecular Hamiltonianu dla $N_2$.
   • Wyjaśnienie inspirowanego chemią i przyjaznego sprzętowi lokalnego unitarnego klastra Jastrowa (LUCJ) [1]
2. Krok 2: Optymalizacja pod kątem docelowego sprzętu
   • Optymalizacja liczby bramek i układu ansatzu pod kątem wykonania na sprzęcie
3. Krok 3: Wykonanie na docelowym sprzęcie
   • Uruchomienie zoptymalizowanego obwodu na rzeczywistym QPU w celu wygenerowania próbek podprzestrzeni.
4. Krok 4: Przetwarzanie końcowe wyników
   • Wprowadzenie samouzgodnionej pętli odzyskiwania konfiguracji [2]
     • Przetwarzanie pełnego zestawu próbek ciągów bitów, przy użyciu wcześniejszej wiedzy o liczbie cząstek i średniej obsadzie orbitalnej obliczonej w ostatniej iteracji.
     • Probabilistyczne tworzenie partii podpróbek z odzyskanych ciągów bitów.
     • Rzutowanie i diagonalizowanie molecular Hamiltonianu w każdej próbkowanej podprzestrzeni.
     • Zapisanie minimalnej znalezionej energii stanu podstawowego we wszystkich partiach i aktualizacja średniej obsady orbitalnej.

W trakcie lekcji użyjemy kilku pakietów oprogramowania.

• PySCF do zdefiniowania cząsteczki i skonfigurowania Hamiltonianu.
• pakiet ffsim do skonstruowania ansatzu LUCJ.
• Qiskit do transpilacji ansatzu pod kątem wykonania na sprzęcie.
• Qiskit IBM Runtime do wykonania obwodu na QPU i zbierania próbek.
• Qiskit addon SQD do odzyskiwania konfiguracji oraz estymacji energii stanu podstawowego przy użyciu rzutowania na podprzestrzeń i diagonalizacji macierzy.

1. Odwzorowanie problemu na obwody i operatory kwantowe

Molecular Hamiltonian

Molecular Hamiltonian ma ogólną postać:

$$\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}$ to fermionowe operatory kreacji/anihilacji powiązane z $p$-tym elementem zbioru bazowego i spinem $\sigma$. $h_{pr}$ i $(pr|qs)$ to jedno- i dwuciałowe całki elektronowe. Korzystając z pySCF, zdefiniujemy cząsteczkę i obliczymy jedno- i dwuciałowe całki Hamiltonianu dla zbioru bazowego 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

W tej lekcji użyjemy transformacji Jordan-Wigner (JW), aby odwzorować funkcję falową fermionową na funkcję falową kubitową, tak aby można ją było przygotować przy użyciu obwodu kwantowego. Transformacja JW odwzorowuje przestrzeń Focka fermionów w M spatial orbitals na przestrzeń Hilberta 2M kubitów, czyli spatial orbital jest dzielony na dwa spin orbitals, jeden związany z elektronem o spinie w górę ($\alpha$), a drugi o spinie w dół ($\beta$). Spin orbital może być obsadzony lub nieobsadzony. Zwykle, gdy mówimy o liczbie orbitali, mamy na myśli liczbę orbitali przestrzennych (spatial). Liczba spin orbitals będzie dwa razy większa. W obwodach kwantowych każdy spin orbital będziemy reprezentować jednym kubitem. W ten sposób jeden zbiór kubitów będzie reprezentować orbitale o spinie w górę lub $\alpha$, a inny zbiór będzie reprezentować orbitale o spinie w dół lub $\beta$. Na przykład cząsteczka $N_2$ dla zestawu bazowego 6-31g ma $16$ spatial orbitals (czyli $16$ $\alpha$ + $16$ $\beta$ = $32$ spin orbitals). Dlatego będziemy potrzebować obwodu kwantowego z $32$ kubitami (możemy potrzebować dodatkowych kubitów pomocniczych, jak omówiono później). Kubity są mierzone w bazie obliczeniowej, aby wygenerować ciągi bitów, które reprezentują konfiguracje elektronowe lub wyznaczniki (Slatera). W całej tej lekcji będziemy używać pojęć ciągów bitów, konfiguracji i wyznaczników zamiennie. Ciągi bitów mówią nam o obsadzeniu elektronów w spin orbitals: $1$ na pozycji bitu oznacza, że odpowiadający spin orbital jest obsadzony, podczas gdy $0$ oznacza, że spin orbital jest pusty. Ponieważ problemy struktury elektronowej zachowują liczbę cząstek, tylko stała liczba spin orbitals musi być obsadzona. Cząsteczka $N_2$ ma $5$ elektronów o spinie w górę ($\alpha$) i $5$ elektronów o spinie w dół ($\beta$). Dlatego każdy ciąg bitów reprezentujący orbitale $\alpha$ i $\beta$ musi mieć po pięć $1\text{-nek}$ dla cząsteczki $N_2$.

1.1 Obwód kwantowy do generowania próbek: ansatz LUCJ

W tej lekcji użyjemy ansatz local unitary coupled cluster Jastrow (LUCJ) \[1\] do przygotowania stanu kwantowego i późniejszego próbkowania. Najpierw wyjaśnimy różne bloki budulcowe pełnego ansatz UCJ oraz przybliżenia wprowadzone w jego wersji lokalnej. Następnie, używając pakietu ffsim, skonstruujemy ansatz LUCJ i zoptymalizujemy go za pomocą transpilera qiskit do wykonania sprzętowego.

Ansatz UCJ ma następującą postać (dla iloczynu $L$ warstw lub powtórzeń operatora UCJ).

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

gdzie $\vert \Phi_{0} \rangle$ jest stanem referencyjnym, zazwyczaj przyjmowanym jako stan Hartree-Fock (HF). Ponieważ stan Hartree-Fock jest zdefiniowany jako mający obsadzone orbitale o najniższych numerach, przygotowanie stanu HF będzie obejmować zastosowanie bramek X do ustawienia kubitów odpowiadających obsadzonym orbitalom na jeden. Na przykład blok przygotowania stanu HF dla 4 spatial orbitals oraz 2 elektronów o spinie w górę i 2 elektronów o spinie w dół może wyglądać następująco: Pojedyncze powtórzenie operatora UCJ ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ składa się z ewolucji diagonalnej Coulomba ($e^{iJ^{(\mu)}}$) otoczonej rotacjami orbital ($e^{K^{(\mu)}}$ oraz $e^{-K^{(\mu)}}$). Bloki rotacji orbital działają na pojedynczym rodzaju spinu ($\alpha$ (spin w górę)/$\beta$ (spin w dół)). Dla każdego rodzaju elektronu rotacja orbital składa się z warstwy jednokubitowych bramek $R_{z}$, po których następuje sekwencja dwukubitowych bramek rotacji Givensa (bramki $XX + YY$).

Bramki dwukubitowe działają na sąsiednich spin-orbitals (najbliższych sąsiadujących kubitach) i dlatego są realizowalne na jednostkach QPU IBM® bez potrzeby stosowania bramek SWAP. $e^{iJ^{(\mu)}}$, znany również jako operator diagonalny Coulomba, składa się z trzech bloków. Dwa z nich działają na sektory o tym samym spinie ($e^{iJ_{\alpha \alpha}^{(\mu)}}$ oraz $e^{iJ_{\beta \beta}^{(\mu)}}$), a jeden działa pomiędzy dwoma sektorami spinowymi ($e^{iJ_{\alpha \beta}^{(\mu)}}$).

Wszystkie bloki w $e^{iJ^{(\mu)}}$ składają się z bramek liczba-liczba $U_{nn}(\phi)$ [1]. Bramka $U_{nn}(\phi)$ może być dalej rozłożona na bramkę $R_{ZZ}(\frac{\phi}{2})$, po której następują dwie jednokubitowe bramki $Rz(-\frac{\phi}{2})$ działające na dwóch oddzielnych kubitach.

Komponenty o tym samym spinie ($J_{\alpha \alpha}$ i $J_{\beta \beta}$) mają bramki $U_{nn}$ pomiędzy wszystkimi możliwymi parami kubitów. Jednakże, ponieważ nadprzewodzące jednostki QPU mają ograniczoną łączność, kubity muszą być zamieniane, aby zrealizować bramki pomiędzy niesąsiadującymi kubitami.

Rozważmy na przykład następujący blok $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (lub $e^{iJ_{\beta \beta}^{(\mu)}}$) dla $N = 4$ spatial orbitals. Przy liniowej łączności kubitów ostatnie trzy bramki nie są bezpośrednio realizowalne, ponieważ działają pomiędzy niesąsiadującymi kubitami (na przykład Q0 i Q2 nie są bezpośrednio połączone). Dlatego potrzebujemy bramek SWAP, aby je zbliżyć (następny rysunek pokazuje przykład z $3$ bramkami SWAP). Następnie, $J_{\alpha \beta}$ realizuje bramki pomiędzy orbitalami o tym samym indeksie z różnych sektorów spinowych (na przykład pomiędzy $0\alpha$ i $0\beta$). Podobnie, jeśli kubity nie są fizycznie sąsiednie na QPU, te bramki również będą wymagały SWAP-ów. Z powyższej dyskusji wynika, że ansatz UCJ napotyka pewne przeszkody przy wykonaniu sprzętowym, ponieważ wymaga bramek SWAP z powodu interakcji niesąsiadujących kubitów. Lokalny wariant ansatz UCJ, LUCJ, rozwiązuje to wyzwanie, usuwając niektóre $U_{nn}$ z operatora diagonalnego Coulomba.

W blokach tych samych rodzajów elektronów, $J_{\alpha \alpha}$ i $J_{\beta \beta}$, zachowujemy tylko bramki $U_{nn}$ kompatybilne z łącznością najbliższych sąsiadów i usuwamy bramki pomiędzy niesąsiadującymi kubitami w wersji LUCJ. Następujący rysunek pokazuje blok LUCJ po usunięciu bramek niesąsiadujących. Następnie, wersja LUCJ bloku $J_{\alpha \beta}$, który działa pomiędzy różnymi rodzajami elektronów, może przybierać różne kształty w zależności od topologii urządzenia.

Również tutaj wersja LUCJ pozbywa się niekompatybilnych bramek. Poniższy rysunek pokazuje warianty bloku $J_{\alpha \beta}$ dla różnych topologii kubitów, w tym siatki, heksagonalnej, heavy-hex i liniowej.

• Kwadratowa: możemy mieć bramki $U_{nn}$ pomiędzy wszystkimi orbitalami $\alpha$ i $\beta$ bez żadnych SWAP-ów, a zatem nie musimy usuwać żadnych bramek $U_{nn}$.
• Heavy-hex: Interakcje $\alpha$-$\beta$ są zachowywane pomiędzy co $4$-tym indeksowanym (takim jak 0-ty, 4-ty i 8-my) spin orbitals i są mediowane przez ancilla, czyli potrzebujemy kubitów pomocniczych (ancilla) pomiędzy łańcuchami liniowymi reprezentującymi orbitale $\alpha$ i $\beta$. Ten układ wymaga ograniczonej liczby SWAP-ów.
• Heksagonalna: Co drugi orbital, taki jak 0-ty, 2-gi i 4-ty indeksowany orbital, staje się najbliższymi sąsiadami, gdy $\alpha$ i $\beta$ są rozłożone w dwóch sąsiadujących łańcuchach liniowych.
• Liniowa: Tylko jeden orbital $\alpha$ i jeden orbital $\beta$ są połączone, co oznacza, że blok $J_{\alpha \beta}$ będzie miał tylko jedną bramkę. Chociaż usunięcie bramek z ansatz UCJ w celu skonstruowania wersji LUCJ czyni go bardziej kompatybilnym ze sprzętem, ansatz traci nieco ekspresywności. Dlatego przy użyciu ansatz LUCJ może być potrzebnych więcej powtórzeń ($L$) zmodyfikowanego operatora UCJ.

1.2 Inicjalizacja ansatz LUCJ

LUCJ jest sparametryzowanym ansatz i musimy zainicjalizować parametry przed wykonaniem sprzętowym. Jednym ze sposobów inicjalizacji ansatz jest użycie amplitud t1 i t2 z klasycznej metody coupled cluster singles and doubles (CCSD), gdzie amplitudy t1 są współczynnikami operatorów pojedynczych wzbudzeń, a t2 amplitudes są dla operatorów podwójnych wzbudzeń.

Zauważ, że chociaż inicjalizacja ansatz LUCJ amplitudami t1 i t2 daje przyzwoite wyniki, parametry ansatz mogą wymagać dalszej optymalizacji.

# 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 Konstruowanie ansatz LUCJ przy użyciu ffsim

Użyjemy pakietu ffsim do utworzenia i zainicjalizowania ansatz amplitudami t1 i t2 obliczonymi powyżej. Ponieważ nasza cząsteczka ma stan Hartree-Fock z zamkniętą powłoką, użyjemy wariantu ansatz UCJ zbalansowanego pod względem spinu, UCJOpSpinBalanced.

Ponieważ sprzęt IBM ma topologię heavy-hex, przyjmiemy wzorzec zig-zag użyty w [1] i wyjaśniony powyżej dla interakcji kubitów. W tym wzorcu orbitale (kubity) o tym samym spinie są połączone topologią liniową (czerwone i niebieskie kółka). Ze względu na topologię heavy-hex orbitale dla różnych spinów mają połączenia pomiędzy co 4-tym orbital, czyli 0-tym, 4-tym, 8-mym itd. (fioletowe kółka).

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 z powtórzonymi warstwami można zoptymalizować, scalając niektóre sąsiednie bloki. Rozważ przypadek dla n_reps=2. Dwa bloki rotacji orbital w środku mogą być scalone w pojedynczy blok rotacji orbital. Pakiet ffsim ma pass manager o nazwie ffsim.qiskit.PRE_INIT do optymalizacji obwodu poprzez scalanie takich sąsiednich bloków.

2. Optymalizacja pod docelowy sprzęt

Najpierw pobieramy wybrany backend. Zoptymalizujemy nasz obwód pod ten backend, a następnie uruchomimy zoptymalizowany obwód na tym samym backendzie, aby wygenerować próbki dla podprzestrzeni.

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

Następnie zalecamy poniższe kroki w celu optymalizacji ansatzu i uczynienia go zgodnym ze sprzętem.

• Wybierz fizyczne kubity (initial_layout) z docelowego sprzętu, które są zgodne ze wzorcem zygzaka (dwa liniowe łańcuchy z kubitem pomocniczym pomiędzy nimi) opisanym powyżej. Ułożenie kubitów w tym wzorcu prowadzi do wydajnego, zgodnego ze sprzętem obwodu z mniejszą liczbą bramek.
• Wygeneruj etapowy pass manager za pomocą funkcji generate_preset_pass_manager z Qiskit z wybranym backend i initial_layout.
• Ustaw etap pre_init swojego etapowego pass managera na ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT zawiera przebiegi transpilera Qiskit, które rozkładają bramki na rotacje orbitalne, a następnie łączą te rotacje orbitalne, co skutkuje mniejszą liczbą bramek w końcowym obwodzie.
• Uruchom pass manager na swoim obwodzie.

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. Wykonanie na docelowym sprzęcie

Po zoptymalizowaniu obwodu pod kątem wykonania sprzętowego jesteśmy gotowi uruchomić go na docelowym sprzęcie i zebrać próbki do oszacowania energii stanu podstawowego. Ponieważ mamy tylko jeden obwód, użyjemy trybu wykonania Job z Qiskit Runtime i wykonamy nasz obwód.

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. Przetwarzanie końcowe wyników

Część post-processing przepływu pracy SQD można podsumować za pomocą następującego diagramu.

Próbkowanie układu LUCJ w bazie obliczeniowej generuje pulę zaszumionych konfiguracji $\tilde{\mathcal{\chi}}$, które są używane w procedurze post-processing. Obejmuje ona metodę zwaną (szczegóły omówione później) odzyskiwaniem konfiguracji (ang. configuration recovery), która probabilistycznie koryguje konfiguracje z błędnymi liczbami elektronów. Konfiguracje tylko z poprawnymi liczbami elektronów $\tilde{\mathcal{\chi}}_{R}$ są następnie podpróbkowywane i dystrybuowane do wielu partii (batches) w oparciu o częstość występowania każdej unikalnej konfiguracji. Każda partia próbek definiuje podprzestrzeń ($\mathcal{S^{(k)}}$). Następnie hamiltonian molekularny, $H$, jest rzutowany na podprzestrzenie:

$$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żdy zrzutowany hamiltonian $H_{\mathcal{S}^{(k)}}$ jest następnie przekazywany do solvera wartości własnych (Eigensolver), gdzie jest diagonalizowany w celu obliczenia wartości własnych i wektorów własnych, aby zrekonstruować eigenstate. W tej lekcji rzutujemy i diagonalizujemy hamiltonian za pomocą pakietu qiskit-addon-sqd, który do diagonalizacji wykorzystuje metodę Davidsona z PySCF.

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

Następnie zbieramy najniższą wartość własną (energię) z partii, a także obliczamy średnie obsadzenie orbitali, $\text{n}$. Informacja o średnim obsadzeniu jest używana w kroku odzyskiwania konfiguracji do probabilistycznej korekty konfiguracji z szumem.

W dalszej części szczegółowo wyjaśniamy samozgodną pętlę odzyskiwania konfiguracji i pokazujemy konkretne przykłady kodu implementujące wspomniane wyżej kroki w celu oszacowania energii stanu podstawowego hamiltonianu $N_2$.

4.1 Odzyskiwanie konfiguracji: przegląd

Każdy bit w bitstring (wyznacznik Slatera) reprezentuje orbital spinowy. Prawa połowa bitstring reprezentuje orbitale ze spinem w górę, a lewa połowa reprezentuje orbitale ze spinem w dół. 1 oznacza, że orbital jest obsadzony przez elektron, a 0 oznacza, że orbital jest pusty. Znamy a priori poprawną liczbę cząstek (zarówno elektronów o spinie w górę, jak i o spinie w dół). Załóżmy, że mamy wyznacznik $x$ z $N_x$ elektronami (tzn. w bitstring znajduje się $N_x$ jedynek). Poprawna liczba cząstek to $N$. Jeśli $N_x \neq N$, to wiemy, że bitstring jest uszkodzony przez szum. Samozgodna procedura odzyskiwania konfiguracji próbuje skorygować bitstring poprzez probabilistyczne odwrócenie $|N_x - N|$ bitów, wykorzystując informacje o średnim obsadzeniu orbitali. Średnie obsadzenie orbitali ($n$) mówi nam, jak prawdopodobne jest, że orbital będzie obsadzony przez elektron. Jeśli $N_x < N$, mamy mniej elektronów i musimy odwrócić niektóre zera na jedynki i vice versa.

Prawdopodobieństwo odwrócenia może wynosić $|x[i] - avg\_occupancy[i]|$ dla i-tego orbitala spinowego. W [2] autorzy użyli ważonego prawdopodobieństwa odwrócenia, wykorzystując zmodyfikowaną funkcję 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}$$

Tutaj $h$ określa położenie „narożnika” funkcji ReLU, a parametr $\delta$ definiuje wartość funkcji ReLU w narożniku. Dla $\delta = 0$ $w$ staje się prawdziwą funkcją ReLU, a dla $\delta >0$ staje się zmodyfikowaną funkcją ReLU. W artykule autorzy użyli $\delta = 0.01$ oraz $h =$ liczba cząstek alfa (lub beta)/liczba orbitali spinowych alfa (lub beta) $= N/M$ (współczynnik wypełnienia).

Średnie obsadzenie orbitali ($n$) nie jest znane a priori. Pierwsza iteracja estymacji stanu podstawowego rozpoczyna się z konfiguracjami zawierającymi tylko poprawne liczby cząstek w obu gatunkach spinowych. Po pierwszej iteracji mamy oszacowanie stanu podstawowego, a korzystając z tego oszacowania, możemy skonstruować pierwsze przybliżenie $n$. To przybliżenie $n$ jest używane do odzyskiwania konfiguracji, uruchomienia kolejnej iteracji estymacji stanu podstawowego i samozgodnego udoskonalania przybliżenia $n$. Proces powtarza się, aż do spełnienia kryterium stopu.

Rozważmy następujący przykład dla $N = 2$ i $x = |1000\rangle$ ($N_x = 1$). Musimy odwrócić jedno z zer na jedynkę, aby skorygować liczbę cząstek, a możliwości to 1100, 1010 i 1001. W oparciu o prawdopodobieństwo odwrócenia, jedna z możliwości zostanie wybrana jako odzyskana konfiguracja (czyli bitstring z poprawną liczbą cząstek).

Załóżmy, że w pierwszej iteracji uruchamiamy dwie partie, a oszacowane stany podstawowe z nich to:

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

Korzystając ze stanów bazy obliczeniowej i ich amplitud, możemy obliczyć prawdopodobieństwo obsadzenia elektronami (w skrócie obsadzenia) na orbital spinowy (Qubit) (zauważ, że prawdopodobieństwo = |amplituda|$^2$). Poniżej przedstawiamy w tabeli obsadzenia poszczególnych kubitów dla każdego bitstring występującego w oszacowanym stanie podstawowym i obliczamy całkowite obsadzenie orbitali dla partii. Zauważ, że zgodnie z konwencją porządkowania w Qiskit, skrajnie prawy bit reprezentuje kubit-0 (Q0), a skrajnie lewy bit reprezentuje Q3.

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

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

Obsadzenie (średnia z partii)

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

Korzystając z obliczonego powyżej średniego obsadzenia orbitali, możemy znaleźć prawdopodobieństwa odwrócenia dla różnych orbitali w konfiguracji $x = \vert 1000 \rangle$. Ponieważ orbital reprezentowany przez Q3 jest już obsadzony i nie trzeba go odwracać, ustawiamy jego p(flip) na $0$. Dla pozostałych orbitali, które są nieobsadzone, prawdopodobieństwo odwrócenia wynosi $\vert x[i] - \text{n}[i] \vert$ dla każdego z nich. Oprócz p(flip) obliczamy również wagę prawdopodobieństwa związanego z odwróceniem, korzystając z opisanej wyżej zmodyfikowanej funkcji ReLU.

Prawdopodobieństwo odwrócenia ($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

Wreszcie, korzystając z powyższych ważonych prawdopodobieństw, możemy odwrócić jeden z nieobsadzonych orbitali Q2, Q1 lub Q0. W oparciu o powyższe wartości, Q0 zostanie odwrócony najbardziej prawdopodobnie, a możliwą odzyskaną konfiguracją może być $\vert \text{1001} \rangle$. Kompletny samozgodny proces odzyskiwania konfiguracji można podsumować następująco:

Pierwsza iteracja: Załóżmy, że bitstring (konfiguracje lub wyznaczniki Slatera) wygenerowane przez komputer kwantowy tworzą zbiór $\widetilde{\chi}$, który obejmuje zarówno konfiguracje z poprawną ($\widetilde{\chi}_{correct}$), jak i niepoprawną ($\widetilde{\chi}_{incorrect}$) liczbą cząstek w każdym sektorze spinowym.

1. Konfiguracje z ($\widetilde{\chi}_{correct}$) są losowo próbkowane w celu utworzenia partii $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ wektorów do rzutowania na podprzestrzeń. Liczba partii i próbek w każdej partii są parametrami zdefiniowanymi przez użytkownika. Im większa liczba próbek w każdej partii, tym większy wymiar podprzestrzeni i bardziej wymagająca obliczeniowo staje się diagonalizacja. Z drugiej strony, zbyt mała liczba próbek może pominąć wektory wspierające stan podstawowy i prowadzić do błędnego oszacowania.
2. Uruchom solver eigenstate (tzn. rzutowanie na podprzestrzeń i diagonalizację) na partiach i uzyskaj przybliżone eigenstates. $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. Na podstawie przybliżonych eigenstates skonstruuj pierwsze przybliżenie $n$.

Kolejne iteracje:

1. Korzystając z $n$, skoryguj konfiguracje z błędną liczbą cząstek w $\widetilde{\chi}_{incorrect}$. Załóżmy, że nazwiemy je $\widetilde{\chi}_{correct\_new}$. Wtedy $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ tworzy nowy zbiór konfiguracji z poprawnymi liczbami cząstek.
2. $\widetilde{\chi}_{R}$ jest próbkowany w celu utworzenia partii $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. Solver eigenstate działa z nowymi partiami i generuje nowe oszacowania stanów podstawowych $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. Na podstawie przybliżonych eigenstates skonstruuj udoskonalone przybliżenie $n$.
5. Jeśli kryterium stopu nie jest spełnione, wróć do kroku 2.1.

4.2 Estymacja stanu podstawowego

Najpierw przekształcimy counts w macierz bitstring i tablicę prawdopodobieństw do post-processingu.

Każdy wiersz w macierzy reprezentuje jeden unikalny bitstring. Ponieważ kubity w Qiskit są indeksowane od prawej strony bitstring, kolumna 0 reprezentuje kubit N-1, a kolumna N-1 reprezentuje kubit 0, gdzie N jest liczbą kubitów.

Orbitale alfa są reprezentowane w zakresie indeksów kolumn (N, N/2] (prawa połowa), a orbitale beta są reprezentowane w zakresie kolumn (N/2, 0] (lewa połowa).

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)

Istnieje kilka opcji kontrolowanych przez użytkownika, które są ważne dla tej techniki:

• iterations: Liczba iteracji self-consistent configuration recovery
• n_batches: Liczba partii konfiguracji używanych przez różne wywołania eigenstate solvera
• samples_per_batch: Liczba unikalnych konfiguracji uwzględnionych w każdej partii
• max_davidson_cycles: Maksymalna liczba cykli Davidsona uruchamianych przez każdy 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 Omówienie wyników

Pierwszy wykres pokazuje, że po kilku iteracjach szacujemy ground state energy z dokładnością do ~24 mH (dokładność chemiczna jest zwykle przyjmowana jako 1 kcal/mol $\approx$ 1,6 mH). Drugi wykres pokazuje średnie obsadzenie każdego orbitalu przestrzennego po ostatniej iteracji. Widzimy, że zarówno elektrony o spinie w górę, jak i w dół z wysokim prawdopodobieństwem zajmują w naszych rozwiązaniach pierwsze pięć orbitali.

Chociaż oszacowana ground state energy jest przyzwoita, nie mieści się w granicach dokładności chemicznej ($\pm \approx 1,6$ mH). Tę rozbieżność można przypisać małemu wymiarowi podprzestrzeni, który wykorzystaliśmy powyżej do projekcji i diagonalizacji. Ponieważ użyliśmy samples_per_batch=500, podprzestrzeń jest rozpięta przez maksymalnie $500$ wektorów, co oznacza brak wektorów należących do nośnika stanu podstawowego. Zwiększenie parametru samples_per_batch powinno poprawić dokładność kosztem większych klasycznych zasobów obliczeniowych i dłuższego czasu działania.

# 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

Ćwiczenie dla czytelnika

Stopniowo zwiększaj parametr samples_per_batch (na przykład od $1000$ do $10000$ z krokiem $1000$; na tyle, na ile pozwala pamięć Twojego komputera) i porównaj oszacowane energie stanu podstawowego.

Bibliografia

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