Algoritmus SqDRIFT pro odhad základního stavu

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

Hledáte verzi v C++?

Tento tutoriál používá Python. Pro implementaci v C++, včetně zdrojového kódu a pokynů pro sestavení, viz tutoriál C++ SqDRIFT.

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

Tento tutoriál je pro výukové účely prezentován jako Python notebook.

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.

Tento tutoriál využívá Qiskit Fermions k vytvoření přirozenějších fermionických obvodů pro algoritmus qDRIFT, následované použitím fermionických passů pro rozvržení a syntézu před zapojením obvodů do tradičního pipeline Qiskitu pro hardwarové vykonání.

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í

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

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

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

# Qiskit Aer
from qiskit_aer import AerSimulator

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

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

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

Příklad se simulátorem

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

Č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 symetrický štítek 1, tedy nevyužívá se žádná bodová grupová symetrie. 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 další buňce.

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, se znovu vygenerovaný soubor 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é buňky 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 = "assets/sqdrift/fcidump_files/N2_sto_3g"

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

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

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é excitovaný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 excitované č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$$

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 pohodlím pro tutoriál, který by měl 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á doprovodný C++ kód, jelikož ten místo obnovy provádí postselekci a liší se i v dalších parametrech.

name = "assets/sqdrift/fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

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

print(len(canon.groups))

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

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

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

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

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

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

for circ in sqdrift_circuits:
    circ.measure_all()

print(len(sqdrift_circuits))

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

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

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

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

shots = 100

sampler = Sampler(mode=backend)

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

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

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

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

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

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

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

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

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

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

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

result_history = []

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

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

computed_energy = result.energy + nuclear_repulsion_energy

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

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

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

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

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

Další kroky

Doporučení

Pokud tě tato práce zaujala, mohly by tě zajímat následující materiály:

• Krylovova kvantová diagonalizace založená na vzorkování pro fermionový mřížkový model - související tutoriál používající obvody časového vývoje místo variačního ansatzu.
• Kvantová diagonalizace založená na vzorkování pro chemický hamiltonián - tutoriál o tom, jak sestavit obvod local unitary cluster Jastrow (LUCJ) pro simulaci kvantové chemie.
• Článek SqDRIFT - literatura, na které je tento tutoriál založen. (Všimni si, že některé optimalizace probírané v tomto článku jsou v současnosti ještě rozpracované a tento tutoriál se může v budoucnu měnit podle vývoje použitých knihoven.)