Algoritmus SqDRIFT pro odhad základního stavu

Odhad využití: 180 sekund na procesoru Heron r3 (POZNÁMKA: Jde pouze o odhad. Tvoje doba běhu se může lišit.)

Výukové cíle

• Nauč se vytvářet obvody s menší hloubkou ve srovnání s Trotterizací

• Projdi kompletním pracovním postupem odhadu základního stavu pomocí qDRIFT a SQD

• Nauč se používat qiskit-fermions společně s dalšími doplňky Qiskit k implementaci takového pracovního postupu

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

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 nutnost volit ansatz, z něhož se vzorkují bitové řetězce, souborem časově evolučních obvodů sestavených přímo z cílového hamiltoniánu. Dosahuje se toho podvzorkováním menších časově evolučních operátorů z hamiltoniánu podle 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 fermionových obvodů pro algoritmus qDRIFT, následovanému použitím fermionového layoutu a syntézních passů před zapojením obvodů do tradičního pipeline Qiskit pro spuštění na hardwaru.

Nechť má hamiltonián tvar:

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

kde bez újmy na obecnosti požadujeme $c_i > 0$ a aby největší vlastní hodnota $h_i$ byla v absolutní hodnotě rovna $1$. Jakýkoli znaménkový nebo komplexní předfaktor je pohlcen do $h_i$, takže koeficienty $c_i$ jsou ryze 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ů navzorkovaných do jednoho obvodu, který je níže značen $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_{th}$ obvod SqDRIFT, definovaný jako:

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

Zde $n$ je počet navzorkovaných operátorů na obvod a $K$ je počet obvodů v souboru. Součin probíhá přes $n$ tahů, ne přes všech $N$ členů hamiltoniánu, a protože se členy losují s vracením, může se stejné $h_i$ objevit ve jediném $V_k$ vícekrát.

Veličina:

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

je norma $L_1$ koeficientů, takže každý z $n$ kroků se vyvíjí po stejnou dobu $\lambda t / n$ bez ohledu na to, který člen byl vylosován. Jednotnost úhlu kroku je charakteristickým rysem qDRIFT: koeficient ovlivňuje výsledek tím, jak často je jeho člen losován, ne tím, jak daleko je tento člen otočen. Indexy se vzorkují z rozdělení:

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

takže řada $(k_1, \ldots, k_n)$ je náhodná posloupnost indexů členů vylosovaných z tohoto rozdělení. Protože $c_i$ jsou kladná a jejich součet je $\lambda$, jde o normalizované rozdělení pravděpodobnosti a střední hodnota výsledného kanálu přes náhodné tahy aproximuje evoluci pod $H$ s chybou, která klesá s rostoucím $n$. Všimni si, že chyba aproximace závisí na $\lambda$, a ne na počtu členů $N$.

(Článek o SqDRIFT značí počet členů jako $\mathcal{N}$ a délku posloupnosti jako $N$; my zde používáme $N$ a $n$, aby byla obě jasně odlišena.)

Tento tutoriál ukazuje, jak vygenerovat soubor takových randomizovaných obvodů. Po vytvoření těchto obvodů, podobně jako 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. Tím je zajištěn vyšší překryv mezi vektory základního stavu a navzorkovaný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í Pythonu (>=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 potřebné balíčky nainstaluješ 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 — installs the packages this notebook needs if they are missing
import importlib.util

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

# Third-party scientific computing
import numpy as np

# PySCF
from pyscf import tools, ao2mo, fci

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

# Qiskit Aer
from qiskit_aer import AerSimulator

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

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

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

Příklad se simulátorem

Krok 1: Převod klasických vstupů na kvantový problém

Načtení a příprava FCIDump

V tomto tutoriálu načteme hamiltonián elektronové struktury dusíku (N2). Fermionové operátory lze vytvářet i jinými způsoby. Viz dokumentace 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 meziatomové 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ů při Jordan-Wignerově transformaci), 14 elektronů ve spinovém singletu, tedy sedm elektronů $\alpha$ a sedm $\beta$. Všem orbitalům je přiřazen symetrický štítek 1, tedy se nevyužívá žádná bodová grupa symetrie. Protože jde o výpis STO-3G v plném prostoru, nejsou žádné orbitaly zmrazeny a korelační prostor je dost malý na to, aby bylo možné klasicky spočítat přesnou referenční energii FCI pro srovnání, jak ukazuje další buňka.

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ávisejí na zkonvergovaných orbitalech SCF, může se znovu vygenerovaný soubor lišit od dodaného fází nebo pořadím orbitalů; celkové energie to neovlivní.

Získání souboru. FCIDump najdeš v tomto repozitáři GitHub. Spuštěním buňky níže ho stáhneš do umístění, které zbytek tutoriálu očekává.

Nejprve použijeme cisolver z pyscf k získání referenční energie. To je skutečná energie základního stavu molekuly, se kterou pracujeme. K tomu nejprve deklarujeme norb a nelec, což jsou po řadě počet orbitalů a počet elektronů. Poté deklarujeme h1e a h2e, což jsou po řadě jednoelektronové a dvouelektronové integrály. Všechny tyto veličiny později použijeme 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}")

Using existing FCIDump at assets/sqdrift/fcidump_files/N2_sto_3g

name = "fcidump_files/N2_sto_3g"

fcidump = tools.fcidump.read(name)

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

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

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

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

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

reference_energy = e_fci

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

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

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

Načtení hamiltoniánu

Když máme potřebná data připravená, 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 zobrazíme hamiltonián do fermionového modelu obvodu pomocí qiskit-fermions, který poskytuje passy transpilátoru a hradla specifická pro fermionové obvody. Ty se později použijí před tradičními passy transpilátoru Qiskit v tomto pracovním postupu.

Seskupení členů

Aby byly výsledky reprodukovatelné, nejprve pomocí canonical_order seřadíme členy pouze podle jejich struktury. Pořadí operátorů v seznamu canon je tedy pevné. To zajišťuje reprodukovatelnost vytvořených operátorů, protože pass QDriftTrotterization, který budeme používat dále, vzorkuje náhodné indexy k vytvoření operátorů qDRIFT.

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

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

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ů

Z hamiltoniánu použitého ke generování obvodů odstraníme diagonální členy, aby se $n$ vzorkovacích slotů qDRIFT využilo na členy, které přesouvají populaci mezi konfiguracemi. Takové členy je nejlepší odfiltrovat z hamiltoniánu právě teď, než se v dalším kroku sestaví hradlo Evolution.

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

• konstantní energetický posun, součin nula operátorů počtu, jehož časový vývoj přispívá jen globální fází;

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

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

Samy o sobě žádný z nich nepřesouvá populaci mezi konfiguracemi obsazovacích čísel; působí jen na fáze konfigurací, které už jsou přítomny. Nejsou však netečné: tyto relativní fáze vstupují do interference vytvářené excitačními členy později v obvodu, takže jejich odfiltrování mění evoluci, která se skutečně generuje, a může změnit vzorkovací rozdělení. Jde o záměrnou aproximaci v kroku generování obvodů, provedenou proto, aby se vzorkování soustředilo na excitační členy, ne o krok, který vzorkované rozdělení ponechává nedotčené. Na rozdíl od seskupení podle symetrie výše, které záruky konvergence qDRIFT zachovává, tento filtr mění operátor, který se vyvíjí. Obvody tedy už neaproximují evoluci pod plným hamiltoniánem a meze chyby qDRIFT platí pro filtrovaný operátor, ne pro původní. Zde je to přijatelné, protože obvody jsou jen vzorkovací heuristikou sloužící k navržení konfigurací: z odhadu energie samotného se neztrácí žádný člen, protože filtr se týká jen hamiltoniánu použitého 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í v navzorkovaném podprostoru bez ohledu na to, jak byly konfigurace navrženy.

Funkce filter_diagonal_terms() takové členy z operátoru odstraní přímo na místě. Rozpozná je podle jejich struktury v normálním uspořádání — multimnožina kreačních módů odpovídá multimnožině anihilačních módů — takže je platná jen na operátoru, který už je v normálním uspořádání. 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ž máme členy hamiltoniánu seskupené, rozhodneme se o následujících parametrech pro generování souboru obvodů:

• Počet obvodů ke generování: num_circuits
• Délka každého obvodu v počtu excitačních skupin: num_exc
• Faktor pro různé evoluční časy: times

Vytvoř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 evolučním časem, který jsme deklarovali dříve. Evolučním operátorem je hamiltonián. Později na těchto obvodech spustíme passy transpilátoru, abychom vytvořili obvody qDRIFT.

Příprava ansatzu

Hartreeho-Fockův stav připravíme pomocí třídy InitializeModes. Pro dusík jde jednoduše o 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 elektronů $\alpha$ a sedm $\beta$ 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 problému pro spuštění na kvantovém hardwaru

Nyní, když máme své obvody, nejprve použijeme passy dostupné v qiskit-fermions k provedení optimalizací na fermionové úrovni a poté transpilujeme obvod pro zvolený backend. Protože jde o experiment se simulátorem, nejprve to uděláme pro AerSimulator. Výpočet vah pro každou skupinu

V tomto kroku provedeme qDRIFT vzorkování členů stochasticky s pravděpodobnostmi úměrnými jejich koeficientům v hamiltoniánu. Udělá to za nás průchod transpilátoru qDRIFT. Nyní můžeme vytvořit mělčí obvody, které lze na hardwaru provádět efektivněji navzdory omezené konektivitě qubitů, a to i tehdy, když hamiltonián obsahuje dalekodosahové vazby a členy vyššího než kvadratického řádu. Po seskupení členů se operátory vzorkují podle jejich vah. Pro každý operátor $h_i$ je váha $W_{h_i}$ definována takto:

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

Protože byly členy seskupeny v kroku 1, je každé $h_i$ zde celá skupina: $c_i$ je střední absolutní hodnota 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 vytváří optimalizovaný finální obvod, který můžeme transpilovat pro běh na našem hardwaru. Jeho výchozí optimalizační fázi nahradíme FermionicPassManager obsahujícím náš průchod QDriftTrotterization:

• Průchod QDriftTrotterization interně používá výpočet vah a vzorkování k vytvoření 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 dočteš v referenci API

Zbývající fáze MultiStagePassManager se spouštějí automaticky a zajišťují kompletní 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 pro mapování instrukcí obvodu založených na fermionech 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

Když jsme nyní hotovi s optimalizacemi na fermionové úrovni, můžeme obvody transpilovat pro provedení na simulátoru.

simulator = AerSimulator()
shots = 100

transpiled_circuits = transpile(sqdrift_circuits, simulator)

Krok 3: Provedení pomocí primitiv Qiskit

Nyní, když máme své obvody, můžeme je spustit pomocí primitiv Qiskit na AerSimulator. Sloučíme všechny počty měření z různých obvodů. Převedeme je na booleovské vektory, než je nakonec zpracujeme pomocí SQD.

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

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

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

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

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

Krok 4: Následné zpracování a vrácení výsledku v požadovaném klasickém formátu

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

Nyní můžeme na vybraných bitových řetězcích spustit schéma diagonalizace, abychom našli nejnižší vlastní hodnotu, která odpovídá energii základního stavu molekuly. Vytvoříme callback funkci, deklarujeme počáteční obsazení a nastavíme parametry, než nakonec schéma diagonalizace spustíme. Callback funkce slouží k vypsání 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 k výsledné energii nuclear_repulsion_energy.

Poznámka: Dimenze podprostoru není napříč iteracemi pevná, a to ani na simulátoru bez šumu — každý podvýběr vytáhne jinou množinu konfigurací a krok obnovy mezi iteracemi přetváří zásobu, takže hlášená dimenze se mezi podvýběry liší. Vzorkování bez šumu samo o sobě dimenzi vybraného podprostoru nezafixuje. Běh na hardwaru však má tendenci dávat systematicky větší podprostory, protože zašuměná měření porušují symetrii zachování počtu částic a obnova konfigurací je mění na další bázové vektory. Proto v hardwarové části zavedeme také další krok pro prořezávání 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ý má běžet rychle, nikoli pevným stropem metody.

Náklady klasického kroku přímo nezávisejí na počtu qubitů. SQD diagonalizuje hamiltonián promítnutý do podprostoru nataženého vzorkovanými konfiguracemi, takže klasické náklady určuje dimenze tohoto vybraného podprostoru — zde ovlivněná hodnotami samples_per_batch, num_batches a tím, kolik různých konfigurací obvody skutečně vytvoří — spolu s řídkou lineární algebrou potřebnou k aplikaci promítnutého hamiltoniánu. Úplný CI prostor roste kombinatoricky s počtem orbitalů a elektronů, ale vybraný podprostor je jeho malý, nastavitelný výřez a jeho velikost přímo řídíme. Počet qubitů a klasická obtížnost se proto mohou měnit poněkud nezávisle: širší orbitalový prostor vzorkovaný do skromného podprostoru může být levnější než menší systém diagonalizovaný ve velmi velkém.

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 řešiči vlastních hodnot. Větší orbitalové prostory obvykle vyžadují větší podprostor k dosažení chemické přesnosti, a právě to nakonec motivuje použití distribuovaných zdrojů — viz qiskit-addon-sqd-hpc pro škálování tohoto kroku. Místo předpokládání pevné hranice je praktickým přístupem sledovat hlášenou dimenzi podprostoru a konvergenci energie napříč iteracemi a zvětšovat podprostor, dokud se energie přestane zlepšovat nebo nevyčerpáš dostupnou paměť.

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

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

Zde se můžeme rozhodnout provést další krok. Když máme všechny bitové řetězce z provedení obvodů, můžeme před spuštěním SQD buď odfiltrovat neplatné bitové řetězce, nebo pokračovat bez prořezávání. Vynechání prořezávání je u běhů na hardwaru obecně výhodnější, protože ponechá měření 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 tato měření rovnou zahodila.

Protože dusík může mít pouze sedm elektronů $\alpha$ a sedm elektronů $\beta$, lze zahodit všechny bitové řetězce, které mají v první a 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 falešné bitové řetězce odfiltrujeme, zbytek se pošle do schématu diagonalizace. Pomocí níže uvedeného příznaku PRUNE můžeš přepínat mezi oběma chováními.

Pamatuj, že prořezávání je jen jednou z několika voleb, které utvářejí výsledný podprostor, vedle počtu obvodů, množiny časů vývoje a filtrování diagonálních členů. Srovnání prořezaného běhu s neprořezaným je informativní jen tehdy, když je vše ostatní drženo beze změny; verze tohoto tutoriálu v C++ to rozebírá podrobněji, protože provádí postselekci místo obnovy a liší se i v těchto 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")

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, mohl by tě zajímat následující materiál:

• 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 klastrového Jastrowa (LUCJ) pro simulaci v kvantové chemii.
• Článek SqDRIFT - literatura, na které tento tutoriál staví. (Upozornění: 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.)