Algoritmus SqDRIFT pro odhad základního stavu

Python

Stáhni si notebook v Pythonu

Stáhni si notebook v Pythonu, abys tento tutoriál mohl spustit sám. Notebook stáhne vstup FCIDump, pokud ještě není k dispozici.

Odhad využití: 180 sekund na procesoru Heron r3 (POZNÁMKA: Toto je pouze odhad. Skutečná doba běhu se může lišit.)

Výukové cíle

• Nauč se, jak vytvořit obvody s menší hloubkou ve srovnání s trotterizací

• Projdi si komplexní postup pro odhad základního stavu pomocí qDRIFT a SQD

• Nauč se, jak používat qiskit-fermions společně s dalšími addony Qiskit k implementaci takového postupu

Předpoklady

• Přečti si přehled Sample-based quantum diagonalization (SQD)

• Přečti si lekci Sample-based Krylov Quantum Diagonalization (SKQD)

Pozadí

SqDRIFT je varianta SKQD, která nahrazuje potřebu volby ansatzu, ze kterého se vzorkují bitové řetězce, souborem obvodů časového vývoje sestavených přímo z cílového hamiltoniánu. Toho je dosaženo subsamplingem menších operátorů časového vývoje z hamiltoniánu na základě jeho koeficientů, což je známo jako trotterizační metoda qDRIFT.

Obě implementace používají Qiskit Fermions k reprezentaci molekulárního hamiltoniánu a k sestavení obvodů pro algoritmus qDRIFT. Implementace v Pythonu používá fermionové průchody transpilátoru, zatímco implementace v C++ explicitně vzorkuje skupiny a mapuje je na qubitové operátory.

Nechť má hamiltonián tvar:

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

kde, bez ztráty obecnosti, požadujeme $c_i > 0$ a aby největší vlastní hodnota $h_i$ byla co do absolutní hodnoty rovna $1$. Jakýkoli znaménkový nebo komplexní prefaktor je absorbován do $h_i$, takže koeficienty $c_i$ jsou striktně kladné váhy, zatímco $h_i$ nesou směr každého členu. Zde $N$ je počet členů (nebo po seskupení počet skupin) v hamiltoniánu; je to vlastnost hamiltoniánu a liší se od počtu operátorů vzorkovaných do jednoho obvodu, dále značeného $n$.

Algoritmus qDRIFT pak pro cílový čas $t$ realizuje nějaký operátor $V_k$, kde $k$ probíhá od $1 \cdots K$ a označuje $k$-tý obvod SqDRIFT, definovaný jako:

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

Zde $n$ je počet vzorkovaných operátorů na obvod a $K$ je počet obvodů v souboru. Součin probíhá přes $n$ tahů, nikoli přes všechny $N$ členy hamiltoniánu, a protože jsou členy vybírány s vracením, stejný $h_i$ se může v jednom $V_k$ objevit vícekrát.

Veličina:

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

je $L_1$ norma koeficientů, takže každý z $n$ kroků probíhá po stejnou dobu $\lambda t / n$ bez ohledu na to, který člen byl vybrán. Jednotnost úhlu kroku je charakteristickým rysem qDRIFT: koeficient ovlivňuje výsledek prostřednictvím toho, jak často je jeho člen vybrán, nikoli tím, o kolik je tento člen otočen. Indexy jsou vzorkovány z rozdělení:

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

takže řada $(k_1, \ldots, k_n)$ je náhodná posloupnost indexů členů vybraných z tohoto rozdělení. Protože $c_i$ jsou kladné a sčítají se na $\lambda$, jde o normalizované pravděpodobnostní rozdělení a očekávaná hodnota výsledného kanálu přes náhodné výběry aproximuje vývoj podle $H$ s chybou, která klesá s rostoucím $n$. Všimni si, že chyba aproximace závisí na $\lambda$, nikoli na počtu členů $N$.

(Článek o SqDRIFT zapisuje počet členů jako $\mathcal{N}$ a délku posloupnosti jako $N$; zde používáme $N$ a $n$, abychom ty dvě veličiny jasně odlišili.)

Tento tutoriál ukazuje, jak vygenerovat soubor takových randomizovaných obvodů. Poté, co jsme tyto obvody vytvořili, podobně jako když vytváříme Krylovův podprostor pro různé operátory, vzorkujeme bitové řetězce z více takových operátorů s různými časovými parametry. To zajišťuje vyšší překryv mezi vektory základního stavu a vzorkovanými bitovými řetězci.

Požadavky

Před zahájením tohoto tutoriálu se ujisti, že máš nainstalováno

• Virtuální prostředí Python (>=3.10)
• pip>=25.1
• qiskit ~= 2.5
• qiskit-fermions==0.1.0 (všimni si, že název je v množném čísle)
• numpy
• pyscf
• qiskit-aer
• qiskit-ibm-runtime
• qiskit-addon-sqd

Všechny požadované balíčky můžeš nainstalovat pomocí:

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

Nastavení

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

Implementace

Krok 1: Namapuj klasické vstupy na kvantový problém

Níže uvedené kroky v Pythonu používají simulátor Aer. Hardwarový příklad následuje poté.

Čtení a příprava FCIDump

Pro tento tutoriál načteme hamiltonián elektronové struktury pro dusík (N2). Existují i jiné způsoby, jak vytvořit fermionické operátory. Podívej se na dokumentaci qiskit_fermions.operators.library.

O tomto FCIDump. Soubor N2_sto_3g popisuje molekulu dusíku ($N_2$) v minimální bázi STO-3G při mezijaderné vzdálenosti 1,09 $\AA$, což je experimentální rovnovážná délka vazby. Jeho hlavička deklaruje NORB=10, NELEC=14 a MS2=0: 10 prostorových orbitalů (tedy 20 spinových orbitalů a 20 qubitů podle Jordan-Wignerovy transformace), 14 elektronů ve spinovém singletu, tedy sedm $\alpha$ a sedm $\beta$ elektronů. Všechny orbitaly mají přiřazen štítek symetrie 1, tedy nevyužívá se žádná symetrie bodové grupy. Jelikož jde o STO-3G výpis celého prostoru, žádné orbitaly nejsou zmrazené a korelační prostor je dostatečně malý na to, aby bylo možné klasicky spočítat přesnou referenční energii FCI pro porovnání, jak je ukázáno v následujícím bloku kódu.

Ekvivalentní soubor lze znovu vygenerovat 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")

Protože integrály závisí na konvergovaných SCF orbitalech, znovu vygenerovaný soubor se může od dodaného lišit ve fázi nebo pořadí orbitalů; celkové energie tím nejsou ovlivněny.

Získání souboru. FCIDump najdeš v tomto GitHub repozitáři. Spuštěním níže uvedeného kódu jej můžeš stáhnout do umístění, které zbytek tutoriálu očekává.

Nejprve použijeme cisolver poskytovaný knihovnou pyscf k získání referenční energie. Jde o skutečnou energii základního stavu molekuly, se kterou pracujeme. Za tímto účelem nejprve deklarujeme norb a nelec, což je počet orbitalů a počet elektronů. Poté deklarujeme h1e a h2e, což jsou jedno- a dvouelektronové integrály. Všechny tyto hodnoty budou později použity i pro 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

Načtení hamiltoniánu

S připravenými potřebnými daty načteme hamiltonián ze souboru FCI ve formátu kompatibilním s qiskit-fermions

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

Fermionové pracovní postupy s qiskit-fermions

Nejprve namapujeme hamiltonián do modelu fermionového obvodu pomocí qiskit-fermions, který poskytuje transpilátorové průchody a hradla specifické pro fermionové obvody. Ty budou později použity před tradičními transpilátorovými průchody Qiskitu v tomto pracovním postupu.

Seskupování členů

Abychom zajistili reprodukovatelnost výsledků, nejprve použijeme canonical_order k seřazení členů pouze na základě jejich struktury. Pořadí operátorů v seznamu canon je tak pevně dané. To zajišťuje reprodukovatelnost vytvořených operátorů, protože průchod QDriftTrotterization, který později použijeme, náhodně vzorkuje indexy pro vytvoření operátorů qDRIFT.

V tomto kroku využíváme mnoha symetrií přítomných v hamiltoniánu elektronové struktury tím, že seskupíme související členy se stejnými koeficienty. I když tím dojde ke změně rozdělení koeficientů operátoru, ze kterého protokol qDRIFT vzorkuje, nemá to vliv na jeho záruky konvergence. Zásadní je, že seskupení členů souvisejících symetrií vede k příznivému rušení Pauliho členů a k celkově kratší hloubce obvodu při časovém vývoji stavu pod jejich působením.

qiskit-fermions poskytuje funkci group_terms_by_electronic_structure, která za nás toto seskupení provede.

Všimni si, že group_terms_by_electronic_structure předpokládá členy v normálním uspořádání.

Filtrování diagonálních členů

Odstraníme diagonální členy z hamiltoniánu použitého k vygenerování obvodů, aby bylo $n$ vzorkovacích slotů qDRIFT využito na členy, které přesouvají populaci mezi konfiguracemi. Takové členy je nejlepší z hamiltoniánu odfiltrovat právě v tomto bodě, ještě před sestavením hradla Evolution v dalším kroku.

Jde o členy, které jsou diagonální v bázi obsazovacích čísel, tedy o součiny číselných operátorů $a^\dagger_i a_i$. Pod tento popis spadají tři druhy členů:

• konstantní posun energie, součin nulového počtu číselných operátorů, jehož časový vývoj přispívá pouze globální fází;

• jednotlivé číselné operátory $n_i$, jejichž časový vývoj se redukuje na jednoqubitové rotace $Z$;

• součiny vyššího řádu, jako je $n_i n_j$.

Samy o sobě žádné z nich nepřesouvají populaci mezi konfiguracemi obsazovacích čísel; působí pouze na fáze již přítomných konfigurací. Nejsou však neutrální: tyto relativní fáze vstupují do interference generované excitačními členy později v obvodu, takže jejich odfiltrování mění skutečně generovaný vývoj a může změnit rozdělení vzorkování. Jde o záměrnou aproximaci v kroku generování obvodu, provedenou proto, aby se vzorkování soustředilo na excitační členy, nikoli o krok, který ponechává vzorkované rozdělení nedotčené. Na rozdíl od výše uvedeného seskupování podle symetrie, které ponechává záruky konvergence qDRIFT nedotčené, tento filtr mění operátor, jehož vývoj se počítá. Obvody proto už neaproximují vývoj pod plným hamiltoniánem a chybové meze qDRIFT platí pro filtrovaný operátor, nikoli pro původní. To je zde přijatelné, protože obvody jsou jen vzorkovací heuristikou používanou k navrhování konfigurací: ze samotného odhadu energie se žádný člen neztratí, protože filtr se vztahuje pouze na hamiltonián použitý k sestavení obvodů, zatímco pozdější klasická diagonalizace používá plný hamiltonián včetně diagonálních členů. Přesnost SQD závisí na tomto klasickém kroku, který zůstává variační ve vzorkovaném podprostoru bez ohledu na to, jak byly konfigurace navrženy.

Funkce filter_diagonal_terms() odstraní takové členy z operátoru na místě (in place). Identifikuje je podle jejich normálně uspořádané struktury — multimnožina kreačních módů odpovídá multimnožině anihilačních módů — takže je platná pouze pro operátor, který je již normálně uspořádaný. Tento předpoklad se za běhu nekontroluje.

# 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

Nyní, když jsme seskupili členy v hamiltoniánu, se rozhodneme pro následující parametry pro vygenerování souboru obvodů:

• Počet obvodů, které se mají vygenerovat: num_circuits
• Délka každého obvodu z hlediska excitovaných skupin: num_exc
• Faktor pro různé časy vývoje: times

Vytváření fermionových obvodů

Nyní vytvoříme fermionové obvody pro každý z časových kroků. Každý obvod bude sestávat z jediného evolučního hradla s časem vývoje, který jsme deklarovali dříve. Operátorem vývoje je hamiltonián. Později na tyto obvody spustíme transpilátorové průchody a vytvoříme obvody qDRIFT.

Příprava ansatzu

Stav Hartree-Fock připravíme pomocí třídy InitializeModes. Pro dusík tento proces spočívá jednoduše v aplikaci hradel X na prvních num_elec_a qubitů a poté na num_elec_b qubitů, přičemž obě hodnoty jsou pro dusík rovny sedmi. Tento stav představuje sedm $\alpha$ a sedm $\beta$ elektronů dusíku.

# 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: Optimalizace úlohy pro spuštění na kvantovém hardwaru

Nyní, když máme naše obvody, nejprve použijeme průchody dostupné v qiskit-fermions k provedení optimalizací na fermionové úrovni a poté náš obvod transpilujeme pro zvolený backend. Jelikož jde o experiment na simulátoru, provedeme to nejprve pro AerSimulator.

Výpočet váhy pro každou skupinu

V tomto kroku provedeme stochastické vzorkování členů metodou qDRIFT s pravděpodobnostmi úměrnými jejich koeficientům v hamiltoniánu. Provede to za nás transpilátorový průchod qDRIFT. Nyní můžeme vytvořit mělčí obvody, které lze efektivněji spouštět na hardwaru i přes omezenou konektivitu qubitů, a to i když hamiltonián obsahuje vazby na dlouhé vzdálenosti a členy vyššího než kvadratického řádu. Po seskupení členů vzorkuje operátory na základě jejich vah. Pro každý operátor $h_i$ je váha $W_{h_i}$ definována následovně:

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

Protože byly členy seskupeny v kroku 1, je každé $h_i$ zde celou skupinou: $c_i$ je střední absolutní koeficient členů ve skupině $i$ a každý člen ve skupině se vyvíjí s koeficientem zredukovaným na jeho znaménko.

Fermionové optimalizace a optimalizace nativní pro hardware

Funkce generate_preset_jw_pass_manager() vrací MultiStagePassManager, který přijímá FermionicCircuit a produkuje optimalizovaný výsledný obvod, jenž můžeme transpilovat pro spuštění na našem hardwaru. Jeho výchozí optimalizační fázi nahradíme za FermionicPassManager obsahující náš průchod QDriftTrotterization:

• Průchod QDriftTrotterization interně používá výpočet vah a vzorkování k vygenerování obvodů, které použijeme pro vzorkování

• Průchod RelabelModes je další optimalizační průchod, který lze použít k permutaci fermionových módů za účelem optimalizace konektivity mezi qubity a snížení hloubky hradel; více se dozvíš v referenci API

Zbývající fáze MultiStagePassManager běží automaticky a zajišťují úplné mapování fermionů na qubity:

• F2QLayout: Přednastavený pass manager aplikuje průchod TrivialF2QLayout, který triviálně mapuje $n$ fermionových bitů na $n$ qubitů.

• F2QSynth: Transpilační průchod, který mapuje instrukce fermionového obvodu na instrukce založené na qubitech.

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

Nyní, když jsme dokončili optimalizace na fermionové úrovni, můžeme obvody transpilovat pro spuštění na simulátoru.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Krok 3: Spuštění pomocí primitiv Qiskitu

Nyní, když máme naše obvody, je můžeme spustit pomocí primitiv Qiskitu na AerSimulator. Spojíme všechny počty (counts) z různých obvodů. Před závěrečným postprocessingem pomocí SQD je převedeme na booleovské vektory.

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: Postprocessing a vrácení výsledku v požadovaném klasickém formátu

Použití bitových řetězců pro SQD

Nyní můžeme spustit diagonalizační schéma na vybraných bitových řetězcích a najít nejnižší vlastní hodnotu, která bude odpovídat energii základního stavu molekuly. Vytvoříme callback funkci, deklarujeme počáteční obsazenosti a nastavíme parametry, než nakonec spustíme diagonalizační schéma. Callback funkce se používá k výpisu aktuální iterace a aktuálního odhadu vlastní hodnoty v každé iteraci.

Nakonec, abychom získali odhad energie základního stavu, přičteme nuclear_repulsion_energy k výsledné energii.

Poznámka: Dimenze podprostoru není napříč iteracemi pevná, a to ani na bezšumovém simulátoru — každý podvzorek vybere jinou sadu konfigurací a krok obnovy (recovery) mezi iteracemi přetváří fond, takže se uváděná dimenze mezi jednotlivými podvzorky liší. Bezšumové vzorkování samo o sobě dimenzi vybraného podprostoru nefixuje. Spuštění na hardwaru však má tendenci dávat systematicky větší podprostory, protože šumové výstřely (shots) narušují symetrii počtu částic a obnova konfigurací je mění na další bázové vektory. Z tohoto důvodu v sekci o hardwaru zavedeme ještě další krok pro prořezávání (pruning) bitových řetězců.

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

Příklad na hardwaru

Tento příklad používá 20 qubitů (10 prostorových orbitalů). Tato volba je dána tím, že tutoriál má proběhnout rychle, nikoli pevným stropem metody.

Náklady na klasický krok nejsou určeny přímo počtem qubitů. SQD diagonalizuje hamiltonián projektovaný do podprostoru rozepnutého vzorkovanými konfiguracemi, takže klasické náklady určuje dimenze tohoto vybraného podprostoru — ta je zde řízena parametry samples_per_batch, num_batches a tím, kolik různých konfigurací obvody skutečně vyprodukují — společně s řídkou lineární algebrou potřebnou k aplikaci projektovaného hamiltoniánu. Plný prostor CI roste kombinatoricky s počtem orbitalů a elektronů, ale vybraný podprostor je jeho malým, laditelným výsekem a jeho velikost řídíme přímo. Počet qubitů a klasickou náročnost je tedy možné do jisté míry měnit nezávisle na sobě: širší prostor orbitalů vzorkovaný do skromného podprostoru může být levnější než menší systém diagonalizovaný přes velmi velký podprostor.

V praxi tedy proveditelná velikost systému závisí na dimenzi podprostoru, kterou potřebuješ pro požadovanou přesnost, a na paměti a jádrech dostupných pro eigensolver. Větší prostory orbitalů typicky vyžadují větší podprostor pro dosažení chemické přesnosti, a to je nakonec to, co motivuje distribuované prostředky — informace o škálování tohoto kroku najdeš v qiskit-addon-sqd-hpc. Namísto předpokladu pevného prahu je praktický přístup sledovat uváděnou dimenzi podprostoru a konvergenci energie napříč iteracemi a zvětšovat velikost podprostoru, dokud se energie přestane zlepšovat nebo dokud nevyčerpáš dostupnou paměť.

Poznámka: Kvůli chybě vzorkování způsobené šumem na hardwaru bude podprostor vytvořený pro diagonalizaci při spuštění na hardwaru větší než ten, který získáme při použití simulátoru. I když se tím zvýší dimenze podprostoru, který chceme diagonalizovat, pracovní postup nám díky odolnosti SQD vůči šumu stále poskytuje přesnou odpověď.

Prořezávání falešných řetězců

Zde se můžeme rozhodnout provést další krok. Jakmile máme všechny bitové řetězce ze spuštění obvodů, můžeme buď před spuštěním SQD odfiltrovat neplatné bitové řetězce, nebo pokračovat bez prořezávání. Vynechání prořezávání je pro spuštění na hardwaru obecně výhodnější, protože ponechává výstřely s porušenou symetrií k dispozici pro obnovu konfigurací, která je může opravit na platné konfigurace a tím rozšířit podprostor, místo aby tyto výstřely rovnou zahodila.

Jelikož dusík může mít pouze sedm $\alpha$ a sedm $\beta$ elektronů, lze zahodit jakékoli bitové řetězce, které mají v první nebo druhé polovině výstupu více nebo méně než sedm jedniček. Definujeme funkci, která kontroluje, zda jsou bitové řetězce platné, a pokud ne, zahodí je. Jakmile odfiltrujeme falešné bitové řetězce, zbytek se pošle do diagonalizačního schématu. Pomocí příznaku PRUNE níže přepínej mezi oběma způsoby chování.

Měj na paměti, že prořezávání je jen jednou z několika voleb, které formují výsledný podprostor, vedle počtu obvodů, sady časů vývoje a filtrování diagonálních členů. Porovnání spuštění s prořezáváním a bez něj je vypovídající pouze tehdy, pokud je vše ostatní ponecháno beze změny; podrobněji to rozebírá karta C++, jelikož ta místo obnovy provádí postselekci a liší se i v dalších parametrech.

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

Získej kód v C++

Zdrojový kód v C++, projekt CMake a pokyny k sestavení jsou v adresáři docs/tutorials/assets/sqdrift/cpp repozitáře dokumentace Qiskit. Vstup FCIDump je v sousedním adresáři fcidump_files. Naklonuj repozitář, nebo si stáhni oba adresáře, abys mohl tento tutoriál sestavit a spustit sám.

Odhad využití: 180 sekund na procesoru Heron r3 (POZNÁMKA: Toto je pouze odhad. Skutečná doba běhu se může lišit.)

Výukové cíle (C++)

• Nauč se, jak vytvořit obvody s menší hloubkou ve srovnání s trotterizací

• Projdi si komplexní postup pro odhad základního stavu pomocí qDRIFT a SQD

• Nauč se, jak používat qiskit-fermions společně s dalšími addony Qiskit k implementaci takového postupu

Předpoklady (C++)

• Přečti si přehled Sample-based quantum diagonalization (SQD)

• Přečti si lekci Sample-based Krylov Quantum Diagonalization (SKQD)

Pozadí (C++)

SqDRIFT je varianta SKQD, která nahrazuje potřebu volby ansatzu, ze kterého se vzorkují bitové řetězce, souborem obvodů časového vývoje sestavených přímo z cílového hamiltoniánu. Toho je dosaženo subsamplingem menších operátorů časového vývoje z hamiltoniánu na základě jeho koeficientů, což je známo jako trotterizační metoda qDRIFT.

Obě implementace používají Qiskit Fermions k reprezentaci molekulárního hamiltoniánu a k sestavení obvodů pro algoritmus qDRIFT. Implementace v Pythonu používá fermionové průchody transpilátoru, zatímco implementace v C++ explicitně vzorkuje skupiny a mapuje je na qubitové operátory.

Nechť má hamiltonián tvar:

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

kde, bez ztráty obecnosti, požadujeme $c_i > 0$ a aby největší vlastní hodnota $h_i$ byla co do absolutní hodnoty rovna $1$. Jakýkoli znaménkový nebo komplexní prefaktor je absorbován do $h_i$, takže koeficienty $c_i$ jsou striktně kladné váhy, zatímco $h_i$ nesou směr každého členu. Zde $N$ je počet členů (nebo po seskupení počet skupin) v hamiltoniánu; je to vlastnost hamiltoniánu a liší se od počtu operátorů vzorkovaných do jednoho obvodu, dále značeného $n$.

Algoritmus qDRIFT pak pro cílový čas $t$ realizuje nějaký operátor $V_k$, kde $k$ probíhá od $1 \cdots K$ a označuje $k$-tý obvod SqDRIFT, definovaný jako:

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

Zde $n$ je počet vzorkovaných operátorů na obvod a $K$ je počet obvodů v souboru. Součin probíhá přes $n$ tahů, nikoli přes všechny $N$ členy hamiltoniánu, a protože jsou členy vybírány s vracením, stejný $h_i$ se může v jednom $V_k$ objevit vícekrát.

Veličina:

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

je $L_1$ norma koeficientů, takže každý z $n$ kroků probíhá po stejnou dobu $\lambda t / n$ bez ohledu na to, který člen byl vybrán. Jednotnost úhlu kroku je charakteristickým rysem qDRIFT: koeficient ovlivňuje výsledek prostřednictvím toho, jak často je jeho člen vybrán, nikoli tím, o kolik je tento člen otočen. Indexy jsou vzorkovány z rozdělení:

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

takže řada $(k_1, \ldots, k_n)$ je náhodná posloupnost indexů členů vybraných z tohoto rozdělení. Protože $c_i$ jsou kladné a sčítají se na $\lambda$, jde o normalizované pravděpodobnostní rozdělení a očekávaná hodnota výsledného kanálu přes náhodné výběry aproximuje vývoj podle $H$ s chybou, která klesá s rostoucím $n$. Všimni si, že chyba aproximace závisí na $\lambda$, nikoli na počtu členů $N$.

(Článek o SqDRIFT zapisuje počet členů jako $\mathcal{N}$ a délku posloupnosti jako $N$; zde používáme $N$ a $n$, abychom ty dvě veličiny jasně odlišili.)

Tento tutoriál ukazuje, jak vygenerovat soubor takových randomizovaných obvodů. Poté, co jsme tyto obvody vytvořili, podobně jako když vytváříme Krylovův podprostor pro různé operátory, vzorkujeme bitové řetězce z více takových operátorů s různými časovými parametry. To zajišťuje vyšší překryv mezi vektory základního stavu a vzorkovanými bitovými řetězci.

Požadavky (C++)

Tato implementace běží na hardwaru IBM Quantum® a pro klasickou diagonalizaci používá externí spustitelný soubor SBD. Neobsahuje příklad se simulátorem v Pythonu.

Podle pokynů k sestavení v C++ nainstaluj překladač, CMake, Rust, Boost a závislosti Qiskit. Rust nainstaluj pomocí rustup: závislosti Qiskit mají připnutou vlastní verzi Rustu, kterou rustup stáhne automaticky, ale Rust nainstalovaný z Homebrew nebo správce balíčků Linuxu může být příliš starý. Závěrečný krok SBD navíc vyžaduje MPI, OpenMP a BLAS/LAPACK. Windows tato implementace aktuálně nepodporuje.

Použij úplný zdrojový kód v C++ a projekt CMake z repozitáře. Při stahování příkladu ponech sousední adresář fcidump_files. Z kopie repozitáře dokumentace sestav projekt následujícími příkazy. První sestavení stáhne a přeloží všechny závislosti a trvá několik minut.

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

Sestavení zapíše spustitelný soubor SqDRIFT do docs/tutorials/assets/sqdrift/cpp (kořen projektu), nikoli do build/. Spusť ./SqDRIFT z kořene projektu, aby našel ../fcidump_files/N2_sto_3g. Všechny ostatní příkazy na této kartě se rovněž spouštějí z kořene projektu, pokud není uvedeno jinak. Pokud ./SqDRIFT selže s chybou chybějící sdílené knihovny, postupuj podle pokynů k cestě ke knihovnám pro tvou platformu v průvodci sestavením; na macOS obvykle spustitelný soubor běží i bez nich.

Nastavení (C++)

Níže uvedené úryvky jsou výňatky z úplného programu SqDRIFT.cpp. Tento soubor sestav a spusť jako jeden program.

Před spuštěním SqDRIFT se ujisti, že je tvůj API token IBM Quantum nakonfigurován pro IBM Cloud Platform.

Pokud máš Python s nainstalovaným qiskit-ibm-runtime, nejsnazší způsob, jak vytvořit správně formátovaný konfigurační soubor, je:

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

Tím se zapíše ~/.qiskit/qiskit-ibm.json (Linux/macOS) nebo %USERPROFILE%\.qiskit\qiskit-ibm.json (Windows) se správnou strukturou. Pokud soubor raději vytvoříš ručně:

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

API klíč získáš na IBM Quantum Platform.

Poznámka: Chybějící nebo chybně formátovaný konfigurační soubor způsobí při spuštění nepříliš nápomocnou paniku Rustu unwrap() místo čisté chybové zprávy. Pokud narazíš na problém, zkontroluj název klíče (default-ibm-quantum-platform) a hodnotu kanálu (ibm_quantum_platform).

Implementace (C++)

Krok 1: Namapuj klasické vstupy na kvantový problém (C++)

Načtení molekulárního hamiltoniánu

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)

Načte hamiltonián molekuly $N_2$ ze souboru FCIDump. Počet qubitů se rovná dvojnásobku počtu orbitalů (spin nahoru a spin dolů).

Vytvoření hamiltoniánu a jeho normální uspořádání

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

Převede data FCIDump na fermionový operátor a použije normální uspořádání, aby zjednodušil strukturu operátoru.

Seskupení členů podle elektronové struktury

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

Seskupí související členy hamiltoniánu podle elektronové struktury před sestavením jejich obvodů vývoje.

Výpočet vzorkovacích vah, normalizace a mapování na qubitové operátory

Na pořadí zde záleží. Normalizace přeškáluje každý koeficient na jednotkovou velikost na místě, takže vzorkovací váhy je nutné vypočítat jako první. Pokud normalizace proběhne první, každé $|c_j|$ se načte zpět jako přesně 1.0 a váha každé skupiny se zredukuje na pouhý počet jejích členů — vzorkovací rozdělení by pak bylo řízeno velikostí skupiny, nikoli fyzikálními koeficienty.

Váha každé skupiny je střední absolutní koeficient jejích členů. Průchod QDriftTrotterization na kartě Python používá stejnou váhu pro seskupený hamiltonián. Použití součtu by způsobilo, že by $\lambda$ rostlo s velikostí skupiny, což by zase zkreslilo čas vývoje $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

Váha $\lambda_i = \frac{1}{|G_i|}\sum_{j \in G_i} |c_{ij}|$ je střední absolutní koeficient přes členy skupiny $G_i$, čtený z původních koeficientů. Normalizace pak před Jordanovou-Wignerovou transformací přeškáluje každý fermionový člen na jednotkovou velikost a ve vyvíjeném operátoru zachová pouze jeho fázi/znaménko. To odpovídá vývoji normalizovaných členů v qDRIFT: velikosti koeficientů určují pravděpodobnosti vzorkování prostřednictvím $\lambda$, zatímco vývoj obvodu používá normalizovaný člen, takže velké koeficienty nejsou započítány dvakrát.

Vzorkování sad operátorů

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

Vylosuje sady operátorů pro 100 náhodných obvodů, z nichž každý obsahuje 10 skupin vzorkovaných z váženého rozdělení. Další krok z těchto sad vytvoří soubor obvodů qDRIFT.

Krok 2: Sestavení obvodů qDRIFT (C++)

Ve verzi v C++ tento krok sestaví obvody ze vzorkovaných sad operátorů. Transpilace pro cílový backend probíhá v kroku 3, bezprostředně před odesláním každého obvodu.

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

Úplný program poté vypíše statistiky obvodů a indexy operátorů obvodu 0, které jsou zde vynechány.

Pro každý vzorkovaný operátor vytvoří obvod Suzukiho-Trotterova vývoje s časem $t_k = \frac{\lambda \tau}{n}$, kde $n$ je počet operátorů vzorkovaných na obvod a $\tau$ je časový krok. Protože každá vzorkovaná skupina byla před mapováním normalizována člen po členu, tento vývoj aplikuje během rotace pouze fázi/znaménko každého členu; velikosti přispívají jen prostřednictvím vzorkovacího rozdělení. Počáteční Hartreeův-Fockův stav se připraví aplikací hradel X na qubity 0–6 a 10–16 (sedm elektronů v každém spinovém sektoru). Instrukce každého obvodu vývoje se ručně připojí k hlavnímu obvodu pomocí API C++.

Krok 3: Spuštění pomocí primitiv Qiskitu (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;

Připojí se k IBM Quantum, vybere Backend podle názvu, transpiluje pro něj každý obvod a odešle jednu úlohu na obvod: 100 samostatných úloh Sampler po 100 měřeních na 20qubitovém registru. Poté každou úlohu dotazuje, dokud se nedokončí, a z result[0] shromáždí bitové řetězce. Každá úloha spotřebuje čas na QPU a ve frontě čeká samostatně, takže počítej s tím, že tento krok potrvá mnohem déle než jedna dávková úloha.

Hardwarové výsledky uvedené na této stránce pocházejí ze spuštění na ibm_pittsburgh (100 obvodů × 100 měření), přičemž backend_name byl změněn z výchozího ibm_fez. Uložené vzorky byly zpracovány s opraveným pořadím bitů determinantů uvedeným v kroku 4. Hardwarový šum se mezi spuštěními liší, takže tvé počty a konečná energie se budou lišit.

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

Změna Backendu: Změň backend_name v SqDRIFT.cpp, abys cílil na jiné zařízení; service.backends() vypíše, k čemu má tvůj účet přístup.

Krok 4: Postprocessing a vrácení výsledku v požadovaném klasickém formátu (C++)

Postselekce bitových řetězců a zápis řetězců 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!

Provede postselekci bitových řetězců se správnou Hammingovou vahou (n_alpha elektronů se spinem nahoru, n_beta se spinem dolů, odvozeno z pole NELEC ve FCIDump), převede je na řetězce Configuration Interaction (CI) pomocí spinové symetrizace a zapíše je do alphadets_from_sqd.txt pro následnou diagonalizaci ve vybrané bázi.

Uvolnění prostředků

// 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;

Řádně uvolní veškerou alokovanou paměť a prostředky. Vektor std::vector<QuantumCircuit> v C++ se uklidí automaticky díky RAII. Ručně uvolní qubitové operátory, fermionové operátory a data FCIDump.

Když jsme nyní vytvořili bázi, do níž budeme promítat náš hamiltonián, můžeme přistoupit k diagonalizaci, abychom získali odhad základního stavu. Následující příkazy spusť z kořene projektu, docs/tutorials/assets/sqdrift/cpp.

Sestavení binárního souboru SBD diag

CMake stáhne SBD pouze jako zdrojový kód do deps/sbd/ (hlavní sestavení ho záměrně nepřekládá). Aplikace SBD používá obyčejný Makefile a samotné SBD vyžaduje MPI, OpenMP a BLAS/LAPACK — takže obyčejné make nestačí, dokud nejsou nainstalovány a dokud soubor Configuration aplikace neodpovídá tvé sadě nástrojů.

Makefile čte příznaky překladače a linkeru ze souboru Configuration ve stejném adresáři. Verze dodávaná upstreamem je vyladěna pro jedno konkrétní uspořádání macOS/Homebrew a na běžném počítači se tak, jak je, nesestaví — zejména -I/opt/homebrew/opt/llvm/include neobsahuje omp.h (ten leží v adresáři zdrojů clangu specifickém pro verzi) a clang++ od Applu příznak -fopenmp rovnou odmítne:

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

Následující postup je ověřen jako funkční na macOS (Apple Silicon, Darwin 25.6, Open MPI 5.x + Homebrew LLVM).

macOS (Apple Silicon) — otestováno

# 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

Poté nahraď deps/sbd/apps/chemistry_tpb_selected_basis_diagonalization/Configuration tímto:

# 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 ../../../..

Úspěšné sestavení vypíše dva řádky překladu/linkování a vedle Makefile zanechá binární soubor diag. Jedno varování ld: warning: ignoring duplicate libraries: '-lomp' je neškodné.

Intel macOS: Prefix Homebrew je /usr/local místo /opt/homebrew, takže ho nahraď v obou předchozích cestách.

Linux (Ubuntu/Debian) — neotestováno

Následující předpoklady pokrývají MPI, OpenMP a BLAS/LAPACK. GCC podporuje -fopenmp nativně, takže není potřeba žádné přepsání překladače. Na Linuxu to nebylo ověřeno:

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

Spuštění SqDRIFT pro vygenerování řetězců CI

Pokud jsi už spustil ./SqDRIFT v krocích 1–4, soubor alphadets_from_sqd.txt už v kořeni projektu existuje a tento krok můžeš přeskočit. Jinak ho spusť z kořene projektu:

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

Spuštění diagonalizace

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: Cesta k souboru FCIDump obsahujícímu hamiltonián
• --adetfile: Cesta k souboru alfa determinantů (řetězce CI ze SqDRIFT)
• --method 0: Davidsonova diagonalizační metoda
• --iteration 100: Maximální počet Davidsonových iterací
• --block 10: Velikost bloku pro Davidsonův algoritmus
• --tolerance 1e-8: Tolerance konvergence energie

 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

Řádek Energy = -107.5108722590187 je odhad základního stavu ve vzorkovaném podprostoru.

Jak se to liší od příkladu v Pythonu

Oba příklady se liší hned v několika ohledech, nejen v tom, jak zacházejí s bitovými řetězci se špatnou symetrií, takže se nepředpokládá, že dají stejné energie:

	Karta C++	Karta Python
Bitové řetězce se špatnou symetrií	Zahozeny funkcí postselect_bitstrings podle Hammingovy váhy	Opraveny iterativní obnovou konfigurací
Struktura smyčky	Jeden průchod: vzorkování → postselekce → jedna diagonalizace	Vnější smyčka: diagonalizace, načtení orbitalových obsazení, obnova konfigurací, opětovná diagonalizace
Počet obvodů	100 obvodů (num_circuits = 100)	400 obvodů (200 na čas vývoje)
Časy vývoje	Jediné $\tau$ (time_step = 1)	Dva časy, times = [1.0, 10.0]
Diagonální členy	Ponechány ve vzorkovaném operátoru	Odstraněny funkcí filter_diagonal_terms, což mění vzorkovací rozdělení i obvody
Řazení členů	Pouze seskupeno	Před seskupením se pro reprodukovatelnost použije canonical_order
Diagonalizační řešič	SBD (diag, externí binární soubor MPI/OpenMP)	Řešič qiskit-addon-sqd běžící v procesu
Backend	Pevně daný názvem (backend_name = "ibm_fez"; běh na této stránce použil ibm_pittsburgh)	Provozní Backend least_busy s dostatkem qubitů
Transpilace	Každý obvod zvlášť, úroveň optimalizace 2	Všechny obvody společně, úroveň optimalizace 3
Odesílání úloh	100 úloh, každá s jedním obvodem	Jedna úloha Sampler obsahující všech 400 obvodů
Běh na simulátoru	Žádný; pouze hardware	Příklad s bezšumovým simulátorem před během na hardwaru

Obě verze používají stejný hamiltonián, seskupují jeho členy podle elektronové struktury, váží každou skupinu středním absolutním koeficientem jejích členů, začínají z Hartreeova-Fockova stavu a vzorkují 10 skupin a 100 měření na obvod.

Každý rozdíl ovlivňuje vzorkovaný podprostor, a tím i energii, nezávisle. Větší soubor obvodů i druhý, delší čas vývoje rozšiřují konfigurace, které příklad v Pythonu prozkoumává; vzorkování ve více než jednom čase dává SqDRIFT jeho krylovovské pokrytí nízkoenergetického prostoru. Filtrování diagonálních členů mění, které operátory může vzorkovač losovat, i jejich váhy, takže obě spuštění vzorkují různá rozdělení.

Postselekce je přísný filtr a měření se špatným počtem částic se rovnou zahodí. V uvedeném běhu to zahodilo většinu dat: kontrolou Hammingovy váhy (7,7) prošlo 1137 z 10000 měření (asi 11 %), která se sloučila do pouhých 49 různých $\alpha$-determinantů. Obnova konfigurací naopak používá průměrná orbitalová obsazení z předchozí diagonalizace k překlápění bitů a takové vzorky opravuje, takže znovu využije to, co by postselekce zahodila. Tato zpětná vazba z ní činí z podstaty iterativní postup: recover_configurations v qiskit-addon-sqd-hpc vyžaduje argument avg_occupancies, který existuje až po proběhnuté diagonalizaci.

Tento běh v C++ dává podprostor 49 $\alpha$-determinantů a energii -107.5109 Ha. Úplný prostor $\alpha$-determinantů pro sedm elektronů v 10 orbitalech má pouze $\binom{10}{7} = 120$ řetězců, takže předložením všech 120 stejnému binárnímu souboru diag získáme přesnou energii FCI -107.6482 Ha. Běh se tedy nachází asi 0.1373 Ha nad přesnou hodnotou, přičemž pokrývá 49 ze 120 $\alpha$-determinantů. Protože je SQD ve vybraném podprostoru variační, menší podprostor může energii jen zvýšit, takže neúplný podprostor je nejpřímějším vysvětlením rozdílu.

Tato kontrola přesnou FCI je možná jen proto, že je příklad malý; není součástí pracovního postupu SqDRIFT a pro velikosti systémů, kde je SQD skutečně potřeba, žádná taková reference neexistuje. Slouží zde pouze ke kvantifikaci toho, jak velkou část prostoru tento běh pokryl.

Srovnání s příkladem v Pythonu není kontrolované, takže by se rozdíl neměl připisovat pouze absenci obnovy konfigurací; všimni si, že počet obvodů, jediný čas vývoje, filtrování diagonálních členů i hardwarový šum se liší současně a všechny formují podprostor.

Abys rozdíl zmenšil, zde je několik možných dalších kroků k vyzkoušení po dokončení tohoto tutoriálu: zvyš shots, num_circuits nebo ops_per_circuit; vzorkuj ve více než jednom time_step; a přidej smyčku obnovy nad Qiskit::addon::sqd::recover_configurations.

Další kroky

Doporučení

• Sample-based Krylov quantum diagonalization of a fermionic lattice model — související tutoriál využívající obvody časového vývoje místo variačního ansatzu.

• Sample-based quantum diagonalization of a chemistry Hamiltonian — tutoriál o tom, jak sestavit obvod lokálního unitárního clusteru Jastrow (LUCJ) pro simulaci kvantové chemie.

• Článek SqDRIFT — literatura, na níž je tento tutoriál založen. (Všimni si, že některé optimalizace diskutované v tomto článku jsou zatím rozpracované a tento tutoriál se může v budoucnu měnit podle vývoje použitých knihoven.)