אלכסון קוונטי Krylov מבוסס דגימה של מודל סריג פרמיוני

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

רקע

מדריך זה מראה כיצד להשתמש באלכסון קוונטי מבוסס דגימה (SQD) כדי להעריך את אנרגיית מצב היסוד של מודל סריג פרמיוני. באופן ספציפי, אנו חוקרים את מודל Anderson חד-ממדי עם זיהום בודד (SIAM), המשמש לתיאור זיהומים מגנטיים המוטמעים במתכות.

מדריך זה עוקב אחר תהליך עבודה דומה למדריך הקשור אלכסון קוונטי מבוסס דגימה של המילטוניאן כימיה. עם זאת, הבדל מרכזי טמון באופן בו המעגלים הקוונטיים נבנים. המדריך האחר משתמש ב-ansatz וריאציוני היוריסטי, שהוא מושך עבור המילטוניאנים כימיים עם אולי מיליוני מונחי אינטראקציה. מצד שני, מדריך זה משתמש במעגלים המקרבים אבולוציית זמן על ידי ההמילטוניאן. מעגלים כאלה יכולים להיות עמוקים, מה שהופך גישה זו לטובה יותר ליישומים למודלי סריג. וקטורי המצב המוכנים על ידי מעגלים אלה מהווים את הבסיס ל-תת-מרחב Krylov, וכתוצאה מכך, האלגוריתם מתכנס באופן מוכח ויעיל למצב היסוד, בהנחות מתאימות.

את הגישה המשמשת במדריך זה ניתן לראות כשילוב של הטכניקות המשמשות ב-SQD ו-אלכסון קוונטי Krylov (KQD). הגישה המשולבת מכונה לעתים אלכסון קוונטי Krylov מבוסס דגימה (SQKD). ראה אלכסון קוונטי Krylov של המילטוניאנים סריג למדריך על שיטת KQD.

מדריך זה מבוסס על העבודה "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", שאליה ניתן לפנות לפרטים נוספים.

מודל Anderson עם זיהום בודד (SIAM)

ההמילטוניאן חד-ממדי SIAM הוא סכום של שלושה מונחים:

$$H = H_{\textrm{imp}}+ H_\textrm{bath} + H_\textrm{hyb},$$

כאשר

$$\begin{align*} H_\textrm{imp} &= \varepsilon \left( \hat{n}_{d\uparrow} + \hat{n}_{d\downarrow} \right) + U \hat{n}_{d\uparrow}\hat{n}_{d\downarrow}, \\ H_\textrm{bath} &= -t \sum_{\substack{\mathbf{j} = 0\\ \sigma\in \{\uparrow, \downarrow\}}}^{L-1} \left(\hat{c}^\dagger_{\mathbf{j}, \sigma}\hat{c}_{\mathbf{j}+1, \sigma} + \hat{c}^\dagger_{\mathbf{j}+1, \sigma}\hat{c}_{\mathbf{j}, \sigma} \right), \\ H_\textrm{hyb} &= V\sum_{\sigma \in \{\uparrow, \downarrow \}} \left(\hat{d}^\dagger_\sigma \hat{c}_{0, \sigma} + \hat{c}^\dagger_{0, \sigma} \hat{d}_{\sigma} \right). \end{align*}$$

כאן, $c^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma}$ הם אופרטורי היצירה/ההשמדה הפרמיוניים עבור אתר המרחץ ה-$\mathbf{j}^{\textrm{th}}$ עם ספין $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ הם אופרטורי יצירה/השמדה עבור מוד הזיהום, ו-$\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$, ו-$V$ הם מספרים ממשיים המתארים את האינטראקציות של דילוג, באתר והיברידיזציה, ו-$\varepsilon$ הוא מספר ממשי המציין את הפוטנציאל הכימי.

שים לב שההמילטוניאן הוא מופע ספציפי של ההמילטוניאן הגנרי של אלקטרון-אינטראקציה,

$$\begin{align*} H &= \sum_{\substack{p, q \\ \sigma}} h_{pq} \hat{a}^\dagger_{p\sigma} \hat{a}_{q\sigma} + \sum_{\substack{p, q, r, s \\ \sigma \tau}} \frac{h_{pqrs}}{2} \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma} \\ &= H_1 + H_2, \end{align*}$$

כאשר $H_1$ מורכב ממונחי גוף-אחד, שהם ריבועיים באופרטורי היצירה וההשמדה הפרמיוניים, ו-$H_2$ מורכב ממונחי שני-גופים, שהם רביעיים. עבור SIAM,

$$H_2 = U \hat{n}_{d\uparrow}\hat{n}_{d\downarrow}$$

ו-$H_1$ מכיל את שאר המונחים בהמילטוניאן. כדי לייצג את ההמילטוניאן באופן תוכניתי, אנו מאחסנים את המטריצה $h_{pq}$ ואת הטנזור $h_{pqrs}$.

בסיסי מיקום ותנע

בשל הסימטריה הטרנסלציונית המשוערת ב-$H_\textrm{bath}$, אנו לא מצפים שמצב היסוד יהיה דליל בבסיס המיקום (בסיס האורביטלים שבו ההמילטוניאן מוגדר למעלה). הביצועים של SQD מובטחים רק אם מצב היסוד דליל, כלומר, יש לו משקל משמעותי רק על מספר קטן של מצבי בסיס חישוביים. כדי לשפר את הדלילות של מצב היסוד, אנו מבצעים את הסימולציה בבסיס האורביטלים שבו $H_\textrm{bath}$ הוא אלכסוני. אנו קוראים לבסיס זה בסיס תנע. מכיוון ש-$H_\textrm{bath}$ הוא המילטוניאן פרמיוני ריבועי, ניתן לאלכסן אותו ביעילות על ידי סיבוב אורביטלי.

קירוב אבולוציית זמן על ידי ההמילטוניאן

כדי לקרב אבולוציית זמן על ידי ההמילטוניאן, אנו משתמשים בפירוק Trotter-Suzuki מסדר שני,

$$e^{-i \Delta t H} \approx e^{-i\frac{\Delta t}{2} H_2} e^{-i\Delta t H_1} e^{-i\frac{\Delta t}{2} H_2}.$$

תחת טרנספורמציית Jordan-Wigner, אבולוציית זמן על ידי $H_2$ מסתכמת ב-CPhase שער בודד בין האורביטלים ספין-למעלה וספין-למטה באתר הזיהום. מכיוון ש-$H_1$ הוא המילטוניאן פרמיוני ריבועי, אבולוציית זמן על ידי $H_1$ מסתכמת בסיבוב אורביטלי.

מצבי בסיס Krylov $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, כאשר $D$ הוא הממד של תת-מרחב Krylov, נוצרים על ידי יישום חוזר של צעד Trotter בודד, כך ש

$$|\psi_k\rangle \approx \left[e^{-i\frac{\Delta t}{2} H_2} e^{-i\Delta t H_1} e^{-i\frac{\Delta t}{2} H_2} \right]^k\ket{\psi_0}.$$

בתהליך העבודה המבוסס-SQD הבא, נדגום מקבוצת המעגלים הזו ונעבד לאחר מכן את קבוצת מחרוזות הביטים המשולבת עם SQD. גישה זו מנוגדת לזו המשמשת במדריך הקשור אלכסון קוונטי מבוסס דגימה של המילטוניאן כימיה, שבו דגימות נלקחו ממעגל וריאציוני היוריסטי בודד.

דרישות

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

• Qiskit SDK v1.0 או מאוחר יותר, עם תמיכת visualization
• Qiskit Runtime v0.22 או מאוחר יותר (pip install qiskit-ibm-runtime)
• SQD Qiskit addon v0.11 או מאוחר יותר (pip install qiskit-addon-sqd)
• ffsim (pip install ffsim)

שלב 1: מיפוי בעיה למעגל קוונטי

ראשית, אנו מייצרים את ההמילטוניאן SIAM בבסיס המיקום. ההמילטוניאן מיוצג על ידי המטריצה $h_{pq}$ והטנזור $h_{pqrs}$. לאחר מכן, אנו מסובבים אותו לבסיס התנע. בבסיס המיקום, אנו ממקמים את הזיהום באתר הראשון. עם זאת, כאשר אנו מסובבים לבסיס התנע, אנו מעבירים את הזיהום לאתר מרכזי כדי להקל על אינטראקציות עם אורביטלים אחרים.

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

import numpy as np

def siam_hamiltonian(
    norb: int,
    hopping: float,
    onsite: float,
    hybridization: float,
    chemical_potential: float,
) -> tuple[np.ndarray, np.ndarray]:
    """Hamiltonian for the single-impurity Anderson model."""
    # Place the impurity on the first site
    impurity_orb = 0

    # One body matrix elements in the "position" basis
    h1e = np.zeros((norb, norb))
    np.fill_diagonal(h1e[:, 1:], -hopping)
    np.fill_diagonal(h1e[1:, :], -hopping)
    h1e[impurity_orb, impurity_orb + 1] = -hybridization
    h1e[impurity_orb + 1, impurity_orb] = -hybridization
    h1e[impurity_orb, impurity_orb] = chemical_potential

    # Two body matrix elements in the "position" basis
    h2e = np.zeros((norb, norb, norb, norb))
    h2e[impurity_orb, impurity_orb, impurity_orb, impurity_orb] = onsite

    return h1e, h2e

def momentum_basis(norb: int) -> np.ndarray:
    """Get the orbital rotation to change from the position to the momentum basis."""
    n_bath = norb - 1

    # Orbital rotation that diagonalizes the bath (non-interacting system)
    hopping_matrix = np.zeros((n_bath, n_bath))
    np.fill_diagonal(hopping_matrix[:, 1:], -1)
    np.fill_diagonal(hopping_matrix[1:, :], -1)
    _, vecs = np.linalg.eigh(hopping_matrix)

    # Expand to include impurity
    orbital_rotation = np.zeros((norb, norb))
    # Impurity is on the first site
    orbital_rotation[0, 0] = 1
    orbital_rotation[1:, 1:] = vecs

    # Move the impurity to the center
    new_index = n_bath // 2
    perm = np.r_[1 : (new_index + 1), 0, (new_index + 1) : norb]
    orbital_rotation = orbital_rotation[:, perm]

    return orbital_rotation

def rotated(
    h1e: np.ndarray, h2e: np.ndarray, orbital_rotation: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
    """Rotate the orbital basis of a Hamiltonian."""
    h1e_rotated = np.einsum(
        "ab,Aa,Bb->AB",
        h1e,
        orbital_rotation,
        orbital_rotation.conj(),
        optimize="greedy",
    )
    h2e_rotated = np.einsum(
        "abcd,Aa,Bb,Cc,Dd->ABCD",
        h2e,
        orbital_rotation,
        orbital_rotation.conj(),
        orbital_rotation,
        orbital_rotation.conj(),
        optimize="greedy",
    )
    return h1e_rotated, h2e_rotated

# Total number of spatial orbitals, including the bath sites and the impurity
# This should be an even number
norb = 20

# System is half-filled
nelec = (norb // 2, norb // 2)
# One orbital is the impurity, the rest are bath sites
n_bath = norb - 1

# Hamiltonian parameters
hybridization = 1.0
hopping = 1.0
onsite = 10.0
chemical_potential = -0.5 * onsite

# Generate Hamiltonian in position basis
h1e, h2e = siam_hamiltonian(
    norb=norb,
    hopping=hopping,
    onsite=onsite,
    hybridization=hybridization,
    chemical_potential=chemical_potential,
)

# Rotate to momentum basis
orbital_rotation = momentum_basis(norb)
h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())
# In the momentum basis, the impurity is placed in the center
impurity_index = n_bath // 2

לאחר מכן, אנו מייצרים את המעגלים לייצור מצבי בסיס Krylov. עבור כל מין ספין, המצב ההתחלתי $\ket{\psi_0}$ ניתן על ידי סופרפוזיציה של כל העירורים האפשריים של שלושת האלקטרונים הקרובים ביותר לרמת Fermi לתוך 4 המודים הריקים הקרובים ביותר החל מהמצב $|00\cdots 0011 \cdots 11\rangle$, ומממש על ידי יישום שבעה XXPlusYYGates. המצבים שעברו אבולוציית זמן מיוצרים על ידי יישומים עוקבים של צעד Trotter מסדר שני.

לתיאור מפורט יותר של מודל זה וכיצד המעגלים מתוכננים, עיין ב-"Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization".

from typing import Sequence

import ffsim
import scipy
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit import CircuitInstruction, Qubit
from qiskit.circuit.library import CPhaseGate, XGate, XXPlusYYGate

def prepare_initial_state(qubits: Sequence[Qubit], norb: int, nocc: int):
    """Prepare initial state."""
    x_gate = XGate()
    rot = XXPlusYYGate(0.5 * np.pi, -0.5 * np.pi)
    for i in range(nocc):
        yield CircuitInstruction(x_gate, [qubits[i]])
        yield CircuitInstruction(x_gate, [qubits[norb + i]])
    for i in range(3):
        for j in range(nocc - i - 1, nocc + i, 2):
            yield CircuitInstruction(rot, [qubits[j], qubits[j + 1]])
            yield CircuitInstruction(
                rot, [qubits[norb + j], qubits[norb + j + 1]]
            )
    yield CircuitInstruction(rot, [qubits[j + 1], qubits[j + 2]])
    yield CircuitInstruction(
        rot, [qubits[norb + j + 1], qubits[norb + j + 2]]
    )

def trotter_step(
    qubits: Sequence[Qubit],
    time_step: float,
    one_body_evolution: np.ndarray,
    h2e: np.ndarray,
    impurity_index: int,
    norb: int,
):
    """A Trotter step."""
    # Assume the two-body interaction is just the on-site interaction of the impurity
    onsite = h2e[
        impurity_index, impurity_index, impurity_index, impurity_index
    ]
    # Two-body evolution for half the time
    yield CircuitInstruction(
        CPhaseGate(-0.5 * time_step * onsite),
        [qubits[impurity_index], qubits[norb + impurity_index]],
    )
    # One-body evolution for the full time
    yield CircuitInstruction(
        ffsim.qiskit.OrbitalRotationJW(norb, one_body_evolution), qubits
    )
    # Two-body evolution for half the time
    yield CircuitInstruction(
        CPhaseGate(-0.5 * time_step * onsite),
        [qubits[impurity_index], qubits[norb + impurity_index]],
    )

# Time step
time_step = 0.2
# Number of Krylov basis states
krylov_dim = 8

# Initialize circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# Generate initial state
for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):
    circuit.append(instruction)
circuit.measure_all()

# Create list of circuits, starting with the initial state circuit
circuits = [circuit.copy()]

# Add time evolution circuits to the list
one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)
for i in range(krylov_dim - 1):
    # Remove measurements
    circuit.remove_final_measurements()
    # Append another Trotter step
    for instruction in trotter_step(
        qubits,
        time_step,
        one_body_evolution,
        h2e_momentum,
        impurity_index,
        norb,
    ):
        circuit.append(instruction)
    # Measure qubits
    circuit.measure_all()
    # Add a copy of the circuit to the list
    circuits.append(circuit.copy())

circuits[0].draw("mpl", scale=0.4, fold=-1)

circuits[-1].draw("mpl", scale=0.4, fold=-1)

שלב 2: אופטימיזציה של הבעיה להרצה קוונטית

כעת, לאחר שיצרנו את המעגלים, אנו יכולים לבצע אופטימיזציה עליהם עבור חומרת יעד. אנו בוחרים את ה-QPU הכי פחות עסוק עם לפחות 127 qubit. בדוק את מסמכי Qiskit IBM® Runtime למידע נוסף.

from qiskit_ibm_runtime import QiskitRuntimeService

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

Using backend ibm_fez

כעת, אנו משתמשים ב-Qiskit כדי לבצע transpile של המעגלים ל-backend היעד.

from qiskit.transpiler import generate_preset_pass_manager

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend
)
isa_circuits = pass_manager.run(circuits)

שלב 3: ביצוע באמצעות primitives של Qiskit

לאחר אופטימיזציה של המעגלים להרצת חומרה, אנו מוכנים להריץ אותם על חומרת היעד ולאסוף דגימות להערכת אנרגיית מצב היסוד. לאחר שימוש ב-primitive Sampler לדגום מחרוזות ביטים מכל מעגל, אנו משלבים את כל התוצאות למילון ספירה אחד ומתווים את 20 מחרוזות הביטים הנדגמות הנפוצות ביותר.

from qiskit.visualization import plot_histogram
from qiskit_ibm_runtime import SamplerV2 as Sampler

# Sample from the circuits
sampler = Sampler(backend)
job = sampler.run(isa_circuits, shots=500)

from qiskit.primitives import BitArray

# Combine the counts from the individual Trotter circuits
bit_array = BitArray.concatenate_shots(
    [result.data.meas for result in job.result()]
)

plot_histogram(bit_array.get_counts(), number_to_keep=20)

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

כעת, אנו מריצים את אלגוריתם SQD באמצעות הפונקציה diagonalize_fermionic_hamiltonian. ראה את תיעוד ה-API להסברים על הארגומנטים לפונקציה זו.

from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
)

# 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}")
        print(
            f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )

rng = np.random.default_rng(24)
result = diagonalize_fermionic_hamiltonian(
    h1e_momentum,
    h2e_momentum,
    bit_array,
    samples_per_batch=100,
    norb=norb,
    nelec=nelec,
    num_batches=3,
    max_iterations=5,
    symmetrize_spin=True,
    callback=callback,
    seed=rng,
)

Iteration 1
	Subsample 0
		Energy: -28.61321893815165
		Subspace dimension: 10609
	Subsample 1
		Energy: -28.628985564542244
		Subspace dimension: 13924
	Subsample 2
		Energy: -28.620151775558114
		Subspace dimension: 10404
Iteration 2
	Subsample 0
		Energy: -28.656893066053115
		Subspace dimension: 34225
	Subsample 1
		Energy: -28.65277622004119
		Subspace dimension: 38416
	Subsample 2
		Energy: -28.670856034959165
		Subspace dimension: 39601
Iteration 3
	Subsample 0
		Energy: -28.684787675404362
		Subspace dimension: 42436
	Subsample 1
		Energy: -28.676984757118426
		Subspace dimension: 50176
	Subsample 2
		Energy: -28.671581704249885
		Subspace dimension: 40804
Iteration 4
	Subsample 0
		Energy: -28.6859683054753
		Subspace dimension: 47961
	Subsample 1
		Energy: -28.69418206537316
		Subspace dimension: 51529
	Subsample 2
		Energy: -28.686083516445752
		Subspace dimension: 51529
Iteration 5
	Subsample 0
		Energy: -28.694665630711178
		Subspace dimension: 50625
	Subsample 1
		Energy: -28.69505984237118
		Subspace dimension: 47524
	Subsample 2
		Energy: -28.6942873883992
		Subspace dimension: 48841

תא הקוד הבא משרטט את התוצאות. העלילה הראשונה מציגה את האנרגיה המחושבת כפונקציה של מספר איטרציות שחזור התצורה, והעלילה השנייה מציגה את תפוסת הממוצע של כל אורביטל מרחבי לאחר האיטרציה הסופית. עבור אנרגיית ההתייחסות, אנו משתמשים בתוצאות של חישוב DMRG שבוצע בנפרד.

import matplotlib.pyplot as plt

dmrg_energy = -28.70659686

min_es = [
    min(result, key=lambda res: res.energy).energy
    for result in result_history
]
min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])

# Data for energies plot
x1 = range(len(result_history))

# 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, min_es, label="energy", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].axhline(
    y=dmrg_energy, color="#BF5700", linestyle="--", label="DMRG energy"
)
axs[0].set_title("Approximated Ground State Energy vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy", 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})

print(f"Reference (DMRG) energy: {dmrg_energy:.5f}")
print(f"SQD energy: {min_e:.5f}")
print(f"Absolute error: {abs(min_e - dmrg_energy):.5f}")
plt.tight_layout()
plt.show()

Reference (DMRG) energy: -28.70660
SQD energy: -28.69506
Absolute error: 0.01154

אימות האנרגיה

האנרגיה המוחזרת על ידי SQD מובטחת להיות חסם עליון לאנרגיית מצב היסוד האמיתית. ניתן לאמת את הערך של האנרגיה מכיוון ש-SQD מחזיר גם את המקדמים של וקטור המצב המקרב את מצב היסוד. ניתן לחשב את האנרגיה מוקטור המצב באמצעות מטריצות הצפיפות המופחתות של חלקיק-1 ו-2, כפי שמודגם בתא הקוד הבא.

rdm1 = result.sci_state.rdm(rank=1, spin_summed=True)
rdm2 = result.sci_state.rdm(rank=2, spin_summed=True)

energy = np.sum(h1e_momentum * rdm1) + 0.5 * np.sum(h2e_momentum * rdm2)

print(f"Recomputed energy: {energy:.5f}")

Recomputed energy: -28.69506

אסמכתאות

• Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization (SKQD paper)
• Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer (SQD paper)
• Diagonalization of large many-body Hamiltonians on a quantum processor (KQD paper)