Algorytm SqDRIFT do estymacji stanu podstawowego

Szacowane użycie: 180 sekund na procesorze Heron r3 (UWAGA: to tylko szacunek. Twój czas wykonania może się różnić.)

Cele nauczania

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

• Przejdź przez kompletny przepływ pracy estymacji stanu podstawowego z użyciem qDRIFT i SQD

• Dowiedz się, jak używać qiskit-fermions razem z innymi dodatkami Qiskit do realizacji takiego przepływu pracy

Ten samouczek jest przedstawiony jako notebook Pythona w celach dydaktycznych.

Wymagania wstępne

• Przeczytaj omówienie Sample-based quantum diagonalization (SQD)

• Przeczytaj lekcję Sample-based Krylov Quantum Diagonalization (SKQD)

Tło

SqDRIFT jest wariantem SKQD, który zastępuje konieczność wyboru ansatzu, z którego próbkowane są ciągi bitów, zespołem obwodów ewolucji w czasie zbudowanych bezpośrednio z docelowego hamiltonianu. Osiąga się to przez podpróbkowanie mniejszych operatorów ewolucji w czasie z hamiltonianu na podstawie jego współczynników, co jest znane jako metoda trotteryzacji qDRIFT.

Ten samouczek korzysta z Qiskit Fermions, aby utworzyć bardziej naturalne obwody fermionowe dla algorytmu qDRIFT, a następnie użyć fermionowych przebiegów układu i syntezy, zanim obwody trafią do tradycyjnego potoku Qiskit w celu wykonania na sprzęcie.

Niech hamiltonian ma postać:

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

gdzie, bez utraty ogólności, wymagamy, aby $c_i > 0$ oraz aby największa wartość własna $h_i$ była równa co do modułu $1$. Każdy czynnik ze znakiem lub zespolony jest wchłaniany w $h_i$, więc współczynniki $c_i$ są ściśle dodatnimi wagami, a $h_i$ niosą kierunek każdego członu. Tu $N$ jest liczbą członów (lub, po grupowaniu, liczbą grup) w hamiltonianie; jest to właściwość hamiltonianu, odrębna od liczby operatorów próbkowanych do pojedynczego obwodu, oznaczanej poniżej przez $n$.

Algorytm qDRIFT realizuje wtedy, dla docelowego czasu $t$, pewien operator $V_k$, gdzie $k$ przebiega od $1 \cdots K$ i oznacza $k_{th}$ obwód SqDRIFT, zdefiniowany jako:

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

Tu $n$ jest liczbą próbkowanych operatorów na obwód, a $K$ liczbą obwodów w zespole. Iloczyn przebiega po $n$ losowaniach, a nie po wszystkich $N$ członach hamiltonianu, a ponieważ człony są losowane ze zwracaniem, ten sam $h_i$ może wystąpić więcej niż raz w jednym $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 $\lambda t / n$ niezależnie od tego, który człon został wylosowany. Jednorodność kąta kroku jest cechą charakterystyczną qDRIFT: współczynnik wpływa na wynik przez to, jak często jego człon jest losowany, a nie przez to, jak daleko ten człon jest obrócony. Indeksy są próbkowane z rozkładu:

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

zatem ciąg $(k_1, \ldots, k_n)$ jest losową sekwencją indeksów członów wylosowanych 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 działaniem $H$, z błędem malejącym 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 oznacza liczbę członów jako $\mathcal{N}$, a długość sekwencji jako $N$; tutaj używamy $N$ i $n$, aby wyraźnie je rozróżnić.)

Ten samouczek pokazuje, jak wygenerować zespół takich losowych obwodów. Po ich utworzeniu, podobnie jak tworzymy podprzestrzeń Kryłowa dla różnych operatorów, próbkujemy ciągi bitów z wielu takich operatorów o różnych parametrach czasu. Zapewnia to większe nakładanie się wektorów stanu podstawowego i próbkowanych ciągów bitów.

Wymagania

Zanim zaczniesz ten samouczek, upewnij się, że masz zainstalowane

• Wirtualne środowisko Pythona (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (zauważ, ż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 — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "qiskit_aer": "qiskit-aer", "qiskit_fermions": "qiskit-fermions", "qiskit_ibm_runtime": "qiskit-ibm-runtime"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

# 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 z symulatorem

Krok 1: Odwzorowanie danych klasycznych na problem kwantowy

Wczytanie i przygotowanie FCIDump

W tym samouczku wczytamy hamiltonian struktury elektronowej dla azotu (N2). Istnieją też inne sposoby tworzenia operatorów fermionowych. Zajrzyj do dokumentacji w qiskit_fermions.operators.library.

O tym 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 długości wiązania w równowadze. Jego nagłówek deklaruje NORB=10, NELEC=14 i MS2=0: 10 orbitali przestrzennych (stąd 20 orbitali spinowych i 20 kubitów w kodowaniu Jordana-Wignera), 14 elektronów w singlecie spinowym, czyli siedem elektronów $\alpha$ i siedem $\beta$. Wszystkie orbitale mają etykietę symetrii 1, to znaczy nie wykorzystuje się symetrii grupy punktowej. Ponieważ jest to zrzut pełnej przestrzeni STO-3G, żadne orbitale nie są zamrożone, a przestrzeń korelacji jest na tyle mała, że dokładną energię odniesienia FCI można obliczyć klasycznie do porównania, jak 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 fazą lub kolejnością orbitali; energie całkowite pozostają niezmienione.

Pobranie pliku. FCIDump znajdziesz w tym repozytorium GitHub. Możesz uruchomić poniższą komórkę, aby pobrać go do lokalizacji, której oczekuje reszta samouczka.

Najpierw używamy cisolver z pyscf, aby uzyskać energię odniesienia. 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 użyte także 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 = "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 = "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

Wczytanie 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

Fermionowe przepływy pracy z qiskit-fermions

Najpierw odwzorujemy hamiltonian na fermionowy model obwodów za pomocą qiskit-fermions, który udostępnia przebiegi transpilera i bramki specyficzne dla obwodów fermionowych. Zostaną one później użyte przed tradycyjnymi przebiegami transpilera Qiskit w tym przepływie pracy.

Grupowanie członów

Aby zapewnić powtarzalność wyników, najpierw używamy canonical_order do posortowania członów wyłącznie na podstawie ich struktury. Kolejność operatorów na liście canon jest więc ustalona. Zapewnia to powtarzalność tworzonych operatorów, ponieważ przebieg QDriftTrotterization, którego użyjemy później, próbkuje losowe indeksy w celu utworzenia operatorów qDRIFT.

W tym kroku wykorzystujemy liczne symetrie obecne w hamiltonianie struktury elektronowej, grupując powiązane człony 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 istotne, grupowanie członów powiązanych symetrią skutkuje korzystnym znoszeniem się członów Pauliego i ogólnie mniejszą głębokością obwodu przy ewolucji stanu w czasie pod ich działaniem.

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

Zauważ, że group_terms_by_electronic_structure zakłada człony w porządku normalnym.

Filtrowanie członów diagonalnych

Usuwamy człony diagonalne z hamiltonianu używanego do generowania obwodów, aby $n$ miejsc próbkowania qDRIFT zostało wykorzystanych na człony przenoszące populację między konfiguracjami. Takie człony najlepiej odfiltrować z hamiltonianu w tym miejscu, zanim w następnym kroku zostanie skonstruowana bramka Evolution.

Omawiane człony to te, które są diagonalne w bazie liczb obsadzeń, czyli iloczyny operatorów liczby $a^\dagger_i a_i$. Pod ten opis podpadają trzy rodzaje członów:

• stały przesunięcie energii, iloczyn zerowej liczby operatorów liczby, którego ewolucja w czasie wnosi tylko globalną fazę;

• pojedyncze operatory liczby $n_i$, których ewolucja w czasie sprowadza się do obrotów $Z$ pojedynczych kubitów;

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

Same w sobie żaden z nich nie przenosi populacji między konfiguracjami liczb obsadzeń; działają tylko na fazy konfiguracji już obecnych. Nie są jednak bierne: te względne fazy wpływają na interferencję generowaną przez człony wzbudzeń dalej 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 obwodów, wykonane po to, by skupić próbkowanie na członach wzbudzeń, a nie krok pozostawiający rozkład próbkowania nietknięty. W odróżnieniu od powyższego grupowania według symetrii, które pozostawia gwarancje zbieżności qDRIFT nienaruszone, ten filtr zmienia ewoluowany operator. Obwody nie przybliżają więc już ewolucji pod działaniem pełnego hamiltonianu, a granice błędu qDRIFT dotyczą przefiltrowanego operatora, a nie oryginalnego. Jest to dopuszczalne, ponieważ obwody są tylko heurystyką próbkowania służącą do proponowania konfiguracji: żaden człon nie jest tracony z samego oszacowania energii, gdyż filtr dotyczy tylko hamiltonianu użytego do zbudowania obwodów, podczas gdy późniejsza klasyczna diagonalizacja używa pełnego hamiltonianu, wraz z członami diagonalnymi. Dokładność SQD zależy od tego klasycznego etapu, który pozostaje wariacyjny w spróbkowanej podprzestrzeni niezależnie od tego, jak zaproponowano konfiguracje.

Funkcja filter_diagonal_terms() usuwa takie człony z operatora w miejscu. Rozpoznaje je po ich strukturze w porządku normalnym — multizbiór modów kreacji zgadza się z multizbiorem modów anihilacji — więc jest poprawna tylko dla operatora, który jest już w porządku normalnym. 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

Skoro zgrupowaliśmy człony w hamiltonianie, zdecydujemy o następujących parametrach służących do wygenerowania zespołu obwodów:

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

Tworzenie obwodów fermionowych

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

Przygotowanie ansatzu

Przygotowujemy stan Hartree-Focka za pomocą klasy InitializeModes. Dla azotu proces polega po prostu na zastosowaniu bramek X do pierwszych num_elec_a kubitów, a następnie do num_elec_b kubitów, z których oba 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 do wykonania na sprzęcie kwantowym

Mając nasze obwody, użyjemy najpierw przebiegów dostępnych w qiskit-fermions do optymalizacji na poziomie fermionowym, a następnie dokonamy transpilacji obwodu dla wybranego backendu. Ponieważ jest to eksperyment na symulatorze, najpierw zrobimy to dla AerSimulator. Obliczenie wag dla każdej grupy

W tym kroku wykonujemy stochastyczne próbkowanie qDRIFT składników z prawdopodobieństwami proporcjonalnymi do ich współczynników w hamiltonianie. Robi to za nas przebieg transpilera qDRIFT. Możemy teraz tworzyć płytsze obwody, które da się wydajniej wykonać na sprzęcie mimo ograniczonej łączności kubitów, nawet gdy hamiltonian zawiera sprzężenia dalekiego zasięgu i składniki wyższego rzędu niż kwadratowe. Po zgrupowaniu składnikó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$$

Ponieważ składniki zostały zgrupowane w kroku 1, każde $h_i$ oznacza tutaj całą grupę: $c_i$ jest średnią wartością bezwzględną współczynników składników w grupie $i$, a każdy składnik w grupie ewoluuje ze współczynnikiem zredukowanym do swojego znaku.

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 poddać transpilacji, aby uruchomić go na naszym sprzęcie. Zastępujemy jego domyślny etap optymalizacji obiektem FermionicPassManager zawierającym nasz przebieg QDriftTrotterization:

• Przebieg QDriftTrotterization wewnętrznie wykorzystuje obliczanie wag i próbkowanie do generowania obwodów, których użyjemy do próbkowania

• Przebieg RelabelModes to kolejny przebieg optymalizacyjny, którego można użyć do permutowania modów fermionowych, aby zoptymalizować łączność między kubitami i zmniejszyć głębokość bramek; więcej informacji znajdziesz w dokumentacji API

Pozostałe etapy MultiStagePassManager działają automatycznie i obsługują pełne odwzorowanie fermionów na kubity:

• F2QLayout: Predefiniowany pass manager stosuje przebieg TrivialF2QLayout, który trywialnie odwzorowuje $n$ bitów fermionowych na $n$ kubitów.

• F2QSynth: Przebieg transpilacji odwzorowują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

Skoro mamy już za sobą optymalizacje na poziomie fermionowym, możemy poddać obwody transpilacji w celu wykonania na symulatorze.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Krok 3: Wykonanie za pomocą prymitywów Qiskit

Mając obwody, możemy uruchomić je za pomocą prymitywów Qiskit na AerSimulator. Połączymy wszystkie zliczenia z różnych obwodów. Przekształcamy je w wektory boolowskie, zanim wreszcie poddamy je przetwarzaniu końcowemu za pomocą SQD.

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 żądanym formacie klasycznym

Użycie ciągów bitów w SQD

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

Na koniec, aby otrzymać oszacowanie stanu podstawowego, dodajemy nuclear_repulsion_energy do uzyskanej energii.

Uwaga: Wymiar podprzestrzeni nie jest stały w kolejnych iteracjach, nawet na bezszumowym symulatorze — każda podpróbka losuje inny zestaw konfiguracji, a krok odzyskiwania zmienia pulę między iteracjami, więc raportowany wymiar zmienia się z jednej podpróbki na drugą. Samo bezszumowe próbkowanie nie ustala wymiaru wybranej podprzestrzeni. Uruchomienie na sprzęcie daje jednak zwykle systematycznie większe podprzestrzenie, ponieważ zaszumione strzały łamią symetrię liczby cząstek, a odzyskiwanie konfiguracji zamienia je w dodatkowe wektory bazowe. Z tego powodu w sekcji o sprzęcie wprowadzimy również dodatkowy krok przycinania ciągów bitów.

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). Ten wybór jest wygodny w tutorialu, który ma działać szybko, a nie stanowi twardego ograniczenia metody.

Koszt kroku klasycznego nie jest wyznaczany bezpośrednio przez liczbę kubitów. SQD diagonalizuje hamiltonian zrzutowany na podprzestrzeń rozpiętą przez próbkowane konfiguracje, więc koszt klasyczny zależy od wymiaru tej wybranej podprzestrzeni — który tutaj wyznaczają samples_per_batch, num_batches oraz to, ile różnych konfiguracji faktycznie wytwarzają obwody — a także od rzadkiej algebry liniowej potrzebnej do zastosowania zrzutowanego hamiltonianu. Pełna przestrzeń CI rośnie kombinatorycznie wraz z liczbą orbitali i elektronów, ale wybrana podprzestrzeń jest jej niewielkim, regulowanym wycinkiem, a jej rozmiar kontrolujemy bezpośrednio. W konsekwencji liczbę kubitów i trudność klasyczną można zmieniać w pewnym stopniu niezależnie: szersza przestrzeń orbitalna próbkowana do skromnej podprzestrzeni może być tańsza niż mniejszy układ diagonalizowany w 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 rdzeni dostępnych dla solwera wartości własnych. Większe przestrzenie orbitalne zwykle wymagają większej podprzestrzeni, aby osiągnąć dokładność chemiczną, i to właśnie ostatecznie uzasadnia użycie zasobów rozproszonych — zobacz qiskit-addon-sqd-hpc, aby skalować ten krok. Zamiast zakładać stałe ograniczenie, praktyczne podejście polega na obserwowaniu raportowanego wymiaru podprzestrzeni i zbieżności energii w kolejnych iteracjach oraz zwiększaniu rozmiaru podprzestrzeni, aż energia przestanie się poprawiać lub wyczerpie się dostępna pamięć.

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

Przycinanie fałszywych ciągów

Tutaj możemy zdecydować się na wykonanie dodatkowego kroku. Gdy mamy już wszystkie ciągi bitów z wykonań obwodów, możemy albo odfiltrować nieprawidłowe ciągi bitów przed uruchomieniem SQD, albo przejść dalej bez przycinania. Rezygnacja z przycinania jest zazwyczaj lepsza w uruchomieniach na sprzęcie, ponieważ pozostawia strzały ze złamaną symetrią dostępne dla odzyskiwania konfiguracji, które może je naprawić do postaci prawidłowych konfiguracji, a przez to poszerzyć podprzestrzeń, zamiast po prostu odrzucać te strzały.

Ponieważ azot może mieć tylko siedem elektronów $\alpha$ i siedem elektronów $\beta$, wszystkie ciągi bitów, które mają więcej lub mniej niż siedem jedynek w pierwszej i drugiej połowie wyniku, można odrzucić. Definiujemy funkcję, która sprawdza, czy ciągi bitów są prawidłowe, a jeśli nie, odrzuca je. Po odfiltrowaniu fałszywych ciągów bitów pozostałe są przekazywane do schematu diagonalizacji. Użyj poniższej flagi PRUNE, aby przełączać się między tymi dwoma zachowaniami.

Pamiętaj, że przycinanie to tylko jeden z kilku wyborów kształtujących końcową podprzestrzeń, obok liczby obwodów, zbioru czasów ewolucji i filtrowania składników diagonalnych. Porównanie uruchomienia z przycinaniem z uruchomieniem bez przycinania jest informatywne tylko wtedy, gdy wszystko inne pozostaje niezmienione; wersja C++ tego tutorialu omawia to dokładniej, ponieważ stosuje postselekcję zamiast odzyskiwania i różni się również pozostałymi parametrami.

name = "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

Dalsze kroki

Zalecenia

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

• Sample-based Krylov quantum diagonalization of a fermionic lattice model - powiązany tutorial wykorzystujący obwody ewolucji w czasie zamiast ansatzu wariacyjnego.
• Sample-based quantum diagonalization of a chemistry Hamiltonian - tutorial o tym, jak skonstruować obwód lokalnego unitarnego klastra Jastrowa (LUCJ) do symulacji chemii kwantowej.
• Artykuł SqDRIFT - literatura, na której opiera się ten tutorial. (Zwróć uwagę, że niektóre optymalizacje omawiane w tym artykule są obecnie w trakcie prac, a ten tutorial może się zmienić w przyszłości w zależności od rozwoju używanych bibliotek.)