Kvantová diagonalizace na základě vzorků pro chemický Hamiltonián

Odhadovaná spotřeba: méně než jedna minuta na procesoru Heron r2 (POZNÁMKA: Jde pouze o odhad. Skutečná doba běhu se může lišit.)

Výsledky učení

• Jak použít SQD Qiskit addon k aproximaci energie základního stavu molekulárního systému pomocí bitových řetězců vzorkovaných z kvantové procesorové jednotky (QPU).

• Jak použít ffsim k sestavení Circuit lokálního unitárního shlukového Jastrow (LUCJ) pro simulaci kvantové chemie.

Předpoklady

• Kvantová chemie a druhá kvantizace

• Použití primitiva Sampler k vzorkování z kvantových obvodů

Pozadí

V tomto tutoriálu ukážeme, jak post-procesovat zašuměné kvantové vzorky k aproximaci základního stavu molekuly dusíku $\text{N}_2$ při rovnovážné délce vazby pomocí SQD Qiskit addonu k implementaci algoritmu kvantové diagonalizace na základě vzorků (SQD). Více podrobností o softwaru najdeš v příslušné dokumentaci, včetně jednoduchého příkladu pro začátek.

Chemie dusíku propojuje tento příklad s udržitelností prostřednictvím výroby amoniaku pro hnojiva. Pochopení toho, jak katalyzátory aktivují dusík, je součástí výzkumu energeticky méně náročných výrobních metod. Výpočet základního stavu izolované molekuly dusíku představuje metody elektronové struktury relevantní pro tuto širší výzvu.

Tento tutoriál je doporučen uživatelům, kteří jsou obeznámeni s kvantovou chemií: konkrétně s hledáním energií základního stavu molekuly. Podrobný průvodce pracovním postupem najdeš v kurzu algoritmu kvantové diagonalizace.

SQD je technika pro hledání vlastních hodnot a vlastních vektorů kvantových operátorů, jako je Hamiltonián kvantového systému, pomocí kombinace kvantového a distribuovaného klasického výpočtu. Klasický distribuovaný výpočet slouží ke zpracování vzorků získaných z kvantového procesoru a k projekci a diagonalizaci cílového Hamiltoniánu v podprostoru jimi generovaném. Pracovní postup SQD se skládá z následujících kroků:

1. Zvol Circuit ansatz a aplikuj ho na kvantovém počítači na referenční stav (v tomto případě stav Hartree-Fock).

2. Vzorkuj bitové řetězce z výsledného kvantového stavu.

3. Spusť proceduru self-consistent configuration recovery na bitových řetězcích, aby ses dostal k aproximaci základního stavu.

Je známo, že SQD funguje dobře, když je cílový vlastní stav řídký: vlnová funkce je podpořena na množině bázových stavů $\mathcal{S} = \{|x\rangle \}$, jejíž velikost neroste exponenciálně s velikostí problému.

Kvantová chemie

Hamiltonián molekulárního systému lze zapsat jako

$$\hat{H} = \sum_{ \substack{pr\\\sigma} } h_{pr} \, \hat{a}^\dagger_{p\sigma} \hat{a}_{r\sigma} + \frac12 \sum_{ \substack{prqs\\\sigma\tau} } h_{prqs} \, \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma},$$

kde $h_{pr}$ a $h_{prqs}$ jsou komplexní čísla nazývaná molekulární integrály, která lze vypočítat ze specifikace molekuly pomocí počítačového programu. V tomto tutoriálu počítáme integrály pomocí softwarového balíku PySCF.

Podrobnosti o tom, jak se odvozuje molekulární Hamiltonián, najdeš v učebnici kvantové chemie (například Modern Quantum Chemistry od Szaba a Ostlunda). Přehledné vysvětlení, jak se problémy kvantové chemie mapují na kvantové počítače, nabízí přednáška Mapping Problems to qubits z Qiskit Global Summer School 2024.

Ansatz lokálního unitárního shlukového Jastrow (LUCJ)

SQD vyžaduje ansatz kvantového obvodu, z něhož se vzorkuje. V tomto tutoriálu použijeme ansatz lokálního unitárního shlukového Jastrow (LUCJ) díky jeho kombinaci fyzikální motivace a přátelskosti k hardwaru. K sestavení Circuit ansatzu použijeme ffsim.

Ansatz LUCJ se přizpůsobuje QPU s omezeným propojením Qubitů. Spinové orbitaly jsou mapovány na Qubity tak, aby ansatz nevyžadoval routing pomocí SWAP Gate. Hardware IBM® má topologii Qubitů buď čtvercové mřížky, nebo heavy-hex mřížky. Pro heavy-hex mřížku můžeme použít vzor „cik-cak", znázorněný níže. V tomto vzoru jsou orbitaly stejného spinu mapovány na Qubity s topologií přímky (červené a modré kroužky) a spojení mezi orbitaly různého spinu je u každého 4. prostorového orbitalu, přičemž spojení je zprostředkováno ancilárním Qubitem (fialové kroužky).

Samokonzistentní obnova konfigurace

Procedura samokonzistentní obnovy konfigurací je navržena k extrakci co největšího signálu ze zašuměných kvantových vzorků. Protože molekulární Hamiltonián zachovává počet částic a spin Z, dává smysl volit Circuit ansatz, který tyto symetrie také zachovává. Při aplikaci na stav Hartree-Fock má výsledný stav v bezhlučném nastavení pevný počet částic a spin Z. Proto by spin-$\alpha$ a spin-$\beta$ části každého bitového řetězce vzorkovaného z tohoto stavu měly mít stejnou Hammingovu váhu jako ve stavu Hartree-Fock. Kvůli přítomnosti šumu v současných kvantových procesorech budou některé naměřené bitové řetězce tuto vlastnost porušovat. Jednoduchá forma postselekce by tyto bitové řetězce zahodila, ale to je plýtvání, protože mohou stále obsahovat určitý signál. Procedura samokonzistentní obnovy se pokouší část tohoto signálu získat zpět v post-processingu. Procedura je iterativní a jako vstup vyžaduje odhad průměrného obsazení každého orbitalu v základním stavu, který se nejprve vypočítá ze surových vzorků. Procedura běží ve smyčce a každá iterace má tyto kroky:

1. Pro každý bitový řetězec, který porušuje zadané symetrie, překlop jeho bity pravděpodobnostní procedurou navrženou tak, aby bitový řetězec přiblížila k aktuálnímu odhadu průměrného obsazení orbitalů, a získej nový bitový řetězec.

2. Shromáždi všechny staré a nové bitové řetězce, které splňují symetrie, a rozděl jejich podvýběr do disjunktních dávek pevné velikosti, zvolené předem.

3. Pro každou dávku bitových řetězců promítni Hamiltonián do podprostoru generovaného odpovídajícími bázovými vektory (popis těchto bázových vektorů viz předchozí sekci) a vypočítej odhad základního stavu promítnutého Hamiltoniánu na klasickém počítači. Každý bitový řetězec se skládá z části se spinem $\alpha$ a části se spinem $\beta$, z nichž každou nazýváme konfigurační interakční (CI) řetězec. Podprostor dávky je tvořen bázovými vektory, které páruje každý CI řetězec se spinem $\alpha$ v dávce s každým CI řetězcem se spinem $\beta$ v dávce.

4. Ořízni každou dávku tak, aby zůstaly pouze CI řetězce s největší váhou v odhadu základního stavu dané dávky. Sluč zbylé CI řetězce ze všech dávek a diagonalizuj Hamiltonián znovu ve sloučeném podprostoru.

5. Aktualizuj odhad průměrného obsazení orbitalů pomocí odhadu základního stavu ze sloučeného podprostoru. Také se přenášejí do podprostorů další iterace CI řetězce s nejvyšší váhou ze sloučeného podprostoru.

Operace ořezu a slučování v kroku 4 tvoří variantu algoritmu trim SQD. Diagonalizace jedné dávky stojí méně než diagonalizace podprostoru, který obsahuje všechny kandidátní CI řetězce. Proto může SQD s ořezáváním vybírat z většího souboru vzorkovaných CI řetězců, než dokáže diagonalizovat najednou.

Diagram pracovního postupu SQD

Následující diagram znázorňuje pracovní postup SQD. Neukazuje operace ořezu a slučování z kroku 4.

Požadavky

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

• Qiskit SDK v1.0 nebo novější, s podporou vizualizace

• Qiskit Runtime v0.22 nebo novější (pip install qiskit-ibm-runtime)

• SQD Qiskit addon v0.14 nebo novější (pip install qiskit-addon-sqd)

• ffsim v0.0.75 nebo novější (pip install ffsim)

Nastavení

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"ffsim": "ffsim", "matplotlib": "matplotlib", "numpy": "numpy", "pyscf": "pyscf", "qiskit": "qiskit", "qiskit_addon_sqd": "qiskit-addon-sqd", "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")

import math

import ffsim
import matplotlib.pyplot as plt
import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.primitives import StatevectorSampler
from qiskit.providers.fake_provider import GenericBackendV2
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

Příklad na malém simulátoru

V tomto tutoriálu najdeme aproximaci základního stavu molekuly dusíku v blízkosti její rovnovážné délky vazby. Nejprve použijeme malou bázi STO-6G, abychom mohli simulovat experiment a ověřit jeho správnost.

Krok 1: Mapování klasických vstupů na kvantový problém

Nejprve specifikujeme molekulu a její vlastnosti.

# Specify molecule properties
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, -0.5)], ["N", (0, 0, 0.5)]],
    basis="sto-6g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)

# Compute exact energy using FCI
reference_energy = cas.run().e_tot

print(f"norb = {norb}")
print(f"nelec = {nelec}")

converged SCF energy = -108.464957764795

CASCI E = -108.595987350986  E(CI) = -32.4115475088426  S^2 = 0.0000000

norb = 8
nelec = (5, 5)

Před sestavením Circuit ansatzu LUCJ nejprve provedeme výpočet CCSD v následující buňce kódu. Amplitudy $t_1$ a $t_2$ z tohoto výpočtu budou použity k inicializaci parametrů ansatzu.

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

E(CCSD) = -108.5933309085007  E_corr = -0.1283731437052356

Nyní použijeme ffsim k vytvoření Circuit ansatzu. Protože naše molekula má Hartree-Fockův stav s uzavřenými slupkami, použijeme spinově vyváženou variantu ansatzu UCJ, UCJOpSpinBalanced. V metodě from_t_amplitudes nastavíme optimize=True, aby bylo povoleno „komprimované" dvojité faktorizování amplitud $t_2$ (viz The local unitary cluster Jastrow (LUCJ) ansatz v dokumentaci ffsim pro podrobnosti).

Protože se ansatz LUCJ přizpůsobuje dostupnému propojení QPU, musíme před vytvořením ansatzu inicializovat backend QPU. Prozatím vytvoříme generický backend s coupling mapou čtvercové mřížky a sadou Gate, do které se ansatz LUCJ přirozeně rozkládá. Poté použijeme ffsim.qiskit.generate_lucj_pass_manager k vytvoření pass manageru specializovaného pro transpilaci ansatzu LUCJ na daný backend podle rozložení „žebřík" (viz vysvětlení v dokumentaci ffsim). Tato funkce používá heuristiku skóre k minimalizaci chyb spojených s vybraným rozložením, což je důležité, pokud je tvůj backend skutečnou QPU nebo simulátorem s modelem šumu. Kromě vrácení pass manageru tato funkce vrací také páry vazeb alpha-beta, které lze na hardwaru implementovat. Pokud nelze implementovat všechny páry, vydá varování.

import warnings

from qiskit.transpiler import CouplingMap

warnings.formatwarning = lambda msg, *args, **kwargs: f"Warning: {msg}\n"

# Set ansatz properties
n_reps = 4
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
rng = np.random.default_rng(12345)
coupling_map = CouplingMap.from_grid(norb, 2)
backend = GenericBackendV2(
    coupling_map.size(),
    coupling_map=coupling_map,
    basis_gates=["cp", "xx_plus_yy", "p", "x", "swap"],
    seed=rng,
)

# Create pass manager
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity="square",
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
    seed_transpiler=int(rng.integers(2**31)),
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it
# to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

Krok 2: Optimalizace pro spuštění na kvantovém hardwaru

Dále optimalizujeme Circuit pro cílový hardware. Obvykle tento krok zahrnuje inicializaci hardwarového backendu a pass manageru pro daný backend. Protože se však ansatz LUCJ přizpůsobuje propojení hardwaru, tyto akce jsme již provedli v předchozím kroku. Zbývá pouze spustit pass manager na obvod a transpilovat ho na ISA Circuit, který lze přímo spustit na QPU.

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

Gate counts: OrderedDict({'xx_plus_yy': 254, 'cp': 88, 'p': 64, 'measure': 16, 'x': 10, 'barrier': 1})

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

Po optimalizaci Circuit pro spuštění na hardwaru jsme připraveni spustit ho na cílovém hardwaru a shromáždit vzorky pro odhad energie základního stavu. Protože máme pouze jeden obvod, použijeme IBM Quantum Compute Service režim spuštění úlohy a spustíme náš obvod.

sampler = StatevectorSampler(seed=rng)
job = sampler.run([isa_circuit], shots=100_000)

primitive_result = job.result()
pub_result = primitive_result[0]

Krok 4: Post-processing a vrácení výsledku v požadovaném klasickém formátu

Užitečnou metrikou pro posouzení kvality výstupu QPU je počet vrácených platných konfigurací. Platná konfigurace má správný počet částic a spin Z, což znamená, že pravá polovina bitového řetězce má Hammingovu váhu rovnou počtu elektronů se spinem nahoru a levá polovina má Hammingovu váhu rovnou počtu elektronů se spinem dolů. Následující buňka vypočítá podíl vzorkovaných konfigurací, které jsou platné.

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
    )

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(
    f"Fraction of sampled configurations that are valid (simulation): {valid_fraction}"
)

Fraction of sampled configurations that are valid (simulation): 1.0

Všechny bitové řetězce jsou platné, protože vzorkujeme Circuit na bezhlučném simulátoru. Při spuštění na zašuměném QPU bude tento podíl nižší než jedna, ale doufejme, že bude vyšší než podíl, který bychom očekávali, kdyby byly bitové řetězce vzorkovány rovnoměrně náhodně, což je vypočítáno v následující buňce.

expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)

Expected fraction of valid configurations from uniformly random bitstrings: 0.0478515625

Nyní odhadneme energii základního stavu Hamiltoniánu pomocí funkce diagonalize_fermionic_hamiltonian. Tato funkce provádí postup self-konzistentního obnovení konfigurací, který iterativně zpřesňuje zašuměné kvantové vzorky za účelem zlepšení odhadu energie. Předáváme callback funkci, abychom mohli uložit mezivýsledky pro pozdější analýzu. Vysvětlení argumentů funkce diagonalize_fermionic_hamiltonian najdeš v dokumentaci API.

Zde použijeme argument initial_occupancies funkce diagonalize_fermionic_hamiltonian k zadání konfigurace Hartree-Fock jako počátečního odhadu obsazení orbitalů v základním stavu. Tento přístup je vhodný pro systémy, kde má základní stav významnou podporu na konfiguraci Hartree-Fock, ale nemusí být vhodný v jiných situacích, ačkoli pokročilejší výpočetní metody mohou v takových případech poskytnout lepší počáteční odhady. Zadání initial_occupancies také umožňuje spuštění obnovy konfigurace, i když nebyly vzorkovány žádné platné konfigurace, což může nastat při vzorkování velkého Circuit na zašuměné QPU. Bez tohoto argumentu by obnova konfigurace selhala a vyvolala chybu, pokud by nebyly poskytnuty žádné platné konfigurace.

Trim SQD zapneme předáním TrimPolicy jako argumentu policy. Každá iterace pak provede num_batches + 1 diagonalizací: jednu pro každou dávku a jednu pro sloučený podprostor. Používáme dvě dávky, takže každá iterace stojí tři diagonalizace. Callback dostává pouze sloučený výsledek, takže každá iterace vypíše jedinou energii. Možnosti TrimPolicy jsou:

• trim_ratio: Podíl CI řetězců každé dávky (na spinový sektor), které ořezem přežijí a vstoupí do sloučeného podprostoru. Sloučený podprostor obsahuje až zhruba num_batches * trim_ratio krát tolik CI řetězců jako jedna dávka, bez těch, které mají dávky společné. Argument max_dim omezuje každý podprostor, včetně sloučeného.

• carryover_threshold a carryover_ratio: Tyto volby určují, které CI řetězce sloučeného podprostoru se přenesou do další iterace. Aby se CI řetězec přenesl, musí se objevit v konfiguraci, jejíž amplituda překračuje carryover_threshold. Z těchto řetězců se ponechá podíl carryover_ratio s největší váhou.

Tento systém je dost malý na to, aby jediná dávka obsahovala většinu možných CI řetězců, takže se dávky silně překrývají. Proto používáme vysoké trim_ratio, aby sloučený podprostor nebyl mnohem menší než dávky. Další podrobnosti o volbě těchto možností najdeš v příručce k trim SQD.

from functools import partial

from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
    solve_sci_batch,
)
from qiskit_addon_sqd.trim import TrimPolicy

# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5
samples_per_batch = 1000
num_batches = 2
max_dim = 50
policy = TrimPolicy(
    trim_ratio=0.8,
    carryover_threshold=1e-4,
    carryover_ratio=0.5,
)
symmetrize_spin = True

# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the sci_solver argument
# in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=200)

# List to capture intermediate results
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"\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    policy=policy,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    max_dim=max_dim,
    initial_occupancies=initial_occupancies,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy (simulation): {final_energy}")
print(f"Final energy error (simulation): {energy_error}")

Iteration 1
	Energy: -108.59585971421083
	Subspace dimension: 1681

Iteration 2
	Energy: -108.59585971421083
	Subspace dimension: 1681
Final energy (simulation): -108.59585971421083
Final energy error (simulation): 0.00012763677520410965

Vizualizace výsledků

Pro každou iteraci ukazuje první graf chybu energie vypočtené ve sloučeném podprostoru. V této simulaci jsme již po první iteraci v rozmezí 1 mH od přesné odpovědi (chemická přesnost se obvykle považuje za 1 kcal/mol $\approx$ 1.6 mH). Jde ale o malý systém a protože jsou vzorky bezšumové, konfigurační obnova (configuration recovery) není potřeba. U většího systému spuštěného na zašuměném QPU by mohlo být potřeba více iterací konfigurační obnovy a výsledná přesnost by mohla být horší. Obecně lze energii zlepšit povolením více iterací konfigurační obnovy nebo zvýšením počtu vzorků na dávku.

Druhý graf zobrazuje průměrnou obsazenost každého prostorového orbitalu po poslední iteraci. Vidíme, že elektrony se spinem nahoru i se spinem dolů obsazují s vysokou pravděpodobností prvních pět orbitalů v našich řešeních.

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]

# Chemical accuracy (+/- 1.6 milli-Hartree)
chem_accuracy = 0.0016

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-5, 1e-2)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

plt.tight_layout()
plt.show()

Příklad na velkém hardwaru

Nyní spustíme větší příklad na skutečném kvantovém hardwaru. Zde odvodíme aktivní prostor pro molekulu dusíku z báze cc-pVDZ. Provedeme také následující změny:

• Snížíme n_reps na 1, abychom omezili šum způsobený chybami gate.
• Nastavíme propojení ansatzu podle rodiny procesorů backendu (heavy-hex pro Heron, čtvercová pro Nighthawk).
• Zvýšíme max_dim na 1000, abychom umožnili větší dimenze podprostoru.
• Snížíme trim_ratio na 0.5. Tento systém má mnohem více možných CI řetězců, než se vejde do podprostoru. Při dvou dávkách poměr 0.5 zajistí, že sloučený podprostor nebude větší než jedna dávka.

Kroky 1–4

Zde shrneme všechny kroky do jednoho pracovního postupu ve větším měřítku, který je poté spuštěn na skutečném kvantovém hardwaru.

# ------------------------------ Step 1 ------------------------------
# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, -0.5)], ["N", (0, 0, 0.5)]],
    basis="cc-pvdz",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
norb = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
n_alpha = (n_electrons + mol.spin) // 2
n_beta = (n_electrons - mol.spin) // 2
nelec = (n_alpha, n_beta)
cas = pyscf.mcscf.CASCI(scf, norb, nelec)
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)

# Store reference energy from SCI calculation performed separately
reference_energy = -109.22802921665716

print(f"norb = {norb}")
print(f"nelec = {nelec}")

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
).run()
t1 = ccsd.t1
t2 = ccsd.t2

# Set ansatz properties
n_reps = 1
pairs_aa = [(p, p + 1) for p in range(norb - 1)]

# Let generate_lucj_pass_manager determine the alpha-beta interactions
pairs_ab = None

# Initialize backend
service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
print(f"Using backend {backend.name}")

# Create pass manager
connectivity_map = {"Heron": "heavy-hex", "Nighthawk": "square"}
pass_manager, pairs_ab = ffsim.qiskit.generate_lucj_pass_manager(
    backend=backend,
    norb=norb,
    connectivity=connectivity_map[
        backend.configuration().processor_type["family"]
    ],
    interaction_pairs=(pairs_aa, pairs_ab),
    optimization_level=3,
)

# Create the LUCJ ansatz operator
ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(pairs_aa, pairs_ab),
    # Setting optimize=True enables the "compressed" factorization
    optimize=True,
)

# create an empty quantum circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it
# to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()

# ------------------------------ Step 2 ------------------------------

isa_circuit = pass_manager.run(circuit)
print(f"Gate counts: {isa_circuit.count_ops()}")

# ------------------------------ Step 3 ------------------------------
sampler = Sampler(mode=backend)
sampler.options.environment.job_tags = ["TUT_SQD"]
job = sampler.run([isa_circuit], shots=100_000)
primitive_result = job.result()
pub_result = primitive_result[0]

# ------------------------------ Step 4 ------------------------------

bit_array = pub_result.data.meas
num_valid = sum(
    is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()
)
valid_fraction = num_valid / bit_array.num_shots
print(
    f"Fraction of sampled configurations that are valid (hardware): {valid_fraction}"
)
expected_fraction_random = (
    math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)
)
print(
    f"Expected fraction of valid configurations from uniformly random bitstrings: "
    f"{expected_fraction_random}"
)
# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5
samples_per_batch = 1000
num_batches = 2
max_dim = 1000
policy = TrimPolicy(
    trim_ratio=0.5,
    carryover_threshold=1e-4,
    carryover_ratio=0.5,
)
symmetrize_spin = True

# Use the Hartree-Fock configuration as an initial guess for the
# orbital occupancies
initial_occupancies = (
    np.array([1] * n_alpha + [0] * (norb - n_alpha)),
    np.array([1] * n_beta + [0] * (norb - n_beta)),
)

# Pass options to the built-in eigensolver. If you just want to use the defaults,
# you can omit this step, in which case you would not specify the
# sci_solver argument in the call to diagonalize_fermionic_hamiltonian below.
sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=200)

# List to capture intermediate results
result_history = []

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=norb,
    nelec=nelec,
    num_batches=num_batches,
    policy=policy,
    energy_tol=energy_tol,
    occupancies_tol=occupancies_tol,
    max_iterations=max_iterations,
    sci_solver=sci_solver,
    symmetrize_spin=symmetrize_spin,
    max_dim=max_dim,
    initial_occupancies=initial_occupancies,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy (hardware): {final_energy}")
print(f"Final energy error (hardware): {energy_error}")

# Data for energies plot
x1 = range(len(result_history))
min_e = [
    min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy
    for result in result_history
]
e_diff = [abs(e - reference_energy) for e in min_e]

# Chemical accuracy (+/- 1.6 milli-Hartree)
chem_accuracy = 0.0016

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-3, 1e-1)
axs[0].axhline(
    y=chem_accuracy,
    color="#BF5700",
    linestyle="--",
    label="chemical accuracy",
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

plt.tight_layout()
plt.show()

converged SCF energy = -108.929838385609
norb = 26
nelec = (5, 5)
E(CCSD) = -109.2177884185544  E_corr = -0.2879500329450045
Using backend ibm_kingston

Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20), (24, 24)].
Removing interaction (24, 24) from the end.
Warning: Backend cannot accommodate pairs_ab=[(0, 0), (4, 4), (8, 8), (12, 12), (16, 16), (20, 20)].
Removing interaction (20, 20) from the end.

Gate counts: OrderedDict({'sx': 7039, 'rz': 6977, 'cz': 1858, 'measure': 52, 'x': 45, 'barrier': 1})
Fraction of sampled configurations that are valid (hardware): 0.00248
Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07
Iteration 1
	Energy: -109.19178247530117
	Subspace dimension: 644809
Iteration 2
	Energy: -109.20055015645298
	Subspace dimension: 339889
Iteration 3
	Energy: -109.20621939005301
	Subspace dimension: 446224
Iteration 4
	Energy: -109.21228961696536
	Subspace dimension: 393129
Iteration 5
	Energy: -109.21328434918014
	Subspace dimension: 403225
Final energy (hardware): -109.21328434918014
Final energy error (hardware): 0.014744867477020307

Další kroky

Doporučení

Pokud tě tato práce zaujala, možná tě budou zajímat i následující materiály:

• Kvantová diagonalizace na základě vzorků metodou Krylov pro fermionový mřížkový model – příbuzný tutoriál využívající obvody s časovým vývojem namísto variačního ansatzu
• Škálování pracovních postupů SQD pro chemii pomocí solveru Dice – stránka ukazující, jak používat efektivnější software Dice pro diagonalizaci
• Dokumentace API SQD addonů – reference pro funkci diagonalize_fermionic_hamiltonian
• Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer – článek, na němž je tento tutoriál založen