Simulazione della diffusione dei neutroni con un flusso di lavoro serverless basato su AQC e dinamica di Trotter

Stima del tempo di esecuzione: 18 minuti su un processore Heron r3 (NOTA: si tratta solo di una stima. (La durata potrebbe variare.)

Risultati di apprendimento

• In che modo lo spettro di diffusione inelastica dei neutroni si ricollega al fattore di struttura dinamico $S(q, \omega)$ di un magnete quantistico di tipo “ 1D ”.
• Come preparare lo stato fondamentale KCuF$_3$ (Heisenberg isotropico) utilizzando il gruppo di rinormalizzazione della matrice di densità (DMRG) e la massimizzazione della fedeltà dello stato di prodotto matriciale (MPS).
• Come eseguire l'evoluzione temporale di Trotter, la compressione dei circuiti tramite compilazione quantistica approssimativa (AQC) e l'esecuzione mitigata con una singola chiamata di funzione.
• Come elaborare le serie temporali “ $\langle \sigma_z \rangle(t)$ ” per ciascun sito in “ $S(q, \omega)$ ” e identificare il continuum a due spinoni.

Prerequisiti

• Conoscenza dei modelli di Qiskit, SparsePauliOp, e dell’evoluzione temporale di Trotter.
• È utile, ma non indispensabile, avere una conoscenza di base dei metodi basati sulle reti di tensori (DMRG e MPS), così come è utile avere familiarità con la libreria qiskit-addon-aqc-tensor utilizzata dalla funzione per comprimere i circuiti di Trotter.

Sfondo

La diffusione anelastica dei neutroni misura il fattore di struttura dinamico $S(q, \omega)$, ovvero la trasformata di Fourier spazio-temporale della funzione di correlazione spin-spin; pertanto, la riproduzione dell’ $S(q, \omega)$ e a partire da un modello di spin microscopico costituisce un test diretto e falsificabile di una simulazione quantistica. Questo tutorial analizz KCuF$_3$, una catena di Heisenberg antiferromagnetica spin- $\frac{1}{2}$, le cui eccitazioni non sono singoli ribaltamenti di spin ma coppie di spinoni frazionati: invece di una dispersione netta dei magnoni, $S(q, \omega)$ mostra un ampio continuum a due spinoni, limitato inferiormente da $\tfrac{\pi}{2}|\sin q|$ e superiormente da $\pi|\sin(q/2)|$. Si tratta delle curve tratteggiate nei grafici che seguono. Gli aspetti fisici nel loro complesso e il confronto con i dati misurati relativi ai neutroni sono trattati nel tutorial originale e in Lee et al., arXiv:2603.15608.

Il flusso di lavoro quantistico rispecchia l'esperimento di diffusione:

1. Preparare lo stato fondamentale della catena $|\psi_0\rangle$.
2. Si inizia con una perturbazione locale nel sito centrale, una rotazione di tipo “ $\pi/2$ ” $Z$, che riproduce il trasferimento di quantità di moto ed energia da parte del neutrone.
3. Evoluzione temporale sotto l’Hamiltoniano di Heisenberg, $e^{-iHt}$, utilizzando una formula di prodotto di Trotter.
4. Misurare la magnetizzazione per sito $\langle \sigma_z^j \rangle(t)$. In funzione del sito $j$ e del tempo $t$, questa è esattamente la funzione di Green ritardata $G^R(j, j_c, t)$, quindi non è necessaria alcuna conversione prima della trasformata di Fourier nel passaggio 5.
5. Trasformata di Fourier: $G^R$ in $S(q, \omega)$.

Possono sorgere problemi nella fase 3, quando i circuiti di Trotter esatti per evoluzioni lunghe diventano troppo complessi per l'hardware. L'AQC con reti tensoriali risolve questo problema comprimendo un blocco di passi di Trotter in un ansatz parametrizzato, fisso e poco profondo, la cui fedeltà dello stato rispetto all'evoluzione esatta viene massimizzata in ambito classico con un simulatore MPS ( arXiv:2301.08609 ). Il template AQC Dynamics racchiude l'intero nucleo quantistico (sintesi di Trotter, compressione AQC ed esecuzione mitigata) in un'unica chiamata:

Le attività specifiche dell'esperimento rimangono qui nel quaderno: preparazione dello stato fondamentale (PRE) e post-elaborazione dell' $S(q, \omega)$ e (POST). Le due fasi ad alta intensità quantistica, la compressione e l'esecuzione, vengono eseguite all'interno della funzione.

Questo tutorial è un complemento al tutorial “Simulazione della diffusione neutronica nei materiali quantistici con circuiti quantistici”, che riproduce lo stesso esperimento passo dopo passo: lo stesso modello KCuF$_3$, la preparazione dello stato fondamentale, il “kick” neutronico e la post-elaborazione, con la sintesi di Trotter, la compressione AQC e l’esecuzione mitigata descritte passo dopo passo. Leggi quel tutorial per scoprire come funziona la compressione AQC. Leggi questo articolo per eseguire lo stesso esperimento tramite un modello di funzione distribuito: il nucleo quantistico diventa una singola chiamata di funzione e la compressione AQC, che richiede diverse ore, viene eseguita all’interno del worker serverless anziché sul tuo computer, quindi non è necessario disporre di un sistema HPC né di un kernel aperto durante l’esecuzione. La stessa chiamata guida anche altri esperimenti sulle dinamiche di 1D.

Requisiti

Prima di iniziare questo tutorial, assicurati di disporre di quanto segue:

• La funzione è stata implementata nel tuo account Qiskit Serverless. Esegui prima il modello di funzione di accompagnamento: distribuisci ed esegui il modello di funzione AQC + dinamica di Trotter. Questa guida illustra come procurarsi i file sorgente e caricare la funzione sul proprio account. Questo tutorial si limita a richiamare la funzione distribuita.

• IBM Quantum® credenziali salvate per QiskitServerless (vedere il modello di funzione). Entrambi gli esempi presenti in questo tutorial richiamano la funzione distribuita, quindi entrambi ne hanno bisogno.

• Qiskit SDK v2.0 o versioni successive (pip install qiskit).

• Il client del catalogo " IBM " di Qiskit (pip install qiskit-ibm-catalog).

• NumPy, SciPy, e Matplotlib (pip install numpy scipy matplotlib). Per l’ottimizzatore COBYQA utilizzato nella preparazione dello stato fondamentale è necessaria la versione 1.14 o successiva di SciPy.

• Lo stack della rete tensoriale AQC, poiché la preparazione dello stato fondamentale nella Fase 1 viene eseguita localmente in questo notebook: pip install 'qiskit-addon-aqc-tensor[quimb-jax]==0.3.1'.

La prima chiamata a una funzione appena distribuita rimane in attesa mentre il worker serverless installa le sue dipendenze, quindi è possibile che si verifichi una latenza aggiuntiva durante quell'esecuzione.

Configura

Importare le librerie e definire gli helper specifici dell’esperimento che verranno utilizzati in seguito: ( build_gs_ansatz l’approccio variazionale hamiltoniano, o HVA, per la preparazione dello stato fondamentale), prepare_ground_state (DMRG con massimizzazione della fedeltà MPS), e get_spectrum, plot_green, e plot_spectrum (la post-elaborazione “ $S(q, \omega)$ ”). Questi contenuti sono tratti dal tutorial originale sulla diffusione dei neutroni.

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.")

Output:

Setup complete - helpers defined.

Carica il modello di funzione

Accedi a Qiskit Serverless e carica il file aqc-dynamics-function. Entrambi gli esempi presenti in questo tutorial richiamano lo stesso fn handle, quindi la funzione viene caricata una sola volta, in questo caso.

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

Esempio di simulatore su piccola scala

Per prima cosa eseguiamo l’intero flusso di lavoro su una piccola catena di 10 sedi utilizzando lo stesso statevector backend. Questo permette di verificare la correttezza della pipeline PRE → FUNCTION → POST prima di impiegare risorse della QPU.

Fase 1: Mappare gli input classici su un problema quantistico

Costruire l' KCuF$_3$ come un Hamiltoniano ( SparsePauliOp Heisenberg isotropico: $XX + YY + ZZ$ con accoppiamento $\tfrac14$ su ciascun legame tra vicini più prossimi; le stringhe sono operatori di Pauli, quindi $\tfrac14$ fornisce l'accoppiamento di spin $\frac{1}{2}$ ). Preparare lo stato fondamentale con il metodo DMRG e la massimizzazione della fedeltà MPS, quindi integrare l’effetto del neutrone: una rotazione di tipo $\pi/2$ $Z$ nel sito centrale. Il circuito preparato è quello che passiamo alla funzione come initial_state. Lasciamo observables il valore predefinito ( $Z$ per sito), che corrisponde esattamente al valore $\langle \sigma_z^j \rangle(t)$ richiesto dal flusso di lavoro dei neutroni.

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}."
)

Output:

DMRG ground-state energy: -4.258035
GS fidelity: 0.9841
Prepared 10-qubit ground state with the neutron kick at site 4.

Passaggi 2 e 3: Comprimere ed eseguire con il modello di funzione

In un flusso di lavoro gestito manualmente, si tratta di due fasi distinte: l'ottimizzazione dei circuiti per l'hardware (Fase 2) e la loro esecuzione (Fase 3). Il modello di funzione raggruppa entrambe le operazioni in un’unica chiamata. Esegue la sintesi di Trotter, la compressione AQC e la transpilazione hardware, quindi esegue i circuiti (in questa fase sul simulatore esatto, successivamente con mitigazione degli errori integrata sull’hardware). I due parametri di configurazione sono aqc_segments (il piano di compressione) e aqc_options (le impostazioni dell'MPS e dell'ottimizzatore). Ogni segmento {"n_steps": k, "ansatz_steps": m} comprime i passi consecutivi k di Trotter in un ansatz costruito a partire da un obiettivo di Trotter di m-passi, mentre eventuali passi successivi vengono sum(n_steps) eseguiti come metodo di Trotter standard. I primi passi, caratterizzati da un basso livello di entanglement, si prestano bene alla compressione in un ansatz poco profondo (ansatz_steps=1), pertanto in questo caso comprimiamo i primi tre passi in un ansatz a singolo strato e i due successivi in un ansatz più profondo a due strati; i restanti cinque dei 10 passi di Trotter vengono eseguiti come Trotter standard. Infatti, aqc_options seguiamo alla lettera il tutorial originale: dimensione del legame MPS max_bond=32, cutoff=1e-8, e un ottimizzatore L-BFGS-B con un limite massimo di 100 iterazioni.

Chiama la funzione caricata in Setup. backend="statevector" esegue il percorso di riferimento esatto: nessun tempo QPU, con i circuiti in esecuzione su un simulatore di vettori di stato esatto all’interno del worker serverless (per richiamarlo è comunque necessario un account Qiskit Serverless salvato). Il contiene initial_state lo stato fondamentale preparato (compreso il kick); observables viene omesso in modo che la funzione misuri l' $Z$ predefinito per sito.

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

Output:

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)

Output:

AQC fidelities: {'1': 1.0, '2': 0.9999, '3': 0.9992, '4': 0.9998, '5': 0.9995}
Green's function shape: (10, 10)

Fase 4: Elaborazione successiva e restituzione del risultato nel formato classico desiderato

Applicare la trasformata di Fourier alla funzione di Green per ottenere l' $S(q, \omega)$ e, applicare la simmetria speculare ed eliminare i valori negativi: la post-elaborazione standard dei dati neutronici. Il mirroring è esatto poiché $S(q, \omega) = S(-q, \omega)$ per questo modello, e i valori negativi che permangono sono artefatti derivanti dalla trasformata di Fourier di una serie temporale finita e campionata in modo discreto, pertanto vengono troncati a zero. In questa piccola simulazione esatta, il continuum a due spinoni è risolto solo in modo approssimativo, ma il meccanismo è identico a quello della simulazione hardware che segue.

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:

Esempio di hardware su larga scala

Lo stesso flusso di lavoro è scalabile senza modificare il codice scientifico: una catena di 30 siti, il doppio della profondità di Trotter (20 passi), un piano di compressione che varia la profondità dell’ansatz (un ansatz più profondo per i passi successivi, caratterizzati da un maggiore entanglement) e l’esecuzione su un processore IBM Quantum con le tecniche di mitigazione degli errori integrate nella funzione (disaccoppiamento dinamico, rotazione di Pauli ed estinzione degli errori di lettura con rotazione (TREX)). Seguiamo gli stessi quattro passaggi dell'esempio del simulatore, riutilizzando l'handle fn creato in "Setup".

Fase 1: Mappare gli input classici su un problema quantistico

Costruisci lo stesso KCuF$_3$ Heisenberg e SparsePauliOp prepara lo stato fondamentale, questa volta con un gs_layers=5 ansatz più profondo per la catena più lunga, quindi applica il kick neutronico $\pi/2$ $Z$ nel sito centrale. È identico alla mappatura su piccola scala, ma all'indirizzo $n = 30$.

Ci si deve aspettare una fedeltà dello stato fondamentale inferiore rispetto alla simulazione a 10 siti: qui circa 0.82 contro 0.98 per la catena più piccola, poiché cinque strati HVA non sono in grado di rappresentare appieno uno stato fondamentale a 30 siti. Si tratta di un comportamento previsto, non di un errore, e il tutorial originale accetta circa 0.65 in 50 siti proprio per lo stesso motivo. Aumentare il limite di iterazioni del gs_layers COBYQA ne migliora le prestazioni, a un costo classico aggiuntivo.

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}."
)

Output:

DMRG ground-state energy: -13.111355
GS fidelity: 0.8201
Prepared 30-qubit ground state with the neutron kick at site 14.

Passaggi 2 e 3: Comprimere ed eseguire con il modello di funzione

La stessa singola chiamata dell'esempio del simulatore, ora con che backend_name punta a un processore IBM Quantum, in modo che la funzione venga transpilata ed eseguita lì. Il piano di compressione varia la profondità dell'ansatz: i primi sei passi di Trotter (a basso entanglement) vengono compressi in un ansatz a singolo strato poco profondo, i quattro successivi in un ansatz a due strati più profondo, mentre i restanti 10 dei 20 passi vengono eseguiti come Trotter standard. aqc_options aumenta la dimensione del legame MPS a max_bond=128 per la catena più lunga e più intrecciata (in linea con l'originale), mantenendo lo stesso ottimizzatore L-BFGS-B con un limite massimo di 100 iterazioni. estimator_options Attivare le funzionalità integrate di mitigazione degli errori: disaccoppiamento dinamico ( XY4 ), “gate twirling” e mitigazione delle misurazioni TREX. I valori predefiniti della funzione corrispondono già a quelli del tutorial originale per tutti questi elementi, ad eccezione del budget di apprendimento TREX (measure_noise_learning). L’intero blocco viene comunque riportato per intero poiché un oggetto fornito dal chiamante estimator_options sostituisce integralmente i valori predefiniti della funzione invece di fondersi con essi; pertanto, l’omissione di una chiave comporterebbe il ricorso al valore predefinito di IBM Quantum Compute anziché a quello della funzione.

# 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_SNSSW"]},
        "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)

Output:

job ID (save this to reconnect later): 43ed8d07-6d7d-4f33-b70a-7f31b765b310

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

Output:

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)

Output:

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}

Fase 4: Elaborazione successiva e restituzione del risultato nel formato classico desiderato

Post-elaborazione identica a quella della simulazione: applicare la trasformata di Fourier alla funzione di Green per ottenere un $S(q, \omega)$ e, applicare la simmetria speculare e tagliare i valori negativi. Grazie alla catena più lunga e all’evoluzione, il continuum a due spinoni risulta molto meglio risolto. Dovrebbe riempire la fascia compresa tra i limiti tratteggiati, con la massima luminosità in prossimità di $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:

Appendice

L'esempio hardware precedente utilizza una singola lunghezza di catena. I tre spettri riportati di seguito provengono da precedenti esecuzioni hardware dello stesso flusso di lavoro ibm_pittsburgh su 10, 20 e 30 siti, con tutte le altre variabili di input mantenute fisse: 20 passi di Trotter a dt = 0.6, il piano di compressione composto da sei passi a un livello più quattro passi a due livelli compressi con AQC, e max_bond = 128. Si tratta di risultati registrati, non dei valori delle celle precedenti.

Poiché per tutte e tre le dimensioni vengono utilizzate le stesse impostazioni, gli spettri sono direttamente comparabili. Regolarli in base alla lunghezza della catena, ad esempio aumentando il numero di livelli dell’ansatz dello stato fondamentale o utilizzando un valore maggiore max_bond, può fornire risultati migliori rispetto a quelli qui riportati.

Passi successivi