Algorytm SqDRIFT do oszacowania stanu podstawowego

Szacowany czas: 180 sekund na procesorze Heron r3 (UWAGA: To jest tylko oszacowanie. Twój rzeczywisty czas wykonania może się różnić.)

Szukasz wersji C++?

Ten samouczek wykorzystuje Python. Implementację w C++, wraz z kodem źródłowym i instrukcjami budowania, znajdziesz w samouczku C++ SqDRIFT.

Efekty uczenia się

• Dowiedz się, jak tworzyć obwody o mniejszej głębokości w porównaniu z trotteryzacją

• Przejdź przez kompletny przepływ pracy do oszacowania stanu podstawowego za pomocą qDRIFT i SQD

• Dowiedz się, jak używać qiskit-fermions w połączeniu z innymi dodatkami Qiskit, aby zaimplementować taki przepływ pracy

Ten samouczek jest przedstawiony jako notatnik Python w celach dydaktycznych.

Wymagania wstępne

• Przeczytaj przegląd Diagonalizacji kwantowej opartej na próbkowaniu (SQD)

• Przeczytaj lekcję Krylowowskiej diagonalizacji kwantowej opartej na próbkowaniu (SKQD)

Kontekst

SqDRIFT to wariant SKQD, który zastępuje potrzebę wyboru ansatzu, z którego próbkuje się ciągi bitowe, zespołem obwodów ewolucji czasowej skonstruowanych bezpośrednio z docelowego Hamiltonianu. Osiąga się to poprzez subsamplowanie mniejszych operatorów ewolucji czasowej z Hamiltonianu na podstawie jego współczynników, co znane jest jako metoda trotteryzacji qDRIFT.

Ten samouczek wykorzystuje Qiskit Fermions do tworzenia bardziej naturalnych obwodów fermionowych dla algorytmu qDRIFT, a następnie wykorzystuje przebiegi fermionowego rozmieszczenia i syntezy przed wprowadzeniem obwodów do tradycyjnego potoku Qiskit do wykonania sprzętowego.

Niech Hamiltonian ma postać:

$$H = \sum_{i=1}^{N} c_i h_i$$

gdzie, bez utraty ogólności, wymagamy $c_i > 0$ oraz aby największa wartość własna $h_i$ była równa, co do wartości bezwzględnej, $1$. Każdy oznaczony lub zespolony współczynnik jest wchłaniany przez $h_i$, więc współczynniki $c_i$ są ściśle dodatnimi wagami, podczas gdy $h_i$ niosą kierunek każdego członu. Tutaj $N$ to liczba członów (lub, po grupowaniu, liczba grup) w Hamiltonianie; jest to właściwość Hamiltonianu i różni się od liczby operatorów próbkowanych do pojedynczego obwodu, oznaczanej dalej jako $n$.

Algorytm qDRIFT następnie realizuje, dla docelowego czasu $t$, pewien operator $V_k$, gdzie $k$ przebiega od $1 \cdots K$ i oznacza $k$-ty obwód SqDRIFT, zdefiniowany jako:

$$V_k = \prod_{j=1}^{n} e^{-i h_{k_j} \lambda t / n }$$

Tutaj $n$ to liczba próbkowanych operatorów na obwód, a $K$ to liczba obwodów w zespole. Iloczyn przebiega po $n$ losowaniach, a nie po wszystkich $N$ członach Hamiltonianu, i ponieważ człony są losowane ze zwracaniem, ten sam $h_i$ może pojawić się więcej niż raz w pojedynczym $V_k$.

Wielkość:

$$\lambda = \sum_{i=1}^{N} c_i$$

jest normą $L_1$ współczynników, więc każdy z $n$ kroków ewoluuje przez ten sam czas trwania $\lambda t / n$, niezależnie od tego, który człon został wylosowany. Jednorodność kąta kroku jest charakterystyczną cechą qDRIFT: współczynnik wpływa na wynik poprzez to, jak często jego człon jest losowany, a nie poprzez to, jak daleko ten człon jest obracany. Indeksy są próbkowane z rozkładu:

$$P[k_i] = \frac{c_i}{\lambda}$$

więc seria $(k_1, \ldots, k_n)$ jest losową sekwencją indeksów członów wylosowaną z tego rozkładu. Ponieważ $c_i$ są dodatnie i sumują się do $\lambda$, jest to znormalizowany rozkład prawdopodobieństwa, a wartość oczekiwana wynikowego kanału po losowych losowaniach przybliża ewolucję pod $H$, z błędem, który maleje wraz ze wzrostem $n$. Zauważ, że błąd przybliżenia zależy od $\lambda$, a nie od liczby członów $N$.

(Artykuł o SqDRIFT zapisuje liczbę członów jako $\mathcal{N}$, a długość sekwencji jako $N$; tutaj używamy $N$ i $n$, aby wyraźnie odróżnić te dwie wielkości.)

Ten samouczek pokazuje, jak wygenerować zespół takich losowych obwodów. Po utworzeniu tych obwodów, podobnie jak przy tworzeniu podprzestrzeni Kryłowa dla różnych operatorów, próbkujemy ciągi bitowe z wielu takich operatorów z różnymi parametrami czasowymi. Zapewnia to wyższe nakładanie się wektorów stanu podstawowego z próbkowanymi ciągami bitowymi.

Wymagania

Przed rozpoczęciem tego samouczka upewnij się, że masz zainstalowane

• środowisko wirtualne Python (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (Zwróć uwagę, że nazwa jest w liczbie mnogiej)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Wszystkie wymagane pakiety możesz zainstalować za pomocą:

pip install "qiskit~=2.5" "qiskit-fermions==0.1.0" qiskit-aer qiskit-ibm-runtime qiskit-addon-sqd pyscf numpy

Konfiguracja

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

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

# Qiskit core
from qiskit import transpile
from qiskit.primitives import BitArray

# Qiskit Aer
from qiskit_aer import AerSimulator

# IBM Quantum Compute Service
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler

# Qiskit Fermions
from qiskit_fermions.operators.library import FCIDump
from qiskit_fermions.operators import FermionOperator
from qiskit_fermions.operators.terms.filtering import filter_diagonal_terms
from qiskit_fermions.operators.terms.grouping import (
    group_terms_by_electronic_structure,
)
from qiskit_fermions.operators.terms.ordering import canonical_order
from qiskit_fermions.circuit import FermionicCircuit
from qiskit_fermions.circuit.library import Evolution
from qiskit_fermions.transpiler import FermionicPassManager
from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
from qiskit_fermions.transpiler.passes import QDriftTrotterization
from qiskit_fermions.circuit.library import InitializeModes

# Qiskit addon SQD
from qiskit_addon_sqd.fermion import (
    diagonalize_fermionic_hamiltonian,
    SCIResult,
)

Przykład symulatora

Krok 1: Zmapuj dane wejściowe klasyczne na problem kwantowy

Odczyt i przygotowanie FCIDump

W tym samouczku wczytamy Hamiltonian struktury elektronowej dla azotu (N2). Istnieją też inne sposoby tworzenia operatorów fermionowych. Zapoznaj się z dokumentacją pod adresem qiskit_fermions.operators.library.

O tym pliku FCIDump. Plik N2_sto_3g opisuje cząsteczkę azotu ($N_2$) w minimalnej bazie STO-3G przy odległości międzyatomowej 1,09 $\AA$, czyli eksperymentalnej równowagowej długości wiązania. W jego nagłówku zadeklarowano NORB=10, NELEC=14 oraz MS2=0: 10 orbitali przestrzennych (a więc 20 orbitali spinowych i 20 kubitów w mapowaniu Jordana-Wignera), 14 elektronów w stanie singletowym spinu, czyli siedem elektronów $\alpha$ i siedem $\beta$. Wszystkim orbitalom przypisano etykietę symetrii 1, co oznacza, że nie wykorzystuje się symetrii grupy punktowej. Ponieważ jest to pełnoprzestrzenny zrzut STO-3G, żadne orbitale nie są zamrożone, a przestrzeń korelacji jest wystarczająco mała, aby można było obliczyć klasycznie dokładną energię referencyjną FCI do porównania, co pokazano w następnej komórce.

Równoważny plik można wygenerować ponownie za pomocą PySCF:

from pyscf import gto, scf, tools

mol = gto.M(atom="N 0 0 0; N 0 0 1.09", basis="sto-3g", symmetry=False)
mf = scf.RHF(mol).run()
tools.fcidump.from_scf(mf, "N2_sto_3g")

Ponieważ całki zależą od zbieżnych orbitali SCF, ponownie wygenerowany plik może różnić się od dostarczonego pliku fazą lub kolejnością orbitali; nie wpływa to na energie całkowite.

Uzyskiwanie pliku. Znajdź plik FCIDump w tym repozytorium GitHub. Możesz uruchomić poniższą komórkę, aby pobrać go do lokalizacji oczekiwanej przez resztę tego samouczka.

Najpierw użyjemy cisolver dostarczonego przez pyscf, aby uzyskać energię referencyjną. Jest to prawdziwa energia stanu podstawowego cząsteczki, z którą pracujemy. W tym celu najpierw zadeklarujemy norb i nelec, czyli odpowiednio liczbę orbitali i liczbę elektronów. Następnie zadeklarujemy h1e i h2e, czyli odpowiednio całki jedno- i dwuelektronowe. Wszystkie one zostaną później wykorzystane również w SQD.

import os
from urllib.request import urlopen

# The FCIDump is stored with this tutorial in the Qiskit documentation repository.
FCIDUMP_URL = "https://raw.githubusercontent.com/Qiskit/documentation/main/docs/tutorials/assets/sqdrift/fcidump_files/N2_sto_3g"
FCIDUMP_PATH = "assets/sqdrift/fcidump_files/N2_sto_3g"

if not os.path.exists(FCIDUMP_PATH):
    os.makedirs(os.path.dirname(FCIDUMP_PATH), exist_ok=True)
    with urlopen(FCIDUMP_URL) as response:
        contents = response.read()
    with open(FCIDUMP_PATH, "wb") as f:
        f.write(contents)
    print(f"Downloaded FCIDump to {FCIDUMP_PATH}")
else:
    print(f"Using existing FCIDump at {FCIDUMP_PATH}")

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha

Wczytywanie hamiltonianu

Mając gotowe niezbędne dane, wczytujemy hamiltonian z pliku FCI w formacie zgodnym z qiskit-fermions

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

Przepływy pracy fermionowe z qiskit-fermions

Najpierw zmapujemy hamiltonian na model obwodu fermionowego za pomocą qiskit-fermions, który dostarcza przebiegi transpilera i bramki specyficzne dla obwodów fermionowych. Będą one później używane przed tradycyjnymi przebiegami transpilera Qiskit w tym przepływie pracy.

Grupowanie wyrazów

Aby zapewnić powtarzalność wyników, najpierw używamy canonical_order, aby posortować wyrazy wyłącznie na podstawie ich struktury. Kolejność operatorów na liście canon jest zatem ustalona. Zapewnia to powtarzalność tworzonych operatorów, ponieważ przebieg QDriftTrotterization, którego użyjemy w dalszej części, próbkuje losowe indeksy, aby utworzyć operatory qDRIFT.

W tym kroku wykorzystujemy liczne symetrie obecne w hamiltonianie struktury elektronowej, grupując powiązane wyrazy o identycznych współczynnikach. Choć zmienia to rozkład współczynników operatora, z którego próbkuje protokół qDRIFT, nie wpływa to na jego gwarancje zbieżności. Co ważne, grupowanie wyrazów powiązanych symetrią prowadzi do korzystnego znoszenia się wyrazów Pauliego oraz do ogólnie krótszej głębokości obwodu podczas ewolucji czasowej stanu pod ich działaniem.

qiskit-fermions udostępnia funkcję group_terms_by_electronic_structure, która wykonuje za nas to grupowanie.

Zwróć uwagę, że group_terms_by_electronic_structure zakłada wyrazy w porządku normalnym.

Filtrowanie wyrazów diagonalnych

Usuwamy wyrazy diagonalne z hamiltonianu używanego do generowania obwodów, tak aby $n$ slotów próbkowania qDRIFT było przeznaczonych na wyrazy przenoszące populację między konfiguracjami. Takie wyrazy najlepiej odfiltrować z hamiltonianu na tym etapie, przed skonstruowaniem bramki Evolution w następnym kroku.

Wyrazy, o których mowa, to te diagonalne w bazie liczby obsadzeń, czyli iloczyny operatorów liczby $a^\dagger_i a_i$. Trzy rodzaje wyrazów wpisują się w ten opis:

• stała przesunięcia energii, iloczyn zera operatorów liczby, którego ewolucja czasowa wnosi jedynie fazę globalną;

• pojedyncze operatory liczby $n_i$, których ewolucja czasowa sprowadza się do jednokubitowych rotacji $Z$;

• iloczyny wyższego rzędu, takie jak $n_i n_j$.

Same w sobie żadne z nich nie przenoszą populacji między konfiguracjami liczby obsadzeń; działają jedynie na fazy już obecnych konfiguracji. Nie są jednak obojętne: te względne fazy wpływają na interferencję generowaną przez wyrazy wzbudzeń później w obwodzie, więc ich odfiltrowanie zmienia faktycznie generowaną ewolucję i może zmienić rozkład próbkowania. Jest to celowe przybliżenie na etapie generowania obwodu, mające na celu skupienie próbkowania na wyrazach wzbudzeń, a nie krok, który pozostawia próbkowany rozkład bez zmian. W przeciwieństwie do grupowania symetrii opisanego powyżej, które pozostawia nienaruszone gwarancje zbieżności qDRIFT, ten filtr zmienia operator poddawany ewolucji. Obwody nie przybliżają więc już ewolucji pod działaniem pełnego hamiltonianu, a granice błędu qDRIFT dotyczą operatora po filtracji, a nie oryginalnego. Jest to tutaj akceptowalne, ponieważ obwody są jedynie heurystyką próbkowania używaną do proponowania konfiguracji: żaden wyraz nie jest tracony z samego oszacowania energii, ponieważ filtr dotyczy wyłącznie hamiltonianu używanego do budowy obwodów, natomiast późniejsza diagonalizacja klasyczna wykorzystuje pełny hamiltonian, wraz z wyrazami diagonalnymi. Dokładność SQD zależy od tego kroku klasycznego, który pozostaje wariacyjny w próbkowanej podprzestrzeni niezależnie od sposobu proponowania konfiguracji.

Funkcja filter_diagonal_terms() usuwa takie wyrazy z operatora w miejscu. Identyfikuje je na podstawie ich struktury w porządku normalnym — multizbioru trybów kreacji odpowiadającego multizbiorowi trybów anihilacji — więc jest ważna tylko dla operatora już uporządkowanego normalnie. To założenie nie jest sprawdzane w czasie wykonania.

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

5060

Teraz, gdy pogrupowaliśmy wyrazy w hamiltonianie, ustalimy następujące parametry, aby wygenerować zespół obwodów:

• Liczbę obwodów do wygenerowania: num_circuits
• Długość każdego obwodu wyrażoną w grupach wzbudzeń: num_exc
• Współczynnik dla różnych czasów ewolucji: times

Tworzenie obwodów fermionowych

Teraz utworzymy obwody fermionowe dla każdego z kroków czasowych. Każdy obwód będzie składał się z pojedynczej bramki ewolucji, z czasem ewolucji zadeklarowanym wcześniej. Operatorem ewolucji jest hamiltonian. Później uruchomimy na tych obwodach przebiegi transpilera, aby utworzyć obwody qDRIFT.

Przygotowanie ansatzu

Przygotowujemy stan Hartree-Focka za pomocą klasy InitializeModes. Dla azotu proces ten polega po prostu na zastosowaniu bramek X do pierwszych num_elec_a kubitów, a następnie do num_elec_b kubitów, przy czym obie te wartości są równe siedem dla azotu. Ten stan reprezentuje siedem elektronów $\alpha$ i siedem $\beta$ azotu.

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []

hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

Krok 2: Optymalizacja problemu pod kątem wykonania na sprzęcie kwantowym

Teraz, gdy mamy nasze obwody, najpierw użyjemy przebiegów dostępnych w qiskit-fermions, aby przeprowadzić optymalizacje na poziomie fermionowym, a następnie transpilujemy nasz obwód dla wybranego backendu. Ponieważ jest to eksperyment symulacyjny, zrobimy to najpierw dla AerSimulator. Obliczanie wag dla każdej grupy

W tym kroku wykonujemy próbkowanie qDRIFT wyrazów stochastycznie, z prawdopodobieństwami proporcjonalnymi do ich współczynników w hamiltonianie. Przebieg transpilera qDRIFT robi to za nas. Możemy teraz tworzyć płytsze obwody, które można wykonywać na sprzęcie bardziej efektywnie, mimo ograniczonej łączności kubitów, nawet gdy hamiltonian zawiera sprzężenia dalekiego zasięgu i wyrazy wyższego niż kwadratowy rzędu. Po pogrupowaniu wyrazów próbkuje on operatory na podstawie ich wag. Dla każdego operatora $h_i$ waga $W_{h_i}$ jest zdefiniowana następująco:

$$W_{h_i} = |c_i| / \lambda$$

Optymalizacje fermionowe i natywne dla sprzętu

Funkcja generate_preset_jw_pass_manager() zwraca MultiStagePassManager, który przyjmuje FermionicCircuit i tworzy zoptymalizowany obwód końcowy, który możemy transpilować do uruchomienia na naszym sprzęcie. Zastępujemy jej domyślny etap optymalizacji przez FermionicPassManager zawierający nasz przebieg QDriftTrotterization:

• Przebieg QDriftTrotterization wewnętrznie wykorzystuje obliczanie wag i próbkowanie, aby wygenerować obwody, których użyjemy do próbkowania

• Przebieg RelabelModes to kolejny przebieg optymalizacyjny, który można wykorzystać do permutowania trybów fermionowych w celu optymalizacji łączności między kubitami i zmniejszenia głębokości bramek; więcej informacji znajdziesz w dokumentacji API

Pozostałe etapy MultiStagePassManager uruchamiają się automatycznie i obsługują pełne mapowanie fermionów na kubity:

• F2QLayout: Ustawiony menedżer przebiegów stosuje przebieg TrivialF2QLayout, który w trywialny sposób mapuje $n$ bitów fermionowych na $n$ kubitów.

• F2QSynth: Przebieg transpilacji mapujący instrukcje obwodu oparte na fermionach na instrukcje oparte na kubitach.

qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

400

Teraz, gdy zakończyliśmy optymalizacje na poziomie fermionowym, możemy transpilować obwody do wykonania na symulatorze.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Krok 3: Wykonanie za pomocą prymitywów Qiskit

Teraz, gdy mamy nasze obwody, możemy je uruchomić za pomocą prymitywów Qiskit na AerSimulator. Połączymy wszystkie zliczenia z różnych obwodów. Przed ostatecznym przetworzeniem za pomocą SQD konwertujemy je na wektory boolowskie.

print(
    f"Executing {len(transpiled_circuits)} circuits with {shots} shots each..."
)

job = simulator.run(transpiled_circuits, shots=shots)
result = job.result()

all_counts = [result.get_counts(i) for i in range(len(transpiled_circuits))]

print(len(all_counts), "length before post processing")

Executing 400 circuits with 100 shots each...
400 length before post processing

Krok 4: Przetwarzanie końcowe i zwrócenie wyniku w pożądanym formacie klasycznym

Wykorzystanie ciągów bitowych do SQD

Możemy teraz uruchomić schemat diagonalizacji na wybranych ciągach bitowych, aby znaleźć najniższą wartość własną, która będzie odpowiadać energii stanu podstawowego cząsteczki. Tworzymy funkcję zwrotną (callback), deklarujemy początkowe obsadzenia i ustawiamy parametry, zanim w końcu uruchomimy schemat diagonalizacji. Funkcja zwrotna służy do wypisywania bieżącej iteracji i bieżącego oszacowania wartości własnej przy każdej iteracji.

Na koniec, aby uzyskać oszacowanie stanu podstawowego, dodajemy nuclear_repulsion_energy do wynikowej energii.

Uwaga: Wymiar podprzestrzeni nie jest ustalony między iteracjami, nawet na bezszumowym symulatorze — każda podpróbka pobiera inny zestaw konfiguracji, a krok odzyskiwania przekształca pulę między iteracjami, więc zgłaszany wymiar różni się między kolejnymi podpróbkami. Samo próbkowanie bezszumowe nie ustala wymiaru wybranej podprzestrzeni. Uruchomienie na sprzęcie ma jednak tendencję do dawania systematycznie większych podprzestrzeni, ponieważ zaszumione strzały łamią symetrię liczby cząstek, a odzyskiwanie konfiguracji przekształca je w dodatkowe wektory bazowe. Z tego powodu w sekcji dotyczącej sprzętu wprowadzimy również kolejny krok przycinania ciągów bitowych.

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

40000
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64767025226178
		Subspace dimension: 5538
	Subsample 1
		Energy: -107.64772799119115
		Subspace dimension: 5670
	Subsample 2
		Energy: -107.64765512281548
		Subspace dimension: 5767
Iteration 2
	Subsample 0
		Energy: -107.64795948524682
		Subspace dimension: 6080
	Subsample 1
		Energy: -107.64806617355072
		Subspace dimension: 6300
	Subsample 2
		Energy: -107.64802260640258
		Subspace dimension: 6308
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999464 0.99999643 0.99584631 0.99332984 0.96684652 0.96686712
 0.99301927 0.0373282  0.0373266  0.00944508]
Orbital occupancies (beta): [0.99999462 0.99999643 0.9958261  0.99332349 0.96684268 0.96686737
 0.99302145 0.03733536 0.03733399 0.0094585 ]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6480661736 Ha
Error:            1.1811817564e-04 Ha

Przykład na sprzęcie

Ten przykład wykorzystuje 20 kubitów (10 orbitali przestrzennych). Wybór ten jest wygodny dla samouczka, który ma działać szybko, a nie stanowi sztywnego ograniczenia dla metody.

Koszt kroku klasycznego nie jest bezpośrednio określany przez liczbę kubitów. SQD diagonalizuje hamiltonian rzutowany na podprzestrzeń rozpiętą przez próbkowane konfiguracje, więc koszt klasyczny zależy od wymiaru tej wybranej podprzestrzeni — determinowanej tu przez samples_per_batch, num_batches oraz liczbę różnych konfiguracji, które obwody faktycznie generują — a także od algebry liniowej rzadkich macierzy potrzebnej do zastosowania rzutowanego hamiltonianu. Pełna przestrzeń CI rośnie kombinatorycznie wraz z liczbą orbitali i elektronów, ale wybrana podprzestrzeń jest jej małym, regulowanym wycinkiem, którego rozmiar kontrolujemy bezpośrednio. W rezultacie liczbę kubitów i trudność klasyczną można w pewnym stopniu zmieniać niezależnie: szersza przestrzeń orbitali próbkowana do skromnej podprzestrzeni może być tańsza niż mniejszy układ diagonalizowany na bardzo dużej podprzestrzeni.

W praktyce więc wykonalny rozmiar układu zależy od wymiaru podprzestrzeni potrzebnego do osiągnięcia pożądanej dokładności oraz od pamięci i liczby rdzeni dostępnych dla solvera własnego. Większe przestrzenie orbitali zwykle wymagają większej podprzestrzeni, aby osiągnąć dokładność chemiczną, i to właśnie ostatecznie motywuje wykorzystanie zasobów rozproszonych — zobacz qiskit-addon-sqd-hpc, aby dowiedzieć się, jak skalować ten krok. Zamiast zakładać stały próg odcięcia, praktycznym podejściem jest obserwowanie zgłaszanego wymiaru podprzestrzeni oraz zbieżności energii między iteracjami i zwiększanie rozmiaru podprzestrzeni, aż energia przestanie się poprawiać lub wyczerpie się dostępna pamięć.

Uwaga: Ze względu na błąd próbkowania wynikający z szumu na sprzęcie, podprzestrzeń utworzona do diagonalizacji w uruchomieniu na sprzęcie będzie większa niż ta, którą otrzymujemy przy użyciu symulatora. Choć zwiększa to wymiar podprzestrzeni, którą chcemy diagonalizować, przepływ pracy nadal daje nam dokładną odpowiedź dzięki odporności SQD na szum.

Przycinanie fałszywych ciągów

Tutaj możemy wybrać wykonanie dodatkowego kroku. Gdy mamy już wszystkie ciągi bitowe z wykonań obwodów, możemy albo odfiltrować nieprawidłowe ciągi bitowe przed uruchomieniem SQD, albo kontynuować bez przycinania. Pominięcie przycinania jest zazwyczaj preferowane w przypadku uruchomień na sprzęcie, ponieważ pozostawia strzały łamiące symetrię dostępne dla odzyskiwania konfiguracji, które może naprawić je do prawidłowych konfiguracji i tym samym poszerzyć podprzestrzeń, zamiast całkowicie odrzucać te strzały.

Ponieważ azot może mieć tylko siedem elektronów $\alpha$ i siedem $\beta$, wszelkie ciągi bitowe zawierające więcej lub mniej niż siedem jedynek w pierwszej i drugiej połowie wyjścia można odrzucić. Definiujemy funkcję sprawdzającą, czy ciągi bitowe są prawidłowe, a jeśli nie, odrzucającą je. Po odfiltrowaniu fałszywych ciągów bitowych reszta jest przekazywana do schematu diagonalizacji. Użyj flagi PRUNE poniżej, aby przełączać się między tymi dwoma zachowaniami.

Pamiętaj, że przycinanie jest tylko jednym z kilku wyborów kształtujących ostateczną podprzestrzeń, obok liczby obwodów, zestawu czasów ewolucji i filtrowania wyrazów diagonalnych. Porównanie uruchomienia z przycinaniem i bez przycinania jest miarodajne tylko wtedy, gdy wszystko inne pozostaje niezmienne; towarzyszący kod w C++ omawia to bardziej szczegółowo, ponieważ dokonuje postselekcji zamiast odzyskiwania, a także różni się pod względem tych innych parametrów.

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

# Extract metadata from the FCIDump header
norb = fcidump["NORB"]  # number of spatial orbitals
nelec = fcidump["NELEC"]  # total number of electrons
e_nuc = fcidump["ECORE"]  # nuclear repulsion / core energy
ms2 = fcidump["MS2"]  # 2S (spin)

num_elec_a = (nelec + ms2) // 2  # alpha electrons
num_elec_b = (nelec - ms2) // 2  # beta  electrons

# Reconstruct full 4-index ERIs from the FCIDump (stored in 8-fold symmetry)
h1e = fcidump["H1"]  # shape (norb, norb)
h2e = ao2mo.restore(  # shape (norb, norb, norb, norb)
    1, fcidump["H2"], norb
)

cisolver = fci.direct_spin1.FCI()
cisolver.max_cycle = 200
cisolver.conv_tol = 1e-12

e_fci, _ = cisolver.kernel(
    h1e,
    h2e,
    norb,
    (num_elec_a, num_elec_b),
    ecore=e_nuc,  # adds nuclear repulsion to the final energy
)

reference_energy = e_fci

print(f"Reference FCI Energy  = {reference_energy:.10f} Ha")

nuclear_repulsion_energy = fcidump["ECORE"]
print(f"Nuclear Repulsion Energy = {nuclear_repulsion_energy:.10f} Ha")

fcidump = FCIDump.from_file(name)
hamiltonian = FermionOperator.from_fcidump(fcidump)
num_modes = 2 * fcidump.norb

# Apply automatic grouping
canon = canonical_order(hamiltonian.normal_ordered().simplify(atol=1e-16))
exit_code = group_terms_by_electronic_structure(
    canon, num_modes, two_body_physicist_order=False
)
filter_diagonal_terms(canon)

print(len(canon.groups))

# SqDRIFT parameters
times = [1.0, 10.0]  # Total evolution times used for the subspace creation
num_exc = 10  # Number of excitation groups per circuit
num_circuits = 200  # Number of circuits to generate

init_circuits = []
hf_gate = InitializeModes.from_hartree_fock(norb, (num_elec_a, num_elec_b))

for time in times:
    evo_gate = Evolution(num_modes, canon, time)
    circ = FermionicCircuit(num_modes)
    circ.append(hf_gate, circ.modes)
    circ.append(evo_gate, circ.modes)
    init_circuits.append(circ)

# Calculate weights for sampling (one per group)
qdrift = QDriftTrotterization(num_exc, rng=19)

pm = generate_preset_jw_pass_manager()
pm.optimization = FermionicPassManager([qdrift])

sqdrift_circuits = []
for circ in init_circuits:
    sqdrift_circuits += (pm.run(circ) for _ in range(num_circuits))

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

# This example assumes you have saved your IBM Quantum Platform account locally.
service = QiskitRuntimeService(channel="ibm_quantum_platform")

# Select backend (choose based on qubit requirements)
backend = service.least_busy(
    operational=True,
    simulator=False,
    min_num_qubits=2 * norb,
)

print(f"Selected backend: {backend.name} ({backend.num_qubits} qubits)")

# Transpile for hardware
transpiled_circuits = transpile(
    sqdrift_circuits,
    backend=backend,
    optimization_level=3,
    seed_transpiler=42,
)

shots = 100

sampler = Sampler(mode=backend)

sampler.options.environment.job_tags = ["TUT-SqDRIFT"]

job = sampler.run(transpiled_circuits, shots=shots)
result = job.result()

# Extract counts from SamplerV2 results
all_counts = [pub_result.data.meas.get_counts() for pub_result in result]

# Set to True to filter out bitstrings that violate electron-number conservation
PRUNE = False

def is_valid_bitstring(
    bitstring: str, norb: int, nelec: tuple[int, int]
) -> bool:
    n_alpha, n_beta = nelec
    return (
        len(bitstring) == 2 * norb
        and bitstring[norb:].count("1") == n_alpha
        and bitstring[:norb].count("1") == n_beta
    )

if PRUNE:
    all_counts_filtered = []
    for counts in all_counts:
        filtered_count = {}
        for key in counts:
            if not is_valid_bitstring(key, norb, (num_elec_a, num_elec_b)):
                continue
            elif key not in filtered_count.keys():
                filtered_count[key] = counts[key]
            else:
                filtered_count[key] += counts[key]
        all_counts_filtered.append(filtered_count)
    all_counts = all_counts_filtered

combined_counts = {}
for counts in all_counts:
    for bitstring, count in counts.items():
        combined_counts[bitstring] = combined_counts.get(bitstring, 0) + count

bit_array = BitArray.from_counts(combined_counts)
print(bit_array.num_shots)

print("Electron configuration:")
print(f"  Total electrons: {nelec}")
print(f"  Alpha electrons: {num_elec_a}")
print(f"  Beta electrons: {num_elec_b}")
print(f"  Number of orbitals: {norb}")
print(f"  Number of spin orbitals (qubits): {2*norb}")
print(f"Integral shapes: h1e={h1e.shape}, h2e={h2e.shape}")

# SQD parameters
samples_per_batch = 300
num_batches = 3
max_iterations = 5

initial_occupancies = (
    np.array([1] * num_elec_a + [0] * (norb - num_elec_a)),  # alpha
    np.array([1] * num_elec_b + [0] * (norb - num_elec_b)),  # beta
)

result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

# Run SQD with configuration recovery
print("\nRunning SQD with configuration recovery...")
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=(num_elec_a, num_elec_b),
    num_batches=num_batches,
    energy_tol=1e-3,
    occupancies_tol=1e-3,
    max_iterations=max_iterations,
    initial_occupancies=initial_occupancies,
    seed=42,
    callback=callback,
)

computed_energy = result.energy + nuclear_repulsion_energy

print("FINAL SQD RESULTS")
print(f"Orbital occupancies (alpha): {result.orbital_occupancies[0]}")
print(f"Orbital occupancies (beta): {result.orbital_occupancies[1]}")

energy_error = abs(computed_energy - reference_energy)
print(f"Reference Energy: {reference_energy:.10f} Ha")
print(f"Computed Energy:  {computed_energy:.10f} Ha")
print(f"Error:            {energy_error:.10e} Ha")

Parsing assets/sqdrift/fcidump_files/N2_sto_3g
Reference FCI Energy  = -107.6481842917 Ha
Nuclear Repulsion Energy = 23.7887003074 Ha
5060
400

Selected backend: ibm_aachen (156 qubits)
40000
Electron configuration:
  Total electrons: 14
  Alpha electrons: 7
  Beta electrons: 7
  Number of orbitals: 10
  Number of spin orbitals (qubits): 20
Integral shapes: h1e=(10, 10), h2e=(10, 10, 10, 10)

Running SQD with configuration recovery...
Iteration 1
	Subsample 0
		Energy: -107.64593072647523
		Subspace dimension: 7221
	Subsample 1
		Energy: -107.6458270048177
		Subspace dimension: 7209
	Subsample 2
		Energy: -107.64007673117075
		Subspace dimension: 7138
Iteration 2
	Subsample 0
		Energy: -107.64757372124944
		Subspace dimension: 9009
	Subsample 1
		Energy: -107.64674060104392
		Subspace dimension: 8245
	Subsample 2
		Energy: -107.64731360491942
		Subspace dimension: 8178
Iteration 3
	Subsample 0
		Energy: -107.64765518770588
		Subspace dimension: 8835
	Subsample 1
		Energy: -107.64767975712016
		Subspace dimension: 8649
	Subsample 2
		Energy: -107.64761634415606
		Subspace dimension: 8648
FINAL SQD RESULTS
Orbital occupancies (alpha): [0.99999504 0.9999964  0.99590318 0.9932359  0.96697158 0.96696295
 0.99298797 0.03728154 0.03728186 0.00938359]
Orbital occupancies (beta): [0.9999946  0.99999641 0.99590413 0.99323077 0.96697361 0.96696174
 0.99298424 0.03728121 0.03728169 0.00939159]
Reference Energy: -107.6481842917 Ha
Computed Energy:  -107.6476797571 Ha
Error:            5.0453460619e-04 Ha

Kolejne kroki

Zalecenia

Jeśli ta praca cię zainteresowała, mogą cię zainteresować następujące materiały:

• Diagonalizacja kwantowa Kryłowa oparta na próbkowaniu dla fermionowego modelu sieciowego - powiązany samouczek wykorzystujący obwody ewolucji czasowej zamiast ansatzu wariacyjnego.
• Diagonalizacja kwantowa oparta na próbkowaniu dla hamiltonianu chemicznego - samouczek o tym, jak skonstruować obwód lokalnego unitarnego klastrowego Jastrowa (LUCJ) do symulacji chemii kwantowej.
• Artykuł SqDRIFT - literatura, na której oparty jest ten samouczek. (Zauważ, że niektóre optymalizacje omówione w tym artykule są obecnie w trakcie realizacji, a ten samouczek może ulec zmianie w przyszłości wraz z rozwojem wykorzystywanych bibliotek.)