דלג לתוכן הראשי

סימולציה של פיזור נויטרונים עם תהליך עבודה Serverless של AQC + דינמיקת Trotter

הערכת שימוש: 18 דקות על מעבד Heron r3 (הערה: זוהי הערכה בלבד. זמן הריצה שלך עשוי להשתנות.)

תוצאות למידה​

  • איך ספקטרום פיזור נויטרונים לא אלסטי ממופה לגורם המבנה הדינמי S(q,ω)S(q, \omega) של מגנט קוונטי חד-ממדי.

  • איך להכין את מצב היסוד של KCuF3_3 (Heisenberg איזוטרופי) באמצעות חבורת הרנורמליזציה של מטריצת הצפיפות (DMRG) ומקסום נאמנות של מצב מכפלת מטריצות (MPS).

  • איך להריץ אבולוציית זמן של Trotter, דחיסת Circuit של קימפול קוונטי מקורב (AQC), וביצוע עם הפחתת שגיאות כקריאת פונקציה בודדת.

  • איך לעבד לאחר את סדרת הזמן ⟨σz⟩(t)\langle \sigma_z \rangle(t) לכל אתר ל-S(q,ω)S(q, \omega) ולזהות את הרצף הדו-ספינוני.

דרישות מוקדמות​

רקע​

פיזור נויטרונים לא אלסטי מודד את גורם המבנה הדינמי S(q,ω)S(q, \omega), טרנספורם פורייה במרחב ובזמן של פונקציית הקורלציה ספין-ספין, כך ששחזור S(q,ω)S(q, \omega) ממודל ספין מיקרוסקופי הוא מבחן ישיר וניתן להפרכה של סימולציה קוונטית. מדריך זה חוקר את KCuF3_3, שרשרת Heisenberg אנטיפרומגנטית עם ספין-12\frac{1}{2} שהעירורים שלה אינם היפוכי ספין בודדים אלא זוגות של ספינונים מפוצלים: במקום דיספרסיית מגנון חדה, S(q,ω)S(q, \omega) מציג רצף דו-ספינוני רחב, תחום מלמטה על ידי π2∣sin⁡q∣\tfrac{\pi}{2}|\sin q| ומלמעלה על ידי π∣sin⁡(q/2)∣\pi|\sin(q/2)|. אלו הן העקומות המקווקוות בתרשימים שיבואו. הפיזיקה במלואה, וההשוואה מול נתוני נויטרונים נמדדים, מכוסות במדריך המקורי וב-Lee et al., arXiv:2603.15608.

תהליך העבודה הקוונטי משקף את ניסוי הפיזור:

  1. הכן את מצב היסוד ∣ψ0⟩|\psi_0\rangle של השרשרת.

  2. תן בעיטה עם הפרעה מקומית באתר המרכזי, סיבוב π/2\pi/2 ZZ, שמחקה את העברת התנע והאנרגיה מהנייטרון.

  3. בצע אבולוציית זמן תחת האמילטוניאן הייזנברג, e−iHte^{-iHt}, עם נוסחת מכפלת טרוטר.

  4. מדוד את המגנטיזציה לכל אתר ⟨σzj⟩(t)\langle \sigma_z^j \rangle(t). כפונקציה של אתר jj וזמן tt, זו בדיוק פונקציית גרין המושהית GR(j,jc,t)G^R(j, j_c, t), כך שלא נדרשת המרה לפני התמרת פורייה בשלב 5.

  5. בצע התמרת פורייה על GRG^R לתוך S(q,ω)S(q, \omega).

בעיות עלולות להתעורר בשלב 3, כאשר מעגלי טרוטר מדויקים לאבולוציות ארוכות נעשים עמוקים מדי עבור החומרה. AQC עם רשתות טנזוריות מטפל בכך על ידי דחיסת בלוק של צעדי טרוטר לתוך אנסאץ פרמטרי קבוע ורדוד שהנאמנות של מצבו לאבולוציה המדויקת ממוקסמת קלאסית עם סימולטור MPS (arXiv:2301.08609). התבנית AQC Dynamics אורזת את כל הליבה הקוונטית הזו (סינתזת טרוטר, דחיסת AQC, וביצוע עם הפחתת שגיאות) מאחורי קריאה אחת:

PRE (המחברת הזו)FUNCTION (aqc-dynamics-function)POST (המחברת הזו)
מצב יסוד מ-DMRG בתוספת מיקסום נאמנות MPS, עם בעיטת הנייטרון אפויה לתוך אותו המעגלסינתזת טרוטר → דחיסת AQC → ביצוע על statevector, fake, או runtime, שמחזיר ⟨σzj⟩(t)\langle \sigma_z^j \rangle(t) לכל אתרS(q,ω)S(q, \omega), גורם המבנה הדינמי

העבודה הספציפית לניסוי נשארת כאן במחברת: הכנת מצב יסוד (PRE) ועיבוד ה-S(q,ω)S(q, \omega) שלאחר מכן (POST). שני השלבים הכבדים קוונטית, הדחיסה והביצוע, רצים בתוך הפונקציה.

מדריך זה הוא מלווה לסימולציית פיזור נייטרונים בחומרים קוונטיים עם מעגלים קוונטיים, שבונה את אותו הניסוי בתוך המחברת עצמה: אותו מודל KCuF3_3, הכנת מצב יסוד, בעיטת נייטרון, ועיבוד לאחר מכן, כאשר סינתזת טרוטר, דחיסת AQC, וביצוע עם הפחתת שגיאות כתובים צעד אחר צעד. קרא את המדריך ההוא כדי ללמוד איך דחיסת AQC עובדת. קרא את זה כדי להריץ את אותו הניסוי דרך תבנית פונקציה שנפרסה: הליבה הקוונטית הופכת לקריאת פונקציה בודדת, ודחיסת ה-AQC בת שעות רבות רצה בתוך עובד ה-Serverless במקום על המחשב שלך, כך שאינך זקוק למערכת HPC או ל-kernel פתוח בזמן שהיא רצה. אותה הקריאה מפעילה גם ניסויי דינמיקה חד-ממדית אחרים.

דרישות​

לפני שתתחיל במדריך זה, ודא שיש לך את הדברים הבאים:

  • הפונקציה פרוסה בחשבון Qiskit Serverless שלך. הרץ תחילה את תבנית הפונקציה המלווה: פרוס והרץ את תבנית פונקציית הדינמיקה AQC + Trotter. המדריך ההוא עובר על השגת קבצי המקור והעלאת הפונקציה לחשבונך. מדריך זה רק קורא לפונקציה הפרוסה.

  • אישורי IBM Quantum® שמורים עבור QiskitServerless (ראה את תבנית הפונקציה). שתי הדוגמאות במדריך זה קוראות לפונקציה הפרוסה, כך ששתיהן זקוקות להם.

  • Qiskit SDK גרסה 2.0 ומעלה (pip install qiskit).

  • לקוח Qiskit IBM Catalog (pip install qiskit-ibm-catalog).

  • NumPy, SciPy, ו-Matplotlib (pip install numpy scipy matplotlib). נדרש SciPy 1.14 ומעלה עבור אופטימייזר COBYQA המשמש בהכנת מצב היסוד.

  • מחסנית הרשתות הטנזוריות של AQC, מכיוון שהכנת מצב היסוד בשלב 1 רצה מקומית במחברת הזו: pip install 'qiskit-addon-aqc-tensor[quimb-jax]==0.3.1'.

הקריאה הראשונה לפונקציה שזה עתה נפרסה ממתינה בזמן שעובד ה-Serverless מתקין את התלויות שלו, כך שיש לצפות לזמן המתנה נוסף בהרצה הזו.

הגדרה​

ייבא את הספריות והגדר את העוזרים הספציפיים לניסוי שישמשו בהמשך: build_gs_ansatz (האנסאץ הוריאציוני של האמילטוניאן, או HVA, להכנת מצב יסוד), prepare_ground_state (DMRG בתוספת מיקסום נאמנות MPS), ו-get_spectrum, plot_green, ו-plot_spectrum (עיבוד ה-S(q,ω)S(q, \omega) שלאחר מכן). אלה מותאמים מהמדריך המקורי לפיזור נייטרונים.

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

טען את תבנית הפונקציה​

התחבר ל-Qiskit Serverless וטען את aqc-dynamics-function הפרוסה. שתי הדוגמאות במדריך זה קוראות לאותו ידית fn, כך שהפונקציה נטענת פעם אחת, כאן.

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

דוגמת סימולטור בקנה מידה קטן​

נתחיל בהרצת זרימת העבודה המלאה על שרשרת קטנה בת 10 אתרים באמצעות שימוש ב-Backend המדויק statevector. זה מאמת את צינור ה-PRE → FUNCTION → POST לפני שמבזבזים זמן QPU כלשהו.

שלב 1: מיפוי קלטים קלאסיים לבעיה קוונטית​

בנה את האמילטוניאן KCuF3_3 בתור SparsePauliOp (הייזנברג איזוטרופי: XX+YY+ZZXX + YY + ZZ בצימוד 14\tfrac14 בכל קשר שכנים קרובים; המחרוזות הן אופרטורי פאולי, כך ש-14\tfrac14 נותן את צימוד הספין-12\frac{1}{2}). הכן את מצב היסוד עם DMRG בתוספת מיקסום נאמנות MPS, ואז אפה את בעיטת הנייטרון: סיבוב π/2\pi/2 ZZ באתר המרכזי. המעגל שהוכן הוא מה שאנו מוסרים לפונקציה כ-initial_state. אנו משאירים את observables בברירת המחדל שלו (ZZ לכל אתר), שהיא בדיוק קריאת ה-⟨σzj⟩(t)\langle \sigma_z^j \rangle(t) שזרימת העבודה של הנייטרון זקוקה לה.

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.

שלבים 2 ו-3: דחוס ובצע עם תבנית הפונקציה​

בזרימת עבודה כתובה ידנית אלה שני שלבים נפרדים: אופטימיזציה של המעגלים לחומרה (שלב 2) וביצועם (שלב 3). תבנית הפונקציה מכווצת את שניהם לקריאה אחת. היא מבצעת סינתזת טרוטר, דחיסת AQC, ותרגום לחומרה, ואז מריצה את המעגלים (כאן על הסימולטור המדויק, מאוחר יותר עם הפחתת שגיאות מובנית על חומרה). שני פרמטרי הכוונון הם aqc_segments (תוכנית הדחיסה) ו-aqc_options (הגדרות ה-MPS והאופטימייזר). כל מקטע {"n_steps": k, "ansatz_steps": m} דוחס k צעדי טרוטר עוקבים לתוך אנסאץ הבנוי מיעד טרוטר בן m צעדים, וכל צעד שמעבר ל-sum(n_steps) רץ כטרוטר רגיל. צעדים מוקדמים ובעלי שזירות נמוכה נדחסים היטב לתוך אנסאץ רדוד (ansatz_steps=1), כך שכאן אנו דוחסים את שלושת הצעדים הראשונים לתוך אנסאץ בעל שכבה אחת ואת השניים הבאים לתוך אנסאץ עמוק יותר בעל שתי שכבות; חמשת הצעדים הנותרים מתוך 10 צעדי טרוטר רצים כטרוטר רגיל. עבור aqc_options אנו משקפים את המדריך המקורי: ממד קשר MPS max_bond=32, cutoff=1e-8, ואופטימייזר L-BFGS-B מוגבל ל-100 איטרציות.

קרא לפונקציה שנטענה בהגדרה. backend="statevector" מריץ את נתיב ההתייחסות המדויק: ללא זמן QPU, כאשר המעגלים רצים על סימולטור וקטור מצב מדויק בתוך עובד ה-serverless (עדיין נדרש חשבון Qiskit Serverless שמור כדי לקרוא לו). ה-initial_state נושא את מצב היסוד המוכן (כולל הבעיטה); observables מושמט כך שהפונקציה מודדת את ברירת המחדל 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)

שלב 4: עבד ולאחר מכן החזר תוצאה בפורמט קלאסי רצוי​

בצע התמרת פורייה של פונקציית גרין לתוך S(q,ω)S(q, \omega), בצע סימטריזציה במראה, וקטום ערכים שליליים: העיבוד הסטנדרטי של הנייטרון לאחר מכן. הסימטריזציה מדויקת מכיוון ש-S(q,ω)=S(−q,ω)S(q, \omega) = S(-q, \omega) עבור מודל זה, והערכים השליליים ששורדים הם ארטיפקטים של התמרת פורייה על סדרת זמן סופית ודגומה באופן בדיד, כך שהם נקטמים לאפס. בהרצה המדויקת הקטנה הזו, הרצף הדו-ספינוני מובחן רק באופן גס, אך המנגנון זהה להרצת החומרה שתבוא בהמשך.

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

דוגמת חומרה בקנה מידה גדול​

אותה זרימת עבודה עולה בקנה מידה מבלי לשנות שום קוד מדעי: שרשרת בת 30 אתרים, כפול עומק טרוטר (20 צעדים), תוכנית דחיסה שמשנה את עומק האנסאץ (אנסאץ עמוק יותר לצעדים המאוחרים, השזורים יותר), וביצוע על מעבד IBM Quantum עם הפחתת השגיאות המובנית של הפונקציה (ניתוק דינמי, טוויירלינג פאולי, והיכחדות שגיאת קריאה בטוויירלינג (TREX)). נעבור על אותם ארבעה שלבים כמו בדוגמת הסימולטור, תוך שימוש חוזר בידית fn מההגדרה.

קנה מידה קטןקנה מידה גדול
קיוביטים1030
צעדי טרוטר1020
צעדים דחוסי-AQC (שכבה 1 + שכבה 2)3 + 2 = 56 + 4 = 10
שכבות אנסאץ מצב יסוד35
ממד קשר מקסימלי MPS32128
BackendstatevectorQPU עם DD, טוויירלינג פאולי, ו-TREX

שלב 1: מיפוי קלטים קלאסיים לבעיה קוונטית​

בנה את אותו SparsePauliOp הייזנברג KCuF3_3 והכן את מצב היסוד, הפעם עם אנסאץ gs_layers=5 עמוק יותר עבור השרשרת הארוכה יותר, ואז אפה את בעיטת הנייטרון π/2\pi/2 ZZ באתר המרכזי. זה זהה למיפוי בקנה המידה הקטן, אך ב-n=30n = 30.

יש לצפות לנאמנות מצב יסוד נמוכה יותר מאשר בהרצה של 10 אתרים: בסביבות 0.82 כאן לעומת 0.98 עבור השרשרת הקטנה יותר, מכיוון שחמש שכבות HVA אינן יכולות ללכוד באופן מלא מצב יסוד של 30 אתרים. זה צפוי ולא כישלון, והמדריך המקורי מקבל בערך 0.65 ב-50 אתרים מאותה הסיבה. הגדלת gs_layers או תקרת האיטרציות של COBYQA משפרת את זה, במחיר קלאסי נוסף.

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.

שלבים 2 ו-3: דחוס ובצע עם תבנית הפונקציה​

אותה קריאה בודדת כמו בדוגמת הסימולטור, הפעם עם backend_name שמצביע על מעבד IBM Quantum, כך שהפונקציה מתרגמת ומבצעת שם. תוכנית הדחיסה משנה את עומק האנסאץ: ששת צעדי הטרוטר הראשונים (בעלי שזירות נמוכה) נדחסים לתוך אנסאץ רדוד בעל שכבה אחת, ארבעת הבאים לתוך אנסאץ עמוק יותר בעל שתי שכבות, ועשרת הצעדים הנותרים מתוך 20 רצים כטרוטר רגיל. aqc_options מעלה את ממד הקשר של ה-MPS ל-max_bond=128 עבור השרשרת הארוכה והשזורה יותר (בהתאמה למקור), תוך שמירה על אותו אופטימייזר L-BFGS-B מוגבל ל-100 איטרציות. ה-estimator_options מפעילים את הפחתת השגיאות המובנית: ניתוק דינמי (XY4), טוויירלינג שערים, והפחתת מדידה TREX. ברירות המחדל של הפונקציה כבר תואמות למדריך המקורי בכל אלה מלבד תקציב הלמידה של TREX (measure_noise_learning). כל הבלוק עדיין כתוב במלואו מכיוון ש-estimator_options שמסופק על ידי הקורא מחליף את ברירות המחדל של הפונקציה במלואן במקום להתמזג לתוכן, כך שהשמטת מפתח הייתה מחזירה לברירת המחדל של IBM Quantum Compute במקום לזו של הפונקציה.

# 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
התחברות מחדש לעבודה ארוכת ריצה

ההרצה בקנה המידה הגדול אינה מהירה, ורוב הזמן הוא קלאסי ולא על ה-QPU. דחיסת ה-AQC רצה בתוך הפונקציה לפני שמשהו מגיע ל-QPU: ב-30 אתרים עם max_bond=128 זה לקח קרוב לארבע שעות בהרצה שלנו, לעומת כ-18 דקות של זמן QPU המצוינות באומדן השימוש בראש מדריך זה. זמן ההמתנה בתור מתווסף על שניהם. אינך צריך לשמור על המחברת או ה-kernel פתוחים בזמן שהיא רצה.

העתק את מזהה העבודה שהודפס על ידי התא הקודם ושמור אותו. שלושת התאים הבאים מאפשרים לך לחזור להרצה מאוחר יותר:

  1. התחבר מחדש, נדרש רק בסשן kernel חדש: הרץ מחדש את תאי ההגדרה כדי לשחזר את serverless, ואז בנה מחדש את ידית ה-job מהמזהה ששמרת. דלג על תא זה אם אתה עדיין באותו הסשן שבו הגשת, מכיוון שהידית כבר פעילה.

  2. בדוק סטטוס: הרץ מחדש עד שהוא מדווח DONE.

  3. אחזר את התוצאה: הרץ רק כאשר הסטטוס הוא DONE.

תא ההתחברות מחדש הבא מכיל מציין מקום. החלף אותו ב-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}

שלב 4: עבד ולאחר מכן החזר תוצאה בפורמט קלאסי רצוי​

עיבוד לאחר מכן זהה להרצת הסימולטור: בצע התמרת פורייה של פונקציית גרין לתוך S(q,ω)S(q, \omega), בצע סימטריזציה במראה, וקטום ערכים שליליים. עם השרשרת והאבולוציה הארוכות יותר, הרצף הדו-ספינוני מובחן טוב הרבה יותר. הוא צריך למלא את הרצועה בין הגבולות המקווקווים, בהיר ביותר קרוב ל-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

נספח​

דוגמת החומרה הקודמת מריצה אורך שרשרת בודד. שלושת הספקטרומים הבאים מגיעים מהרצות חומרה קודמות של אותה זרימת עבודה על ibm_pittsburgh ב-10, 20, ו-30 אתרים, כאשר כל קלט אחר נשמר קבוע: 20 צעדי טרוטר ב-dt = 0.6, תוכנית הדחיסה של שישה צעדים בעלי שכבה אחת בתוספת ארבעה צעדים דחוסי-AQC בעלי שתי שכבות, ו-max_bond = 128. אלה תוצאות מתועדות, לא פלט מהתאים הקודמים.

אותן הגדרות משמשות בכל שלושת הגדלים, כך שהספקטרומים ניתנים להשוואה ישירה. כוונון שלהם לפי אורך שרשרת, עם יותר שכבות אנסאץ מצב יסוד או max_bond גדול יותר, למשל, יכול לתת תוצאות טובות יותר מכל אלה שמוצגות כאן.

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

צעדים הבאים​

המלצות
מקור: IBM Quantum docs — עודכן בתאריך 9 בספט׳ 2026
גרסה באנגלית ב-doQumentation — עודכן בתאריך 9 בספט׳ 2026
בעיה באתר או בתרגום?בעיה בתוכן? ערוך ב-IBM Quantum docs