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 hluč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 addonů 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 se vyskytuje každé 4. prostorové orbitaly, přičemž spojení je zprostředkováno ancilárním Qubitem (fialové kroužky).

Samokonzistentní obnova konfigurace

Procedura self-consistent configuration recovery je navržena k extrakci co největšího signálu z hluč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 postselektion by tyto bitové řetězce zahodila, ale to je plýtvání, protože mohou stále obsahovat určitý signál. Procedura self-consistent recovery 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 vzorkuj podmnožiny pevné velikosti, zvolené předem.

3. Pro každou podmnožinu 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.

4. Aktualizuj odhad průměrného obsazení orbitalů pomocí odhadu základního stavu s nejnižší energií.

Diagram pracovního postupu SQD

Pracovní postup SQD je znázorněn v následujícím diagramu:

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.12 nebo novější (pip install qiskit-addon-sqd)

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

Nastavení

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

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.464957764796
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.5933309085008  E_corr = -0.1283731437052352

Nyní použijeme ffsim k vytvoření Circuit ansatzu. Protože naše molekula má uzavřený stav Hartree-Fock, použijeme spin-vyvážený variant 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
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"],
)

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

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

rng = np.random.default_rng(12345)
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: {valid_fraction}")

Fraction of sampled configurations that are valid: 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 hlučné 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 hlučné QPU. Bez tohoto argumentu by obnova konfigurace selhala a vyvolala chybu, pokud by nebyly poskytnuty žádné platné konfigurace.

from functools import partial

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

# SQD options
energy_tol = 1e-3
occupancies_tol = 1e-3
max_iterations = 5

# Eigenstate solver options
num_batches = 3
samples_per_batch = 1000
max_dim = 50
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# 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=max_cycle)

# 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"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy + nuclear_repulsion_energy}")
        print(
            f"\t\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,
    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,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

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

Iteration 1
	Subsample 0
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 1
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 2
		Energy: -108.59589540253083
		Subspace dimension: 2500
Iteration 2
	Subsample 0
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 1
		Energy: -108.59589540253083
		Subspace dimension: 2500
	Subsample 2
		Energy: -108.59589540253083
		Subspace dimension: 2500
Final energy: -108.59589540253083
Final energy error: 9.19484552071026e-05

Vizualizace výsledků

První graf ukazuje, že v této simulaci jsme po první iteraci již 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 vylepšit povolením více iterací konfigurační obnovy nebo zvýšením počtu vzorků na batch.

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 ze základní sady 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 na základě rodiny procesorů backendu (heavy-hex pro Heron, čtvercová pro Nighthawk)
• Zvýšíme max_dim na 1000, abychom umožnili větší dimenze podprostoru.

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: {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

# Eigenstate solver options
num_batches = 3
samples_per_batch = 1000
max_dim = 1000
symmetrize_spin = True
carryover_threshold = 1e-4
max_cycle = 200

# 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=max_cycle)

# 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,
    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,
    carryover_threshold=carryover_threshold,
    callback=callback,
    seed=rng,
)

final_energy = result.energy + nuclear_repulsion_energy
energy_error = final_energy - reference_energy
print(f"Final energy: {final_energy}")
print(f"Final energy error: {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.2879500329450041
Using backend ibm_boston

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': 7036, 'rz': 6965, 'cz': 1858, 'measure': 52, 'x': 50, 'barrier': 1})
Fraction of sampled configurations that are valid: 0.01677
Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07
Iteration 1
	Subsample 0
		Energy: -109.17922180754668
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.17713541195594
		Subspace dimension: 962361
	Subsample 2
		Energy: -109.17407767345108
		Subspace dimension: 968256
Iteration 2
	Subsample 0
		Energy: -109.1873108379026
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.1876970140745
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.18255107152576
		Subspace dimension: 1000000
Iteration 3
	Subsample 0
		Energy: -109.19052148187927
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19316375313544
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19087348337145
		Subspace dimension: 1000000
Iteration 4
	Subsample 0
		Energy: -109.19374847699348
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19517039328333
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19422416032815
		Subspace dimension: 1000000
Iteration 5
	Subsample 0
		Energy: -109.19580768221108
		Subspace dimension: 1000000
	Subsample 1
		Energy: -109.19608553581007
		Subspace dimension: 1000000
	Subsample 2
		Energy: -109.19548131833729
		Subspace dimension: 1000000
Final energy: -109.19608553581007
Final energy error: 0.03194368084709254

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í Circuit 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