Přeskočit na hlavní obsah

Simulace neutronového rozptylu pomocí Serverless workflow AQC + Trotter dynamics

Odhadovaný čas: 18 minut na procesoru Heron r3 (POZNÁMKA: Jde pouze o odhad. Tvůj čas se může lišit.)

Výsledky učení​

  • Jak se spektrum neelastického rozptylu neutronů mapuje na dynamický strukturní faktor S(q,ω)S(q, \omega) 1D kvantového magnetu.

  • Jak připravit základní stav KCuF3_3 (izotropní Heisenbergův model) pomocí renormalizační grupy matice hustoty (DMRG) a maximalizace fidelity matrix product state (MPS).

  • Jak spustit časový vývoj Trotter, kompresi circuitu pomocí aproximované kvantové kompilace (AQC) a zmírněné (mitigated) spuštění v jediném volání funkce.

  • Jak post-processovat časovou řadu ⟨σz⟩(t)\langle \sigma_z \rangle(t) pro jednotlivá místa do S(q,ω)S(q, \omega) a identifikovat two-spinon kontinuum.

Předpoklady​

Souvislosti​

Neelastický rozptyl neutronů měří dynamický strukturní faktor S(q,ω)S(q, \omega), časoprostorovou Fourierovu transformaci korelační funkce spin-spin, takže reprodukce S(q,ω)S(q, \omega) z mikroskopického spinového modelu je přímý, falzifikovatelný test kvantové simulace. Tento tutoriál studuje KCuF3_3, antiferomagnetický Heisenbergův řetězec se spinem 12\frac{1}{2}, jehož excitace nejsou jednoduché překlopení spinu, ale páry frakcionalizovaných spinonů: místo ostré magnonové disperze vykazuje S(q,ω)S(q, \omega) široké two-spinon kontinuum, ohraničené zdola π2∣sin⁡q∣\tfrac{\pi}{2}|\sin q| a shora π∣sin⁡(q/2)∣\pi|\sin(q/2)|. To jsou přerušované křivky na následujících grafech. Celá fyzika a srovnání s naměřenými neutronovými daty jsou popsány v původním tutoriálu a v Lee et al., arXiv:2603.15608.

Kvantový pracovní postup napodobuje rozptylový experiment:

  1. Připrav základní stav řetězce ∣ψ0⟩|\psi_0\rangle.

  2. Kopni do něj lokální poruchou na centrálním místě, rotací ZZ o π/2\pi/2, napodobující přenos hybnosti a energie od neutronu.

  3. Časově vyviň pod Heisenbergovým hamiltoniánem, e−iHte^{-iHt}, pomocí Trotter produktového vzorce.

  4. Změř magnetizaci na jednotlivých místech ⟨σzj⟩(t)\langle \sigma_z^j \rangle(t). Jako funkce místa jj a času tt jde přesně o zpožděnou Greenovu funkci GR(j,jc,t)G^R(j, j_c, t), takže před Fourierovou transformací v kroku 5 není potřeba žádná konverze.

  5. Fourierovsky transformuj GRG^R na S(q,ω)S(q, \omega).

Problémy mohou nastat v kroku 3, kdy se přesné Trotter circuits pro dlouhé evoluce stanou pro hardware příliš hluboké. AQC s tenzorovými sítěmi tento problém řeší kompresí bloku Trotter kroků do pevného, mělkého parametrizovaného ansatzu, jehož fidelita stavu vůči přesné evoluci je klasicky maximalizována pomocí MPS simulátoru (arXiv:2301.08609). Šablona AQC Dynamics zabaluje celé toto kvantové jádro (Trotter syntézu, AQC kompresi a zmírněné spuštění) za jediné volání:

PRE (tento notebook)FUNKCE (aqc-dynamics-function)POST (tento notebook)
Základní stav z DMRG plus maximalizace MPS fidelity, s neutronovým kopnutím zapečeným ve stejném circuituTrotter syntéza → AQC komprese → spuštění na statevector, fake, nebo runtime, vracející ⟨σzj⟩(t)\langle \sigma_z^j \rangle(t) pro každé místoS(q,ω)S(q, \omega), dynamický strukturní faktor

Práce specifická pro daný experiment zůstává zde v notebooku: příprava základního stavu (PRE) a post-processing S(q,ω)S(q, \omega) (POST). Dva na kvantum náročné kroky, komprese a spuštění, běží uvnitř funkce.

Tento tutoriál je doplňkem k Simulaci neutronového rozptylu v kvantových materiálech pomocí kvantových circuits, který sestavuje stejný experiment inline: stejný model KCuF3_3, přípravu základního stavu, neutronové kopnutí a post-processing, s Trotter syntézou, AQC kompresí a zmírněným spuštěním vypsanými krok za krokem. Přečti si onen tutoriál, abys pochopil, jak funguje AQC komprese. Přečti si tento, abys spustil stejný experiment přes nasazenou funkční šablonu: kvantové jádro se stane jediným voláním funkce a vícehodinová AQC komprese běží uvnitř Serverless workeru místo na tvém počítači, takže při jejím běhu nepotřebuješ HPC systém ani otevřené jádro (kernel). Stejné volání pohání i další experimenty s 1D dynamikou.

Požadavky​

Než začneš s tímto tutoriálem, ujisti se, že máš následující:

  • Funkci nasazenou na svém účtu Qiskit Serverless. Nejprve spusť doprovodnou funkční šablonu: Nasazení a spuštění funkční šablony AQC + Trotter dynamics. Tato příručka tě provede získáním zdrojových souborů a nahráním funkce do tvého účtu. Tento tutoriál pouze volá nasazenou funkci.

  • Přihlašovací údaje IBM Quantum® uložené pro QiskitServerless (viz funkční šablona). Oba příklady v tomto tutoriálu volají nasazenou funkci, takže je oba potřebují.

  • Qiskit SDK v2.0 nebo novější (pip install qiskit).

  • Klient Qiskit IBM Catalog (pip install qiskit-ibm-catalog).

  • NumPy, SciPy a Matplotlib (pip install numpy scipy matplotlib). SciPy 1.14 nebo novější je potřeba pro optimalizátor COBYQA použitý při přípravě základního stavu.

  • Zásobník tenzorových sítí AQC, protože příprava základního stavu v kroku 1 běží lokálně v tomto notebooku: pip install 'qiskit-addon-aqc-tensor[quimb-jax]==0.3.1'.

První volání nově nasazené funkce čeká, dokud si Serverless worker nenainstaluje své závislosti, takže při tomto běhu očekávej vyšší latenci.

Nastavení​

Importuj knihovny a definuj pomocné funkce specifické pro experiment, které se použijí později: build_gs_ansatz (Hamiltonian variational ansatz, neboli HVA, pro přípravu základního stavu), prepare_ground_state (DMRG plus maximalizace MPS fidelity) a get_spectrum, plot_green a plot_spectrum (post-processing S(q,ω)S(q, \omega)). Tyto jsou upraveny z původního tutoriálu o neutronovém rozptylu.

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-addon-aqc-tensor qiskit-ibm-catalog quimb scipy
from functools import partial

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize

import quimb.tensor as qtn
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit_addon_aqc_tensor.simulation import tensornetwork_from_circuit
from qiskit_addon_aqc_tensor.simulation.quimb import QuimbSimulator
from qiskit_ibm_catalog import QiskitServerless
#  Dynamical structure factor via discrete Fourier transform

def get_spectrum(n, Gjjc, dt, time_steps, q_steps, w_steps):
    """Compute the dynamical structure factor from the retarded Green's function.

    Uses the center-site approximation and a discrete Fourier transform.
    """
    green = Gjjc / 4  # sigma -> S=1/2
    omega_max = np.pi / dt
    qpoints = np.arange(0, 2 * np.pi, 2 * np.pi / q_steps)
    omegas = np.arange(0, omega_max, omega_max / w_steps)
    green_map = np.zeros((omegas.shape[0], qpoints.shape[0]))
    center = n // 2 - 1
    for iw, w in enumerate(omegas):
        exponent = np.exp(1j * w * dt * np.arange(1, time_steps + 1))
        S_w = np.dot(green.T, exponent) * dt
        for iq, q in enumerate(qpoints):
            q_matrix = np.exp(-1j * q * np.arange(-center, center + 2, 1))
            green_map[iw, iq] = np.imag(np.dot(S_w, q_matrix))
    return green_map

#  Plotting helpers

def plot_spectrum(
    dsf,
    dt,
    q_steps,
    w_steps,
    lower_bound=False,
    upper_bound=False,
    title=None,
):
    """Heat-map of the dynamical structure factor."""
    omega_max = np.pi / dt
    qpoints = np.arange(0, 2 * np.pi, 2 * np.pi / q_steps)
    omegas = np.arange(0, omega_max, omega_max / w_steps)
    x, y = np.meshgrid(qpoints, omegas)
    fig, ax = plt.subplots(figsize=(8, 5))
    c = ax.pcolormesh(x, y, dsf / np.max(dsf), cmap="viridis", shading="auto")
    fig.colorbar(c, ax=ax, label="Normalized intensity")
    if lower_bound:
        ax.plot(
            qpoints,
            np.pi * np.abs(np.sin(qpoints)) / 2,
            "--",
            color="white",
            lw=1.5,
            label="Lower bound",
        )
    if upper_bound:
        ax.plot(
            qpoints,
            np.pi * np.abs(np.sin(qpoints / 2)),
            "--",
            color="red",
            lw=1.5,
            label="Upper bound",
        )
    ax.set_ylim(0, 3.6)
    ax.set_xlim(0, 2 * np.pi - 2 * np.pi / q_steps)
    ax.set_xlabel(r"$q$", fontsize=16)
    ax.set_ylabel(r"$\tilde{\omega} = \omega / J$", fontsize=16)
    ax.set_xticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi])
    ax.set_xticklabels(["0", r"$\pi/2$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"])
    if lower_bound or upper_bound:
        ax.legend(loc="upper right", fontsize=11)
    if title:
        ax.set_title(title, fontsize=14)
    plt.tight_layout()
    plt.show()

def plot_green(n, Gjjc, time_steps, dt, title=None):
    """Heat-map of the retarded Green's function in real space and time."""
    fig, ax = plt.subplots(figsize=(8, 6))
    t_axis = np.arange(1, time_steps + 1) * dt
    site_axis = np.arange(n)
    x, y = np.meshgrid(t_axis, site_axis)
    c = ax.pcolormesh(
        x,
        y,
        np.real(Gjjc).T,
        cmap="RdBu",
        vmax=0.5,
        vmin=-0.5,
        shading="auto",
    )
    fig.colorbar(c, ax=ax, label=r"Re $G^R(j, j_c, t)$")
    ax.set_xlabel(r"Time  ($t / J^{-1}$)", fontsize=16)
    ax.set_ylabel("Site index $j$", fontsize=16)
    if title:
        ax.set_title(title, fontsize=14)
    plt.tight_layout()
    plt.show()

#  Variational ground-state ansatz (HVA)

def _apply_xxz_pair_gate(qc, q0, q1, theta):
    """Apply the parameterized XXZ-type two-qubit gate used in the HVA."""
    qc.cx(q0, q1)
    qc.rz(theta, q1)
    qc.h(q0)
    qc.rz(theta + np.pi / 2, q0)
    qc.cx(q0, q1)
    qc.rz(-theta, q1)
    qc.h(q1)
    qc.cx(q1, q0)
    qc.rz(np.pi / 2, q1)
    qc.rz(-np.pi / 2, q0)
    qc.h(q1)
    qc.h(q0)

def build_gs_ansatz(n, params, layers):
    """Build the Hamiltonian variational ansatz (HVA) circuit for
    ground-state preparation of the 1D Heisenberg model.

    Starts from a product of singlet pairs and applies alternating
    odd/even layers of parameterized XXZ gates. For layer r,
    params[2 * r] is the odd-layer (inter-pair) angle and
    params[2 * r + 1] is the even-layer (intra-pair) angle.
    """
    qc = QuantumCircuit(n)
    # Initial singlet product state
    for i in range(n // 2):
        qc.x(2 * i)
        qc.x(2 * i + 1)
        qc.h(2 * i + 1)
        qc.cx(2 * i + 1, 2 * i)
    # Variational layers
    for r in range(layers):
        for i in range(1, (n + 1) // 2):  # odd layer
            _apply_xxz_pair_gate(qc, 2 * i - 1, 2 * i, params[2 * r])
        for i in range(n // 2):  # even layer
            _apply_xxz_pair_gate(qc, 2 * i, 2 * i + 1, params[2 * r + 1])
    return qc

def prepare_ground_state(n, gs_layers=5, max_bond=128, cutoff=1e-8):
    """Prepare the KCuF3 (isotropic Heisenberg) ground state as a QuantumCircuit.

    Runs DMRG (quimb MPO + DMRG2) to get the chain's ground state, then optimizes
    the HVA angles to maximize the MPS overlap |<psi_ansatz|psi_DMRG>|^2. No exact
    diagonalization, so it scales to larger n.
    """
    J = Jz = 1.0
    builder = qtn.SpinHam1D(S=1 / 2)
    builder += J * 0.5, "+", "-"
    builder += J * 0.5, "-", "+"
    builder += Jz, "Z", "Z"
    H_mpo = builder.build_mpo(L=n)
    dmrg = qtn.DMRG2(H_mpo)
    dmrg.solve(tol=1e-8, verbosity=0)

    gs_sim = QuimbSimulator(
        quimb_circuit_factory=partial(
            qtn.CircuitMPS, gate_opts=dict(cutoff=cutoff, max_bond=max_bond)
        ),
        autodiff_backend="jax",
    )

    def gs_infidelity(params):
        psi = tensornetwork_from_circuit(
            build_gs_ansatz(n, params, gs_layers), gs_sim
        ).psi
        return 1 - abs(psi.H @ dmrg.state) ** 2

    # Seed and optimizer match the original tutorial. Each layer starts at
    # [0, pi/2]: an odd-layer angle of 0 makes the inter-pair gate the identity,
    # and an even-layer angle of pi/2 makes the intra-pair gate a SWAP (since
    # 0.5 * (XX + YY + ZZ) = SWAP - I/2). That puts the seed at the singlet-pair
    # product limit, which is already a decent approximation to the Heisenberg
    # ground state, so the optimizer only has to refine it. The small jitter
    # (fixed RNG seed, so runs are reproducible) breaks the exact symmetry
    # between layers; COBYQA then runs for up to 100 iterations.
    rng = np.random.default_rng(12345)
    x0 = np.tile([0.0, np.pi / 2], gs_layers) + rng.normal(
        scale=0.1, size=2 * gs_layers
    )
    result_gs = scipy.optimize.minimize(
        gs_infidelity, x0, method="COBYQA", options={"maxiter": 100}
    )
    print(f"DMRG ground-state energy: {dmrg.energy:.6f}")
    print(f"GS fidelity: {1 - result_gs.fun:.4f}")
    return build_gs_ansatz(n, result_gs.x, gs_layers)

print("Setup complete - helpers defined.")
Setup complete - helpers defined.

Načtení funkční šablony​

Připoj se ke Qiskit Serverless a načti nasazenou funkci aqc-dynamics-function. Oba příklady v tomto tutoriálu volají stejný handle fn, takže funkce se zde načte jen jednou.

# Credentials are read from the account saved once via QiskitServerless.save_account(...)
serverless = QiskitServerless()
fn = serverless.load("aqc-dynamics-function")

Příklad s malým simulátorem​

Nejprve spustíme celý workflow na malém 10místném řetězci pomocí přesného backendu statevector. Tím ověříme pipeline PRE → FUNCTION → POST ještě předtím, než strávíme jakýkoli čas na QPU.

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

Sestav Hamiltonián KCuF3_3 jako SparsePauliOp (izotropní Heisenberg: XX+YY+ZZXX + YY + ZZ s vazbou 14\tfrac14 na každé vazbě nejbližších sousedů; řetězce jsou Pauliho operátory, takže 14\tfrac14 dává vazbu pro spin-12\frac{1}{2}). Připrav základní stav pomocí DMRG a maximalizace MPS fidelity, poté zapracuj neutronový kop: π/2\pi/2 ZZ-rotaci na centrálním místě. Připravený obvod je to, co předáváme funkci jako initial_state. Parametr observables ponecháváme na výchozí hodnotě (per-site ZZ), což je přesně čtení ⟨σzj⟩(t)\langle \sigma_z^j \rangle(t), které neutronový workflow potřebuje.

n = 10
dt = 0.6  # physical time per Trotter step (also the omega-axis unit in POST)
time_steps = 10
center = n // 2 - 1

# MPS-simulator settings, shared by the ground-state prep here and the AQC
# compression inside the function (matches the original tutorial).
mps_max_bond = 32
mps_cutoff = 1e-8

# 1D isotropic Heisenberg (KCuF3) Hamiltonian on n qubits
H = SparsePauliOp.from_sparse_list(
    [(p, [i, i + 1], 0.25) for i in range(n - 1) for p in ("XX", "YY", "ZZ")],
    num_qubits=n,
)

# Ground state (DMRG + fidelity max) + neutron kick baked into the same circuit
gs_circuit = prepare_ground_state(
    n, gs_layers=3, max_bond=mps_max_bond, cutoff=mps_cutoff
)
gs_circuit.rz(
    np.pi / 2, center
)  # exp(-i (pi/2)/2 Z_center): the neutron perturbation
print(
    f"Prepared {n}-qubit ground state with the neutron kick at site {center}."
)
DMRG ground-state energy: -4.258035
GS fidelity: 0.9841
Prepared 10-qubit ground state with the neutron kick at site 4.

Kroky 2 a 3: Zkomprimuj a spusť pomocí function template​

V ručně psaném workflow jsou to dvě samostatné fáze: optimalizace obvodů pro hardware (krok 2) a jejich spuštění (krok 3). Function template obě fáze slučuje do jednoho volání. Provádí Trotterovu syntézu, AQC kompresi a transpilaci pro hardware, a poté obvody spustí (zde na přesném simulátoru, později na hardwaru se zabudovanou mitigací chyb). Dva ladicí parametry jsou aqc_segments (kompresní plán) a aqc_options (nastavení MPS a optimalizátoru). Každý segment {"n_steps": k, "ansatz_steps": m} zkomprimuje k po sobě jdoucích Trotterových kroků do ansatzu postaveného z m-krokového Trotterova cíle, a jakékoli kroky nad rámec sum(n_steps) běží jako čistý Trotter. Rané kroky s nízkou provázaností se dobře komprimují do mělkého ansatzu (ansatz_steps=1), takže zde komprimujeme první tři kroky do jednovrstvého ansatzu a další dva do hlubšího dvouvrstvého ansatzu; zbývajících pět z 10 Trotterových kroků běží jako čistý Trotter. Pro aqc_options zrcadlíme původní tutoriál: MPS bond dimension max_bond=32, cutoff=1e-8 a optimalizátor L-BFGS-B omezený na 100 iterací.

Zavolej funkci načtenou v nastavení. backend="statevector" spouští přesnou referenční cestu: žádný čas na QPU, obvody běží na přesném statevector simulátoru uvnitř serverless workeru (uložený účet Qiskit Serverless je pro její volání stále potřeba). Parametr initial_state nese připravený základní stav (včetně kopu); observables je vynechán, takže funkce měří výchozí per-site ZZ.

job = fn.run(
    t_steps=time_steps,
    aqc_segments=[
        {
            "n_steps": 3,
            "ansatz_steps": 1,
        },  # early steps -> shallow 1-layer ansatz
        {
            "n_steps": 2,
            "ansatz_steps": 2,
        },  # later steps -> deeper 2-layer ansatz
    ],
    aqc_options={
        "max_bond": mps_max_bond,  # MPS bond dimension for AQC compression
        "cutoff": mps_cutoff,
        "optimizer_settings": {
            "method": "L-BFGS-B",
            "jac": True,
            "options": {"maxiter": 100},
        },
    },
    dt=dt,
    hamiltonian=H,
    initial_state=gs_circuit,  # prepared ground state including the neutron kick
    # observables omitted -> default per-site Z (the neutron sigma_z readout)
    backend="statevector",
)
print(job.status())  # rerun this cell until status says DONE
DONE
# The per-site <sigma_z>(t) the function returns is the retarded Green's function
# G(j, j_c, t). The workflow samples t = 1..time_steps, so drop the t = 0 row (the
# prepared+kicked state before any evolution) before post-processing.
result = job.result()
print(
    "AQC fidelities:",
    {k: round(v, 4) for k, v in result["metadata"]["aqc_fidelities"].items()},
)

ev = np.array(result["expectation_values"])
Gjjc = ev[1:]  # shape (time_steps, n)
print("Green's function shape:", Gjjc.shape)
AQC fidelities: {'1': 1.0, '2': 0.9999, '3': 0.9992, '4': 0.9998, '5': 0.9995}
Green's function shape: (10, 10)

Krok 4: Zpracuj výsledek a vrať ho v požadovaném klasickém formátu​

Fourierovsky transformuj Greenovu funkci na S(q,ω)S(q, \omega), zrcadlově symetrizuj a ořízni záporné hodnoty: standardní neutronové postprocessing. Zrcadlení je přesné, protože pro tento model platí S(q,ω)=S(−q,ω)S(q, \omega) = S(-q, \omega), a přežívající záporné hodnoty jsou artefakty Fourierovy transformace konečné, diskrétně vzorkované časové řady, takže se ořezávají na nulu. Na tomto malém přesném běhu je dvouspinonové kontinuum rozlišeno jen hrubě, ale mechanismus je totožný s následujícím hardwarovým během.

q_res, w_res = 100, 100
spectrum = get_spectrum(n, Gjjc, dt, time_steps, q_res, w_res)
spectrum = -(spectrum + spectrum[:, ::-1]) / 2  # mirror symmetry
spectrum = np.clip(spectrum, a_min=0, a_max=None)  # clip negatives

plot_green(
    n,
    Gjjc,
    time_steps,
    dt,
    title=f"Retarded Green's function - {n} qubits (AQC, statevector)",
)
plot_spectrum(
    spectrum,
    dt,
    q_res,
    w_res,
    lower_bound=True,
    upper_bound=True,
    title=f"Dynamical structure factor - {n} qubits (AQC, statevector)",
)

Output of the previous code cell

Output of the previous code cell

Rozsáhlý hardwarový příklad​

Stejný workflow se škáluje nahoru bez jakékoli změny vědeckého kódu: 30místný řetězec, dvojnásobná hloubka Trotterova rozvoje (20 kroků), kompresní plán, který mění hloubku ansatzu (hlubší ansatz pro pozdější, více provázané kroky), a spuštění na procesoru IBM Quantum se zabudovanou mitigací chyb funkce (dynamické odsdružení, Pauliho twirling a twirled readout error extinction (TREX)). Projdeme si stejné čtyři kroky jako v příkladu se simulátorem a znovu použijeme handle fn z nastavení.

Malé měřítkoVelké měřítko
Qubity1030
Trotterovy kroky1020
AQC-komprimované kroky (1vrstvý + 2vrstvý)3 + 2 = 56 + 4 = 10
Vrstvy ansatzu základního stavu35
Maximální bond dimension MPS32128
BackendstatevectorQPU s DD, Pauliho twirlingem a TREX

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

Sestav stejný Heisenbergův SparsePauliOp pro KCuF3_3 a připrav základní stav, nyní s hlubším ansatzem gs_layers=5 pro delší řetězec, a poté zapracuj π/2\pi/2 ZZ neutronový kop na centrálním místě. Toto je totožné s mapováním pro malé měřítko, ale při n=30n = 30.

Očekávej nižší fidelitu základního stavu než u 10místného běhu: zhruba 0,82 zde oproti 0,98 pro menší řetězec, protože pět vrstev HVA nedokáže plně zachytit základní stav 30místného řetězce. To je očekávané, nikoli selhání, a původní tutoriál akceptuje ze stejného důvodu zhruba 0,65 u 50 míst. Zvýšení gs_layers nebo limitu iterací COBYQA to zlepší za cenu dalších klasických výpočtů.

n = 30
dt = 0.6
time_steps = 20
center = n // 2 - 1

# Same MPS settings as the original large-scale run: a larger bond for the
# longer, more-entangled chain (shared by GS prep and AQC compression).
mps_max_bond = 128
mps_cutoff = 1e-8

# Same KCuF3 Hamiltonian and ground-state prep, on a larger chain
H = SparsePauliOp.from_sparse_list(
    [(p, [i, i + 1], 0.25) for i in range(n - 1) for p in ("XX", "YY", "ZZ")],
    num_qubits=n,
)
gs_circuit = prepare_ground_state(
    n, gs_layers=5, max_bond=mps_max_bond, cutoff=mps_cutoff
)
gs_circuit.rz(np.pi / 2, center)  # neutron kick at the center site
print(
    f"Prepared {n}-qubit ground state with the neutron kick at site {center}."
)
DMRG ground-state energy: -13.111355
GS fidelity: 0.8201
Prepared 30-qubit ground state with the neutron kick at site 14.

Kroky 2 a 3: Zkomprimuj a spusť pomocí function template​

Stejné jednotlivé volání jako v příkladu se simulátorem, nyní s backend_name odkazujícím na procesor IBM Quantum, takže funkce tam obvody transpiluje a spouští. Kompresní plán mění hloubku ansatzu: prvních šest (nízce provázaných) Trotterových kroků se komprimuje do mělkého jednovrstvého ansatzu, dalších čtyři do hlubšího dvouvrstvého ansatzu a zbývajících 10 z 20 kroků běží jako čistý Trotter. aqc_options zvyšuje MPS bond dimension na max_bond=128 pro delší, více provázaný řetězec (v souladu s originálem), přičemž zachovává stejný optimalizátor L-BFGS-B omezený na 100 iterací. estimator_options zapíná zabudovanou mitigaci chyb: dynamické odsdružení (XY4), gate twirling a mitigaci měření TREX. Výchozí hodnoty funkce se již shodují s původním tutoriálem pro všechny tyto vlastnosti kromě rozpočtu na učení TREX (measure_noise_learning). Celý blok je přesto vypsán, protože volajícím dodané estimator_options nahrazují výchozí hodnoty funkce úplně, místo aby se s nimi slučovaly, takže vynechání klíče by vedlo k výchozí hodnotě IBM Quantum Compute místo výchozí hodnoty funkce.

# Steps 2 + 3: the function compresses (varied ansatz) and executes on hardware.
job = fn.run(
    t_steps=time_steps,
    aqc_segments=[
        {
            "n_steps": 6,
            "ansatz_steps": 1,
        },  # early steps -> shallow 1-layer ansatz
        {
            "n_steps": 4,
            "ansatz_steps": 2,
        },  # later steps -> deeper 2-layer ansatz
    ],
    aqc_options={
        "max_bond": mps_max_bond,  # 128 for the longer chain
        "cutoff": mps_cutoff,
        "optimizer_settings": {
            "method": "L-BFGS-B",
            "jac": True,
            "options": {"maxiter": 100},
        },
    },
    dt=dt,
    hamiltonian=H,
    initial_state=gs_circuit,
    backend_name="ibm_pittsburgh",
    # Mitigation settings from the original tutorial. Only the two
    # measure_noise_learning values differ from the function's defaults; the rest
    # restates them, because a caller-supplied estimator_options dict replaces the
    # function's defaults wholesale rather than merging into them.
    estimator_options={
        "environment": {"job_tags": ["TUT-SNS"]},
        "dynamical_decoupling": {"enable": True, "sequence_type": "XY4"},
        "twirling": {
            "enable_gates": True,
            "num_randomizations": 1000,
            "shots_per_randomization": 128,
        },
        "resilience": {
            "measure_mitigation": True,
            "measure_noise_learning": {
                "num_randomizations": 32,
                "shots_per_randomization": 100,
            },
        },
    },
)
print("job ID (save this to reconnect later):", job.job_id)
job ID (save this to reconnect later): 43ed8d07-6d7d-4f33-b70a-7f31b765b310
Znovupřipojení k dlouho běžící úloze

Rozsáhlý běh není rychlý a většina času se spotřebuje klasicky, nikoli na QPU. AQC komprese běží uvnitř funkce ještě předtím, než cokoli dorazí na QPU: při 30 místech s max_bond=128 to v našem běhu trvalo téměř čtyři hodiny, oproti zhruba 18 minutám času na QPU uvedeným v odhadu využití na začátku tohoto tutoriálu. K oběma se navíc přičítá čekání ve frontě. Tento notebook ani kernel nemusíš mít po dobu běhu otevřený.

Zkopíruj ID úlohy vypsané předchozí buňkou a ulož si ho. Následující tři buňky ti umožní se k úloze později vrátit:

  1. Znovupřipojení, potřebné pouze v nové kernel session: znovu spusť buňky Nastavení, aby se znovu vytvořil serverless, a poté znovu sestav handle job z uloženého ID. Tuto buňku přeskoč, pokud jsi stále v session, ve které jsi úlohu odeslal, protože handle je již aktivní.

  2. Zkontroluj stav: spouštěj znovu, dokud nenahlásí DONE.

  3. Získej výsledek: spusť pouze poté, co je stav DONE.

Následující buňka pro znovupřipojení obsahuje zástupný symbol. Nahraď ho vlastním job_id:

# Reconnect to a previously submitted job by its ID. Only needed in a NEW kernel
# session; if you are still in the session where you submitted, the `job` handle
# from the preceding cell is already live, so skip this cell. Replace the ID that follows with your own.
job = serverless.get_job_by_id("<your job ID>")
# Check where the job is. Re-run this until it reports DONE before fetching the
# result in the following cell: QUEUED -> INITIALIZING -> RUNNING: OPTIMIZING_FOR_HARDWARE ->
# RUNNING: WAITING_FOR_QPU -> RUNNING: EXECUTING_QPU -> RUNNING: POST_PROCESSING
# -> DONE.
print(job.status())
DONE
# Run this only once the preceding status cell reports DONE. result() blocks until
# the job finishes, so calling it earlier just waits (possibly for hours).
result = job.result()
print(
    "AQC fidelities:",
    {k: round(v, 4) for k, v in result["metadata"]["aqc_fidelities"].items()},
)

ev = np.array(result["expectation_values"])
Gjjc = ev[1:]  # drop the t = 0 row -> shape (time_steps, n)
AQC fidelities: {'1': 1.0, '2': 0.9994, '3': 0.9944, '4': 0.9853, '5': 0.9747, '6': 0.959, '7': 0.9495, '8': 0.9542, '9': 0.9533, '10': 0.9451}

Krok 4: Zpracuj výsledek a vrať ho v požadovaném klasickém formátu​

Postprocessing totožný s během na simulátoru: Fourierovsky transformuj Greenovu funkci na S(q,ω)S(q, \omega), zrcadlově symetrizuj a ořízni záporné hodnoty. S delším řetězcem a delší evolucí je dvouspinonové kontinuum mnohem lépe rozlišeno. Mělo by vyplnit pásmo mezi přerušovanými mezemi, nejjasněji poblíž q=πq = \pi.

n = result["metadata"]["n"]
q_res, w_res = 100, 100
spectrum = get_spectrum(n, Gjjc, dt, time_steps, q_res, w_res)
spectrum = -(spectrum + spectrum[:, ::-1]) / 2  # mirror symmetry
spectrum = np.clip(spectrum, a_min=0, a_max=None)  # clip negatives

plot_green(
    n,
    Gjjc,
    time_steps,
    dt,
    title=f"Retarded Green's function - {n} qubits (AQC, hardware)",
)
plot_spectrum(
    spectrum,
    dt,
    q_res,
    w_res,
    lower_bound=True,
    upper_bound=True,
    title=f"Dynamical structure factor - {n} qubits (AQC, hardware)",
)

Output of the previous code cell

Output of the previous code cell

Dodatek​

Předchozí hardwarový příklad spouští jednu konkrétní délku řetězce. Následující tři spektra pocházejí z dřívějších hardwarových běhů téhož workflow na ibm_pittsburgh při 10, 20 a 30 místech, se všemi ostatními vstupy pevně danými: 20 Trotterových kroků při dt = 0.6, kompresní plán šesti jednovrstvých plus čtyř dvouvrstvých AQC-komprimovaných kroků a max_bond = 128. Jde o zaznamenané výsledky, nikoli o výstup z předchozích buněk.

Ve všech třech velikostech se používá stejné nastavení, takže jsou spektra přímo porovnatelná. Ladění pro každou délku řetězce zvlášť, například s více vrstvami ansatzu základního stavu nebo větším max_bond, může přinést lepší výsledky, než jaké jsou ukázány zde.

Dynamical structure factor at 10 sites, a single sharp bright peak at q = pi near the lower bound

Dynamical structure factor at 20 sites, spectral weight filling the band between the two dashed two-spinon bounds

Dynamical structure factor at 30 sites, the continuum resolved more finely with fainter contrast and some weight outside the bounds

Další kroky​

Doporučení
Zdroj: IBM Quantum docs — aktualizováno 9. 9. 2026
Anglická verze na doQumentation — aktualizováno 9. 9. 2026
Problém s webem nebo překladem?Problém v obsahu? Upravit v IBM Quantum docs