Algorytm SqDRIFT do oszacowania stanu podstawowego

Python

Pobierz notatnik Pythona

Pobierz notatnik Pythona, aby samodzielnie uruchomić ten samouczek. Notatnik pobiera dane wejściowe FCIDump, jeśli nie są jeszcze obecne.

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

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

Wymagania wstępne

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

• Przeczytaj lekcję Kwantowa diagonalizacja Kryłowa oparta 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.

Obie implementacje używają Qiskit Fermions do reprezentowania hamiltonianu molekularnego i konstruowania obwodów dla algorytmu qDRIFT. Implementacja w Pythonie używa fermionowych przebiegów transpilatora, natomiast implementacja w C++ jawnie próbkuje grupy i odwzorowuje je na operatory kubitowe.

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 czynnik ze znakiem (np. ujemny) lub zespolony zostaje włączony do $h_i$, więc współczynniki $c_i$ są ściśle dodatnimi wagami, podczas gdy $h_i$ niosą kierunek każdego wyrazu. Tutaj $N$ to liczba wyrazó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 wyraz został wylosowany. Jednorodność kąta kroku jest charakterystyczną cechą qDRIFT: współczynnik wpływa na wynik poprzez to, jak często jego wyraz jest losowany, a nie poprzez to, jak daleko ten wyraz 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 wyrazó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 wyborach 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 wyrazó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

# 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,
)

Implementacja

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

Poniższe kroki w Pythonie używają symulatora Aer. Przykład sprzętowy następuje później.

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ępnym bloku kodu.

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ższy kod, 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 = "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}")

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

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 modów kreacji odpowiadającego multizbiorowi modó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$$

Ponieważ wyrazy zostały zgrupowane w kroku 1, każde $h_i$ jest tu całą grupą: $c_i$ to średni współczynnik bezwzględny wyrazów w grupie $i$, a każdy wyraz w grupie jest ewoluowany z współczynnikiem zredukowanym do jego 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 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 modó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: Predefiniowany 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; zakładka 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 = "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")

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

C++

Pobierz kod C++

Kod źródłowy C++, projekt CMake i instrukcje budowania znajdują się w katalogu docs/tutorials/assets/sqdrift/cpp repozytorium dokumentacji Qiskit. Dane wejściowe FCIDump są w sąsiednim katalogu fcidump_files. Sklonuj repozytorium lub pobierz oba katalogi, aby samodzielnie zbudować i uruchomić ten samouczek.

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

Efekty uczenia się (C++)

• 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

Wymagania wstępne (C++)

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

• Przeczytaj lekcję Kwantowa diagonalizacja Kryłowa oparta na próbkowaniu (SKQD)

Kontekst (C++)

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.

Obie implementacje używają Qiskit Fermions do reprezentowania hamiltonianu molekularnego i konstruowania obwodów dla algorytmu qDRIFT. Implementacja w Pythonie używa fermionowych przebiegów transpilatora, natomiast implementacja w C++ jawnie próbkuje grupy i odwzorowuje je na operatory kubitowe.

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 czynnik ze znakiem (np. ujemny) lub zespolony zostaje włączony do $h_i$, więc współczynniki $c_i$ są ściśle dodatnimi wagami, podczas gdy $h_i$ niosą kierunek każdego wyrazu. Tutaj $N$ to liczba wyrazó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 wyraz został wylosowany. Jednorodność kąta kroku jest charakterystyczną cechą qDRIFT: współczynnik wpływa na wynik poprzez to, jak często jego wyraz jest losowany, a nie poprzez to, jak daleko ten wyraz 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 wyrazó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 wyborach 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 wyrazó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 (C++)

Ta implementacja działa na sprzęcie IBM Quantum® i używa zewnętrznego pliku wykonywalnego SBD do klasycznej diagonalizacji. Nie zawiera przykładu z symulatorem w Pythonie.

Postępuj zgodnie z instrukcją budowania C++, aby zainstalować kompilator, CMake, Rust, Boost i zależności Qiskit. Zainstaluj Rust za pomocą rustup: zależności Qiskit przypinają własną wersję Rusta, którą rustup pobiera automatycznie, ale Rust zainstalowany przez Homebrew lub menedżer pakietów Linuksa może być zbyt stary. Końcowy krok SBD wymaga dodatkowo MPI, OpenMP oraz BLAS/LAPACK. Windows nie jest obecnie obsługiwany przez tę implementację.

Użyj kompletnego kodu źródłowego C++ i projektu CMake z repozytorium. Zachowaj sąsiedni katalog fcidump_files podczas pobierania przykładu. Z kopii repozytorium dokumentacji zbuduj projekt następującymi poleceniami. Pierwsza kompilacja pobiera i kompiluje wszystkie zależności i trwa kilka minut.

cd docs/tutorials/assets/sqdrift/cpp
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --config Release

Kompilacja zapisuje plik wykonywalny SqDRIFT w docs/tutorials/assets/sqdrift/cpp (katalog główny projektu), a nie w build/. Uruchom ./SqDRIFT z katalogu głównego projektu, aby mógł znaleźć ../fcidump_files/N2_sto_3g. Wszystkie pozostałe polecenia w tej zakładce także uruchamiaj z katalogu głównego projektu, chyba że podano inaczej. Jeśli ./SqDRIFT zakończy się błędem brakującej biblioteki współdzielonej, postępuj zgodnie z instrukcjami dotyczącymi ścieżki bibliotek dla danej platformy w przewodniku budowania; w macOS plik wykonywalny zwykle działa bez nich.

Konfiguracja (C++)

Poniższe fragmenty to wyimki z kompletnego programu SqDRIFT.cpp. Zbuduj i uruchom ten plik jako pojedynczy program.

Przed uruchomieniem SqDRIFT upewnij się, że twój token API IBM Quantum jest skonfigurowany dla IBM Cloud Platform.

Jeśli masz Pythona z zainstalowanym qiskit-ibm-runtime, najłatwiejszy sposób utworzenia poprawnie sformatowanego pliku konfiguracyjnego to:

from qiskit_ibm_runtime import QiskitRuntimeService
QiskitRuntimeService.save_account(channel="ibm_quantum_platform", token="YOUR_IBM_QUANTUM_TOKEN_HERE")

To zapisuje ~/.qiskit/qiskit-ibm.json (Linux/macOS) lub %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) w odpowiednim kształcie. Jeśli wolisz utworzyć plik ręcznie:

{
  "default-ibm-quantum-platform": {
    "channel": "ibm_quantum_platform",
    "token": "YOUR_IBM_QUANTUM_TOKEN_HERE"
  }
}

Uzyskaj swój klucz API z IBM Quantum Platform.

Uwaga: Brakujący lub nieprawidłowy plik konfiguracyjny powoduje przy uruchomieniu mało pomocną panikę Rusta unwrap() zamiast czytelnego błędu. Jeśli napotkasz problem, sprawdź dokładnie nazwę klucza (default-ibm-quantum-platform) i wartość kanału (ibm_quantum_platform).

Implementacja (C++)

Krok 1: Zmapuj dane wejściowe klasyczne na problem kwantowy (C++)

Wczytaj hamiltonian molekularny

char filename[] = "../fcidump_files/N2_sto_3g";
QfFCIDump* fcidump = qf_fcidump_from_file(filename);
if (fcidump == nullptr) {
    std::cerr << "Failed to load FCIDump!" << std::endl;
    return 1;
}

uint32_t norb = qf_fcidump_norb(fcidump);
uint32_t nelec = qf_fcidump_nelec(fcidump);
uint32_t n_alpha = nelec / 2;
uint32_t n_beta  = nelec - n_alpha;
uint32_t num_modes = 2 * norb;
uint32_t num_qubits = num_modes;
std::cout << "✓ Loaded N2 molecule (" << norb << " orbitals, "
          << num_qubits << " qubits, " << nelec << " electrons)" << std::endl;

Output:

✓ Loaded N2 molecule (10 orbitals, 20 qubits, 14 electrons)

Wczytuje hamiltonian cząsteczki $N_2$ z pliku FCIDump. Liczba kubitów jest równa podwojonej liczbie orbitali (spin w górę i spin w dół).

Utwórz hamiltonian i uporządkuj go normalnie

// Create and normal order Hamiltonian
QfFermionOperator* hamiltonian = qf_ferm_op_from_fcidump(fcidump);
QfFermionOperator* normal_ordered = qf_ferm_op_normal_ordered(hamiltonian, nullptr);

Konwertuje dane FCIDump na operator fermionowy i stosuje uporządkowanie normalne, aby uprościć strukturę operatora.

Pogrupuj wyrazy według struktury elektronowej

// Group terms by electronic structure
QfExitCode exit_code = qf_ferm_op_group_terms_by_electronic_structure(
    normal_ordered, num_modes, false
);
if (exit_code != QfExitCode_Success) {
    std::cerr << "Failed to group terms!" << std::endl;
    return 1;
}

uint32_t num_groups = qf_ferm_op_num_groups(normal_ordered);
std::cout << "Grouped into " << num_groups << " groups" << std::endl;

// Split into group operators
QfFermionOperator** group_ops = new QfFermionOperator*[num_groups];
qf_ferm_op_split_out_groups(normal_ordered, nullptr, 0, group_ops);

Output:

Grouped into 1902 groups

Grupuje powiązane wyrazy hamiltonianu według struktury elektronowej przed skonstruowaniem ich obwodów ewolucji.

Oblicz wagi próbkowania, a następnie znormalizuj i odwzoruj na operatory kubitowe

Kolejność ma tu znaczenie. Normalizacja w miejscu przeskalowuje każdy współczynnik do jednostkowej wielkości, więc wagi próbkowania trzeba obliczyć wcześniej. Jeśli normalizacja zadziała pierwsza, każde $|c_j|$ odczytuje się jako dokładnie 1.0, a waga każdej grupy zwija się do zwykłej liczby jej wyrazów — rozkład próbkowania byłby wtedy napędzany rozmiarem grupy, a nie fizycznymi współczynnikami.

Waga każdej grupy to średni współczynnik bezwzględny jej wyrazów. Przebieg QDriftTrotterization w zakładce Python używa tej samej wagi dla zgrupowanego hamiltonianu. Użycie sumy sprawiłoby, że $\lambda$ skalowałoby się z rozmiarem grupy, co z kolei zaburzałoby czas ewolucji $t_k = \lambda\tau/n$.

// Calculate sampling weights BEFORE normalization.
std::vector<double> weights(num_groups);
double total_weight = 0.0;

for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    double abs_sum = 0.0;
    for (uint64_t j = 0; j < num_terms; j++) {
        abs_sum += std::sqrt(coeffs[j].re * coeffs[j].re +
                             coeffs[j].im * coeffs[j].im);
    }
    weights[i] = (num_terms > 0) ? (abs_sum / static_cast<double>(num_terms)) : 0.0;
    total_weight += weights[i];
}

std::cout << "Total weight (λ): " << total_weight << std::endl;

// Normalize the evolution operators term-by-term, now that the sampling
// weights above have already captured the original magnitudes.
for (uint32_t i = 0; i < num_groups; i++) {
    QkComplex64* coeffs;
    uint64_t num_terms;
    qf_ferm_op_get_coeffs(group_ops[i], &coeffs, &num_terms);

    for (uint64_t j = 0; j < num_terms; j++) {
        const double magnitude = std::sqrt(coeffs[j].re * coeffs[j].re +
                                           coeffs[j].im * coeffs[j].im);
        if (magnitude > 0.0) {
            coeffs[j].re /= magnitude;
            coeffs[j].im /= magnitude;
        }
    }
}

QkObs** qubit_ops = new QkObs*[num_groups];
for (uint32_t i = 0; i < num_groups; i++) {
    QfExitCode jw_exit = qf_ferm_op_jordan_wigner(group_ops[i], num_qubits, &qubit_ops[i]);
    if (jw_exit != QfExitCode_Success) {
        std::cerr << "Failed to map group " << i << std::endl;
        return 1;
    }
}
std::cout << "Mapped all " << num_groups << " normalized groups to qubit operators" << std::endl;

Output:

Total weight (λ): 337.261
Mapped all 1902 normalized groups to qubit operators

Waga $\lambda_i = \frac{1}{|G_i|}\sum_{j \in G_i} |c_{ij}|$ jest średnim współczynnikiem bezwzględnym wyrazów grupy $G_i$, odczytanym z oryginalnych współczynników. Normalizacja przeskalowuje następnie każdy wyraz fermionowy do jednostkowej wielkości przed transformacją Jordana-Wignera, zachowując w ewoluowanym operatorze tylko jego fazę/znak. Odpowiada to znormalizowanej ewolucji wyrazów w qDRIFT: wielkości współczynników określają prawdopodobieństwa próbkowania poprzez $\lambda$, a ewolucja obwodu używa znormalizowanego wyrazu, więc duże współczynniki nie są liczone dwukrotnie.

Spróbkuj zestawy operatorów

// SqDRIFT Sampling: Create operator sets
    const int num_circuits = 100;    // Number of circuits to create
    const int ops_per_circuit = 10;  // Operators per circuit
    const double time_step = 1;    // Time step for evolution
    
    std::cout << "\n Generating " << num_circuits << " circuits with " 
              << ops_per_circuit << " operators each..." << std::endl;
    
    std::mt19937 gen(42);
    std::discrete_distribution<> dist(weights.begin(), weights.end());
    
    // Store sampled operator indices
    std::vector<std::vector<int>> operator_sets(num_circuits);
    for (int i = 0; i < num_circuits; i++) {
        operator_sets[i].resize(ops_per_circuit);
        for (int j = 0; j < ops_per_circuit; j++) {
            operator_sets[i][j] = dist(gen);
        }
    }

Output:

 Generating 100 circuits with 10 operators each...

Losuje zestawy operatorów dla 100 losowych obwodów, z których każdy zawiera 10 grup spróbkowanych z rozkładu ważonego. Następny krok zamienia te zestawy w zespół obwodów qDRIFT.

Krok 2: Zbuduj obwody qDRIFT (C++)

W wersji C++ ten krok buduje obwody ze spróbkowanych zestawów operatorów. Transpilacja dla docelowego Backendu odbywa się w kroku 3, tuż przed przesłaniem każdego obwodu.

    // Create Suzuki-Trotter circuits by composing all operators
    std::vector<Qiskit::circuit::QuantumCircuit> circuits;
    circuits.reserve(num_circuits);

    const uint32_t trotter_order = 1;
    const uint32_t trotter_reps = 1;
    const bool preserve_order = false;
    const bool insert_barriers = false;

    // Gate name mapping is static — compute once before all loops
    auto name_map = Qiskit::circuit::get_standard_gate_name_mapping();

    for (int circuit_idx = 0; circuit_idx < num_circuits; circuit_idx++) {
        // Create initial circuit with state preparation using C++ API
        Qiskit::circuit::QuantumCircuit qc(num_qubits, num_qubits);

        // State preparation: apply X gates to spin-up qubits 0..n_alpha-1
        // and spin-down qubits norb..norb+n_beta-1 (derived from FCIDump NELEC)
        for (uint32_t q = 0; q < n_alpha && q < num_qubits; q++) {
            qc.x(q);
        }
        for (uint32_t q = norb; q < norb + n_beta && q < num_qubits; q++) {
            qc.x(q);
        }

        // For each operator in this circuit's set
        for (int op_idx = 0; op_idx < ops_per_circuit; op_idx++) {
            int group_idx = operator_sets[circuit_idx][op_idx];
            double evolution_time = (total_weight * time_step) / ops_per_circuit;
            
            // Create Suzuki-Trotter evolution circuit
            QkCircuit* evolution_circuit = qk_circuit_library_suzuki_trotter(
                qubit_ops[group_idx], trotter_order, trotter_reps,
                evolution_time, preserve_order, insert_barriers
            );
            
            if (evolution_circuit != nullptr) {
                // Get number of instructions
                uint32_t num_ops = qk_circuit_num_instructions(evolution_circuit);
                
                // Manually append each instruction
                for (uint32_t i = 0; i < num_ops; i++) {
                    QkCircuitInstruction inst;
                    qk_circuit_get_instruction(evolution_circuit, i, &inst);

                    // Prepare qubit vector
                    std::vector<uint32_t> qubits(inst.num_qubits);
                    for (uint32_t j = 0; j < inst.num_qubits; j++) {
                        qubits[j] = inst.qubits[j];
                    }

                    // Get operation kind
                    QkOperationKind kind = qk_circuit_instruction_kind(evolution_circuit, i);

                    // Get mutable pointer to circuit
                    QkCircuit* mutable_circuit = qc.get_rust_circuit().get();

                    // Append based on operation type
                    if (kind == QkOperationKind_Gate) {
                        // Convert gate name string to QkGate enum
                        std::string gate_name(inst.name);
                        auto gate_it = name_map.find(gate_name);
                        if (gate_it != name_map.end()) {
                            QkGate gate_enum = gate_it->second.gate_map();
                            qk_circuit_parameterized_gate(
                                mutable_circuit,
                                gate_enum,           // Use QkGate enum, not string
                                qubits.data(),
                                inst.params
                            );
                        } else {
                            std::cerr << "Unknown gate in evolution circuit: " << gate_name << std::endl;
                            qk_circuit_instruction_clear(&inst);
                            qk_circuit_free(evolution_circuit);
                            return 1;
                        }
                    } else if (kind == QkOperationKind_Barrier) {
                        qk_circuit_barrier(
                            mutable_circuit,
                            qubits.data(), 
                            inst.num_qubits
                        );
                    } else if (kind == QkOperationKind_Reset) {
                        qk_circuit_reset(
                            mutable_circuit,
                            qubits[0]
                        );
                    }

                    qk_circuit_instruction_clear(&inst);
                }

                qk_circuit_free(evolution_circuit);
            }
        }

        circuits.push_back(qc);

        // Display progress
        if ((circuit_idx + 1) % 10 == 0) {
            std::cout << "  Created " << (circuit_idx + 1) << "/" 
                    << num_circuits << " circuits" << std::endl;
        }
    }

    std::cout << "✓ Created all " << num_circuits << " Suzuki-Trotter circuits" << std::endl;

Output:

  Created 10/100 circuits
  Created 20/100 circuits
  Created 30/100 circuits
  Created 40/100 circuits
  Created 50/100 circuits
  Created 60/100 circuits
  Created 70/100 circuits
  Created 80/100 circuits
  Created 90/100 circuits
  Created 100/100 circuits
✓ Created all 100 Suzuki-Trotter circuits

Kompletny program następnie wypisuje statystyki obwodów i indeksy operatorów obwodu 0, które tutaj pominięto.

Dla każdego spróbkowanego operatora tworzy obwód ewolucji Suzukiego-Trottera z czasem $t_k = \frac{\lambda \tau}{n}$, gdzie $n$ to liczba operatorów spróbkowanych na obwód, a $\tau$ to krok czasowy. Ponieważ każda spróbkowana grupa została znormalizowana wyraz po wyrazie przed odwzorowaniem, ta ewolucja stosuje podczas rotacji tylko fazę/znak każdego wyrazu; wielkości wnoszą wkład wyłącznie przez rozkład próbkowania. Stan początkowy Hartree-Focka jest przygotowywany przez zastosowanie bramek X do kubitów 0-6 i 10-16 (siedem elektronów w każdym sektorze spinowym). Instrukcje każdego obwodu ewolucji są ręcznie dołączane do obwodu głównego za pomocą API C++.

Krok 3: Wykonanie za pomocą prymitywów Qiskit (C++)

std::cout << "\n Connecting to IBM Quantum Cloud..." << std::endl;

// Initialize IBM Quantum Compute Service using qiskit-cpp
Qiskit::service::QiskitRuntimeService service;
std::cout << "✓ Connected to IBM Quantum" << std::endl;

// Get backend
const std::string backend_name = "ibm_fez";
auto backend = service.backend(backend_name);

if (backend.name().empty()) {
    std::cerr << "Backend " << backend_name
              << " is not available to this account." << std::endl;
    return 1;
}

std::cout << "✓ Selected backend: " << backend.name() << std::endl;

// Create sampler primitive
const int32_t shots = 100;  // Number of shots per circuit
Qiskit::primitives::BackendSamplerV2 sampler(backend, shots);

std::cout << "\n Submitting " << num_circuits << " circuits ("
            << shots << " shots each)..." << std::endl;

// Submit one job per circuit.
std::vector<std::shared_ptr<Qiskit::primitives::BasePrimitiveJob>> jobs;

for (int i = 0; i < num_circuits; i++) {
    for (uint32_t q = 0; q < num_qubits; q++) {
        circuits[i].measure(q, q);
    }
    Qiskit::circuit::QuantumCircuit transpiled_qc =
        Qiskit::compiler::transpile(circuits[i], backend, 2, 1.0, 42);
    std::vector<Qiskit::primitives::SamplerPub> pubs;
    pubs.push_back(Qiskit::primitives::SamplerPub(transpiled_qc, shots));
    jobs.push_back(sampler.run(pubs));
}

std::cout << "\n⏳ Waiting for " << num_circuits << " jobs to complete..." << std::endl;

// Collect results — poll each job and extract bitstrings
std::vector<boost::dynamic_bitset<>> all_bitstrings;

for (int i = 0; i < num_circuits; i++) {
    while (!jobs[i]->in_final_state()) {
        sleep(5);
    }
    if (jobs[i]->done()) {
        auto result = jobs[i]->result();
        if (result.size() > 0) {
            auto& pub_result = result[0];
            auto& bit_array  = pub_result.data();
            for (size_t s = 0; s < bit_array.num_shots(); s++) {
                auto sample = bit_array[s].to_string();
                boost::dynamic_bitset<> bs(num_qubits);
                for (uint32_t bit = 0; bit < num_qubits; bit++) {
                    bs[bit] = (sample[num_qubits - 1 - bit] == '1');
                }
                all_bitstrings.push_back(bs);
            }
        }
    } else {
        std::cout << "  Circuit " << (i+1) << " did not complete successfully" << std::endl;
    }
}

std::cout << "\n✓ Collected " << all_bitstrings.size()
            << " total bitstrings from all circuits" << std::endl;

Łączy się z IBM Quantum, wybiera Backend po nazwie, transpiluje każdy obwód dla tego Backendu i przesyła jedno zadanie na obwód: 100 oddzielnych zadań Sampler po 100 strzałów na rejestrze 20 kubitów. Następnie odpytuje każde zadanie do jego zakończenia i zbiera ciągi bitów z result[0]. Każde zadanie zużywa czas QPU i oddzielnie czeka w kolejce, więc spodziewaj się, że ten krok potrwa znacznie dłużej niż pojedyncze zadanie wsadowe.

Wyniki sprzętowe przytoczone na tej stronie pochodzą z uruchomienia na ibm_pittsburgh (100 obwodów × 100 strzałów), z backend_name zmienionym z domyślnego ibm_fez. Zapisane próbki przetworzono z poprawionym porządkiem bitów determinant pokazanym w kroku 4. Szum sprzętowy zmienia się między uruchomieniami, więc twoje liczności i końcowa energia będą się różnić.

Output:

 Connecting to IBM Quantum Cloud...
✓ Connected to IBM Quantum
✓ Selected backend: ibm_pittsburgh

 Submitting 100 circuits (100 shots each)...

⏳ Waiting for 100 jobs to complete...

✓ Collected 10000 total bitstrings from all circuits

Zmiana Backendu: Zmień backend_name w SqDRIFT.cpp, aby wskazać inne urządzenie; service.backends() wypisuje to, do czego twoje konto ma dostęp.

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

Postselekcja ciągów bitów i zapis ciągów CI

// Postselect bitstrings by Hamming weight (derived from FCIDump NELEC)
std::vector<double> bitstring_weights(all_bitstrings.size(), 1.0);
auto [filtered_bitstrings, filtered_weights] = Qiskit::addon::sqd::postselect_bitstrings(
    all_bitstrings,
    bitstring_weights,
    Qiskit::addon::sqd::MatchesRightLeftHamming<uint32_t>(n_alpha, n_beta)
);
std::ignore = filtered_weights;

std::cout << "✓ Postselected " << filtered_bitstrings.size()
          << " bitstrings with Hamming weight ("
          << n_alpha << "," << n_beta << ")" << std::endl;

// Convert bitstrings to CI strings using SQD addon
auto ci_strings = Qiskit::addon::sqd::bitstrings_to_ci_strings_symmetrize_spin(
    filtered_bitstrings,
    std::nullopt  // No dimension limit
);

std::cout << "✓ Generated " << ci_strings.size() << " CI strings" << std::endl;

// Write CI strings to file for SBD
std::ofstream alpha_file("alphadets_from_sqd.txt");
for (const auto& ci_string : ci_strings) {
    std::string bitstr;
    // SBD expects the highest orbital index first (bit zero on the right).
    boost::to_string(ci_string, bitstr);
    alpha_file << bitstr << "\n";
}
alpha_file.close();

std::cout << "✓ Wrote " << ci_strings.size()
          << " CI strings to alphadets_from_sqd.txt" << std::endl;
std::cout << " Ready for SBD diagonalization!" << std::endl;

Output:

✓ Postselected 1137 bitstrings with Hamming weight (7,7)
✓ Generated 49 CI strings
✓ Wrote 49 CI strings to alphadets_from_sqd.txt
 Ready for SBD diagonalization!

Dokonuje postselekcji ciągów bitów o właściwej wadze Hamminga (n_alpha elektronów ze spinem w górę, n_beta ze spinem w dół, wyprowadzonych z pola NELEC pliku FCIDump), konwertuje je na ciągi Configuration Interaction (CI) z użyciem symetryzacji spinowej i zapisuje je do alphadets_from_sqd.txt dla późniejszej Selected Basis Diagonalization.

Zwolnij zasoby

// Cleanup
// circuits vector will be automatically cleaned up (RAII)
// Free qubit operators
for (uint32_t i = 0; i < num_groups; i++) {
    qk_obs_free(qubit_ops[i]);
}
delete[] qubit_ops;
// Free fermionic operators
for (uint32_t i = 0; i < num_groups; i++) {
    qf_ferm_op_free(group_ops[i]);
}
delete[] group_ops;
qf_ferm_op_free(normal_ordered);
qf_ferm_op_free(hamiltonian);
qf_fcidump_free(fcidump);
return 0;

Prawidłowo zwalnia całą zaalokowaną pamięć i zasoby. std::vector<QuantumCircuit> w C++ jest czyszczony automatycznie dzięki RAII. Ręcznie zwalnia operatory kubitowe, operatory fermionowe i dane FCIDump.

Skoro mamy już bazę, na którą rzutujemy nasz hamiltonian, możemy przejść do procesu diagonalizacji, aby uzyskać oszacowanie stanu podstawowego. Uruchamiaj poniższe polecenia z katalogu głównego projektu, docs/tutorials/assets/sqdrift/cpp.

Zbuduj plik binarny SBD diag

CMake pobiera SBD wyłącznie jako kod źródłowy do deps/sbd/ (główna kompilacja celowo pomija jego budowanie). Aplikacja SBD używa zwykłego Makefile, a sam SBD wymaga MPI, OpenMP oraz BLAS/LAPACK — więc zwykłe make nie wystarczy, dopóki ich nie zainstalujesz, a plik Configuration aplikacji nie będzie odpowiadał twojemu zestawowi narzędzi.

Makefile odczytuje flagi kompilatora i konsolidacji z pliku Configuration w tym samym katalogu. Wersja dostarczana w upstreamie jest dostrojona do jednego konkretnego układu macOS/Homebrew i nie zbuduje się bez zmian na typowej maszynie — w szczególności -I/opt/homebrew/opt/llvm/include nie zawiera omp.h (znajduje się w katalogu zasobów clanga specyficznym dla wersji), a clang++ firmy Apple odrzuca -fopenmp wprost:

clang++: error: unsupported option '-fopenmp'

Poniższy przepis jest zweryfikowany jako działający w macOS (Apple Silicon, Darwin 25.6, Open MPI 5.x + Homebrew LLVM).

macOS (Apple Silicon) — przetestowano

# 1. Install MPI and an OpenMP-capable compiler.
#    Apple clang cannot do OpenMP, so Homebrew's LLVM provides both
#    clang++ and libomp. BLAS/LAPACK come from Apple's Accelerate
#    framework, which is already part of macOS — no openblas needed.
brew install open-mpi llvm

Następnie zastąp deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization/Configuration poniższą zawartością:

# Path to the SBD library
SBD_PATH=../..

# MPI C++ compiler. OMPI_CXX points Open MPI's wrapper at Homebrew's clang++,
# which (unlike Apple clang) supports -fopenmp.
CCCOM=OMPI_CXX=/opt/homebrew/opt/llvm/bin/clang++ mpicxx

# Build flags. omp.h lives in a version-specific clang resource dir, so let
# Homebrew's clang++ find it itself rather than hardcoding an -I path.
CCFLAGS= -std=c++17 -fopenmp -O3

# Link flags: OpenMP runtime from Homebrew LLVM, BLAS/LAPACK from the
# Accelerate framework that ships with macOS (no openblas install needed).
SYSLIB= -L/opt/homebrew/opt/llvm/lib -lomp -framework Accelerate

# 2. Build
cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization
make
cd ../../../..

Pomyślna kompilacja wypisuje dwie linie kompilacji/konsolidacji i zostawia plik binarny diag obok Makefile. Jedno ostrzeżenie ld: warning: ignoring duplicate libraries: '-lomp' jest nieszkodliwe.

Intel macOS: Prefiks Homebrew to /usr/local, a nie /opt/homebrew, więc podstaw go w obu poprzednich ścieżkach.

Linux (Ubuntu/Debian) — nieprzetestowano

Poniższe wymagania wstępne obejmują MPI, OpenMP oraz BLAS/LAPACK. GCC natywnie obsługuje -fopenmp, więc nie trzeba zastępować kompilatora. Nie zostało to zweryfikowane w Linuksie:

sudo apt install -y libopenmpi-dev libomp-dev libblas-dev liblapack-dev

SBD_PATH=../..
CCCOM=mpicxx
CCFLAGS= -std=c++17 -fopenmp -O3
SYSLIB= -llapack -lblas

Uruchom SqDRIFT, aby wygenerować ciągi CI

Jeśli już uruchomiłeś ./SqDRIFT w krokach 1–4, alphadets_from_sqd.txt już istnieje w katalogu głównym projektu i możesz pominąć ten krok. W przeciwnym razie uruchom go z katalogu głównego projektu:

./SqDRIFT
# Writes alphadets_from_sqd.txt in the project root

Uruchom diagonalizację

cd deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization

# The FCIDump lives in the shared sqdrift/fcidump_files/ directory, one level
# above this tutorial's project root, hence five levels up from here.
ln -sf ../../../../../fcidump_files/N2_sto_3g fcidump.txt

# alphadets_from_sqd.txt is written by ./SqDRIFT into the project root itself,
# which is four levels up.
ln -sf ../../../../alphadets_from_sqd.txt alphadets.txt

./diag \
--fcidump fcidump.txt \
--adetfile alphadets.txt \
--method 0 \
--iteration 100 \
--block 10 \
--tolerance 1e-8 \
| tee ../../../../sbd_output.txt

cd ../../../..

Parametry:

• --fcidump: ścieżka do pliku FCIDump zawierającego hamiltonian
• --adetfile: ścieżka do pliku determinant alfa (ciągi CI z SqDRIFT)
• --method 0: metoda diagonalizacji Davidsona
• --iteration 100: maksymalna liczba iteracji Davidsona
• --block 10: rozmiar bloku dla algorytmu Davidsona
• --tolerance 1e-8: tolerancja zbieżności energii

 Elapsed time for helper construction 0.001512 (sec)
 Elapsed time for init 1e-06 (sec)
 Davidson iteration 0.0 (tol=0.1857493937033264): -107.493531425221
 Davidson iteration 0.1 (tol=0.0455699976495391): -107.5098150493379 -105.93368397386
 Davidson iteration 0.2 (tol=0.009291923190816917): -107.5108356028703 -106.1245977311231 -105.6823366883346
 Davidson iteration 0.3 (tol=0.001488674426419356): -107.5108711084406 -106.3219246016818 -105.7376912083761 -105.4392974334491
 Davidson iteration 0.4 (tol=0.0003188697654737439): -107.5108722237947 -106.5935152393201 -106.1108017566682 -105.4932621297472
 Davidson iteration 0.5 (tol=3.080273116687244e-05): -107.5108722587463 -106.7245930484005 -106.2313070085243 -105.6423237790271
 Davidson iteration 0.6 (tol=2.981303155550973e-06): -107.5108722590151 -106.8008116251279 -106.330411722324 -105.69489630447
 Davidson iteration 0.7 (tol=4.565075388900609e-07): -107.5108722590186 -106.8090569201597 -106.3345525819528 -106.0578325555555
 Davidson iteration 0.8 (tol=4.870689101320012e-08): -107.5108722590186 -106.8098939189322 -106.3973756614146 -106.1617813840531
 Davidson iteration 0.9 (tol=6.224330513865847e-09): -107.5108722590187 -106.8159869231416 -106.4315891597535 -106.2222021964013
 Elapsed time for davidson 0.031276 (sec)
 Elapsed time for diagonalization 0.031286 (sec)
 Elapsed time for mult 0.002511 (sec)
 Energy = -107.5108722590187
 Elapsed time for measurement 4.6e-05 (sec)
 Sample-based diagonalization: Energy = -107.5108722590187
 Sample-based diagonalization: density = [1.999989340109303,1.99999662423008,1.999199121034164,1.981735389885314,1.999825046204946,1.999847567448985,1.995347419300126,0.01136512984757364,0.009844780833945199,0.002849581105572152
 Sample-based diagonalization: carryover bitstrings = [], size = 0

Linia Energy = -107.5108722590187 to oszacowanie stanu podstawowego w spróbkowanej podprzestrzeni.

Czym to się różni od przykładu w Pythonie

Oba przykłady różnią się naraz w kilku aspektach, nie tylko sposobem traktowania ciągów bitów o złej symetrii, więc nie należy oczekiwać identycznych energii:

	Zakładka C++	Zakładka Python
Ciągi bitów o złej symetrii	Odrzucane przez postselect_bitstrings według wagi Hamminga	Naprawiane przez iteracyjne odzyskiwanie konfiguracji
Struktura pętli	Jeden przebieg: próbkowanie → postselekcja → jednokrotna diagonalizacja	Pętla zewnętrzna: diagonalizacja, odczyt obsadzeń orbitali, odzyskanie konfiguracji, ponowna diagonalizacja
Liczba obwodów	100 obwodów (num_circuits = 100)	400 obwodów (200 na czas ewolucji)
Czasy ewolucji	Pojedyncze $\tau$ (time_step = 1)	Dwa czasy, times = [1.0, 10.0]
Wyrazy diagonalne	Zachowane w spróbkowanym operatorze	Usuwane przez filter_diagonal_terms, co zmienia rozkład próbkowania i obwody
Porządek wyrazów	Tylko grupowanie	canonical_order stosowane przed grupowaniem dla powtarzalności
Solver diagonalizacji	SBD (diag, zewnętrzny plik binarny MPI/OpenMP)	solver qiskit-addon-sqd działający w tym samym procesie
Backend	Ustalony po nazwie (backend_name = "ibm_fez"; uruchomienie na tej stronie użyło ibm_pittsburgh)	Operacyjny Backend least_busy z wystarczającą liczbą kubitów
Transpilacja	Każdy obwód osobno, poziom optymalizacji 2	Wszystkie obwody razem, poziom optymalizacji 3
Przesyłanie zadań	100 zadań, po jednym obwodzie	Jedno zadanie Sampler zawierające wszystkie 400 obwodów
Uruchomienie na symulatorze	Brak; tylko sprzęt	Przykład na bezszumowym symulatorze przed uruchomieniem na sprzęcie

Obie wersje używają tego samego hamiltonianu, grupują jego wyrazy według struktury elektronowej, ważą każdą grupę średnim współczynnikiem bezwzględnym jej wyrazów, startują ze stanu Hartree-Focka i próbkują 10 grup oraz 100 strzałów na obwód.

Każda różnica niezależnie wpływa na spróbkowaną podprzestrzeń, a więc i na energię. Większy zespół oraz drugi, dłuższy czas ewolucji poszerzają konfiguracje badane przez przykład w Pythonie; próbkowanie w więcej niż jednym czasie daje SqDRIFT jego podobne do Kryłowa pokrycie przestrzeni niskich energii. Filtrowanie wyrazów diagonalnych zmienia, które operatory może wylosować próbnik, oraz ich wagi, więc oba uruchomienia próbkują różne rozkłady.

Postselekcja jest ścisłym filtrem, a strzał z niewłaściwą liczbą cząstek jest po prostu odrzucany. W powyższym uruchomieniu odrzuciło to większość danych: 1137 z 10000 strzałów przeszło kontrolę wagi Hamminga (7,7) (około 11%), zapadając się do zaledwie 49 różnych determinant $\alpha$. Odzyskiwanie konfiguracji natomiast używa średnich obsadzeń orbitali z wcześniejszej diagonalizacji, aby odwracać bity i naprawiać takie strzały, odzyskując to, co postselekcja by odrzuciła. To sprzężenie zwrotne czyni je z natury iteracyjnym: recover_configurations w qiskit-addon-sqd-hpc wymaga argumentu avg_occupancies, który istnieje dopiero po przeprowadzeniu diagonalizacji.

To uruchomienie C++ daje podprzestrzeń 49 determinant $\alpha$ i energię -107.5109 Ha. Pełna przestrzeń determinant $\alpha$ dla siedmiu elektronów na 10 orbitalach to tylko $\binom{10}{7} = 120$ ciągów, więc podanie wszystkich 120 temu samemu plikowi binarnemu diag daje dokładną energię FCI -107.6482 Ha. Uruchomienie leży więc około 0.1373 Ha powyżej dokładnej wartości, obejmując 49 z 120 determinant $\alpha$. Ponieważ SQD jest wariacyjne w wybranej podprzestrzeni, mniejsza podprzestrzeń może tylko podnieść energię, więc niepełna podprzestrzeń jest najbardziej bezpośrednim wyjaśnieniem różnicy.

Ta kontrola dokładnym FCI jest możliwa tylko dlatego, że przykład jest mały; nie jest częścią przepływu pracy SqDRIFT, a dla rozmiarów układów, w których SQD jest naprawdę potrzebne, takie odniesienie nie istnieje. Służy tu wyłącznie do oszacowania, jaką część przestrzeni to uruchomienie odzyskało.

Porównanie z przykładem w Pythonie nie jest kontrolowane, więc różnicy nie należy przypisywać samemu brakowi odzyskiwania konfiguracji; zauważ, że liczba obwodów, pojedynczy czas ewolucji, filtrowanie wyrazów diagonalnych i szum sprzętowy różnią się naraz i wszystkie kształtują podprzestrzeń.

Aby zmniejszyć tę różnicę, oto kilka możliwych kolejnych kroków do wypróbowania po ukończeniu tego samouczka: zwiększ shots, num_circuits lub ops_per_circuit; próbkuj w więcej niż jednym time_step; i dodaj pętlę odzyskiwania na wierzchu Qiskit::addon::sqd::recover_configurations.

Kolejne kroki

Zalecenia

• Oparta na próbkach kwantowa diagonalizacja Kryłowa fermionowego modelu sieciowego — powiązany samouczek używający obwodów ewolucji w czasie zamiast ansatzu wariacyjnego.

• Oparta na próbkach kwantowa diagonalizacja hamiltonianu chemicznego — samouczek o tym, jak skonstruować obwód lokalnego unitarnego klastra Jastrowa (LUCJ) do symulacji chemii kwantowej.

• Artykuł SqDRIFT — literatura, na której opiera się ten samouczek. (Zauważ, że niektóre optymalizacje omawiane w tym artykule są obecnie w trakcie prac, a ten samouczek może się zmienić w przyszłości w zależności od rozwoju używanych bibliotek.)