Stichprobenbasierte Krylov-Quantendiagonalisierung von an fermionischn Gittermodell

Nutzungsschätzung: Nein Sekunden auf an Heron r2 Prozessor (HINWEIS: Des is nur a Schätzung. Dei Laufzeit ko variirn.)

Hintergrund

Des Tutorial zoagt, wia mia de stichprobenbasierte Quantendiagonalisierung (SQD) verwenden, um de Grundzustandsenergie von an fermionischn Gittermodell z schätzen. Konkret untersuchen mia des eindimensionale Einzelstörstellen-Anderson-Modell (SIAM), des zur Beschreibung magnetischer Störstellen in Metallen verwendet wird.

Des Tutorial folgt an ähnlichn Arbeitsablauf wia des verwandte Tutorial Stichprobenbasierte Quantendiagonalisierung eines Chemie-Hamiltonians. A wesentlicher Unterschied liegt aber darin, wia de Quantenschaltkreise aufgebaut wern. Des andere Tutorial verwendet an heuristischen Variationsansatz, der für Chemie-Hamiltonians mit potenziell Millionen von Wechselwirkungstermen attraktiv is. Des Tutorial hingegen verwendet Schaltkreise, de de Zeitentwicklung durch den Hamiltonian approximiern. Solche Schaltkreise kinntn tief sei, was desn Ansatz besser für Anwendungen auf Gittermodelle macht. De Zustandsvektoren, de von diesen Schaltkreisen präpariert wern, bilden de Basis für an Krylov-Unterraum, und als Ergebnis konvergiert da Algorithmus nachweislich und effizient zum Grundzustand unter geeigneten Annahmen.

Da in dem Tutorial verwendete Ansatz ko als a Kombination von den in SQD und Krylov-Quantendiagonalisierung (KQD) verwendetn Techniken betrachtet wern. Da kombinierte Ansatz wird manchmal als stichprobenbasierte Krylov-Quantendiagonalisierung (SQKD) bezeichnet. Schau da Krylov-Quantendiagonalisierung von Gitter-Hamiltonians für a Tutorial zur KQD-Methode.

Des Tutorial basiert auf da Arbeit "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", auf de für weitere Details verwiesen wern ko.

Einzelstörstellen-Anderson-Modell (SIAM)

Da eindimensionale SIAM-Hamiltonian is a Summe aus drei Termen:

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

wobei

$$\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*}$$

Do san $c^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma}$ de fermionischen Erzeugungs-/Vernichtungsoperatoren für de $\mathbf{j}^{\textrm{te}}$ Bad-Stell mit Spin $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ san Erzeugungs-/Vernichtungsoperatoren für den Störstellenmodus, und $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$ und $V$ san reelle Zahlen, de de Hüpf-, Vor-Ort- und Hybridisierungswechselwirkungen beschreibn, und $\varepsilon$ is a reelle Zahl, de des chemische Potenzial spezifiziert.

Beacht, dass da Hamiltonian a spezifische Instanz vum generischen Wechselwirkungs-Elektronen-Hamiltonian is,

$$\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*}$$

wobei $H_1$ aus Einkörpertermen besteht, de quadratisch in den fermionischen Erzeugungs- und Vernichtungsoperatoren san, und $H_2$ aus Zweikörpertermen besteht, de quartisch san. Für des SIAM gilt

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

und $H_1$ enthält de restlichn Terme im Hamiltonian. Um den Hamiltonian programmatisch darzustellen, speichern mia de Matrix $h_{pq}$ und den Tensor $h_{pqrs}$.

Orts- und Impulsbasen

Wegen da approximativen Translationssymmetrie in $H_\textrm{bath}$ erwarten mia ned, dass da Grundzustand in da Ortsbasis (da Orbitalbasis, in da da Hamiltonian obn spezifiziert is) dünn besetzt is. De Leistung von SQD is nur garantiert, wenn da Grundzustand dünn besetzt is, des hoaßt, wenn er signifikantes Gewicht auf nur ana kleinen Anzahl von Rechenbasis-Zuständen ham. Um de Dünnbesetztheit vum Grundzustands z verbessern, führen mia de Simulation in da Orbitalbasis durch, in da $H_\textrm{bath}$ diagonal is. Mia nennen diese Basis de Impulsbasis. Da $H_\textrm{bath}$ a quadratischer fermionischer Hamiltonian is, ko er effizient durch a Orbitalrotation diagonalisiert wern.

Approximative Zeitentwicklung durch den Hamiltonian

Um de Zeitentwicklung durch den Hamiltonian z approximiern, verwenden mia a Trotter-Suzuki-Zerlegung zweiter Ordnung,

$$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}.$$

Unter da Jordan-Wigner-Transformation entspricht de Zeitentwicklung durch $H_2$ an einzelnen CPhase-Gate zwischn den Spin-up- und Spin-down-Orbitalen an da Störstell. Da $H_1$ a quadratischer fermionischer Hamiltonian is, entspricht de Zeitentwicklung durch $H_1$ ana Orbitalrotation.

De Krylov-Basiszustände $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, wobei $D$ de Dimension vum Krylov-Unterraums is, wern durch wiederholte Anwendung von an einzelnen Trotter-Schritt gebildet, sodass

$$|\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}.$$

Im folgndn SQD-basierten Arbeitsablauf wern mia Stichproben aus dem Satz von Schaltkreisen ziagn und den kombiniertn Satz von Bitfolgen mit SQD nachbearbeiten. Dieser Ansatz steht im Gegensatz zu dem im verwandtn Tutorial Stichprobenbasierte Quantendiagonalisierung eines Chemie-Hamiltonians verwendetn, wo Stichproben aus an einzelnen heuristischen Variationsschaltkreis gzogn worn san.

Anforderungen

Stell sicher, bevor du des Tutorial anfangst, dass du Folgendes installiert hast:

• Qiskit SDK v1.0 oder höher, mit Unterstützung für Visualisierung
• Qiskit Runtime v0.22 oder höher (pip install qiskit-ibm-runtime)
• SQD Qiskit Addon v0.11 oder höher (pip install qiskit-addon-sqd)
• ffsim (pip install ffsim)

Schritt 1: Problem auf an Quantenschaltkreis abbilden

Zuerst erzeugen mia den SIAM-Hamiltonian in da Ortsbasis. Da Hamiltonian wird durch de Matrix $h_{pq}$ und den Tensor $h_{pqrs}$ dargestellt. Danach rotieren mia ihn in de Impulsbasis. In da Ortsbasis platzieren mia de Störstell an da ersten Stell. Wenn mia aber zur Impulsbasis rotieren, verschieben mia de Störstell an a zentrale Stell, um Wechselwirkungen mit anderen Orbitalen z erleichtern.

# Added by doQumentation — installs packages not in the Binder environment
!pip install -q ffsim qiskit-addon-sqd

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

Als Nächstes erzeugen mia de Schaltkreise zur Erzeugung von de Krylov-Basiszuständen. Für jede Spin-Spezies is da Anfangszustand $\ket{\psi_0}$ durch de Superposition aller möglichn Anregungen von de drei Elektronen, de dem Fermi-Niveau am nächstn san, in de 4 nächstn leeren Moden gegeben, ausgehend vom Zustand $|00\cdots 0011 \cdots 11\rangle$, und wird durch de Anwendung von sieben XXPlusYYGates realisiert. De zeitentwickeltn Zustände wern durch sukzessive Anwendungen von an Trotter-Schritt zweiter Ordnung erzeugt.

Für a detailliertere Beschreibung von dem Modell und wia de Schaltkreise designt worn san, schau "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)

Schritt 2: Problem für de Quantenausführung optimiern

Nachdem mia de Schaltkreise erstellt ham, kinna mia sie für a Ziel-Hardware optimiern. Mia wählen de am wenigsten ausgelastete QPU mit mindestens 127 Qubits. Weitere Informationen findst in da Qiskit IBM® Runtime-Dokumentation.

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

Jetzt verwenden mia Qiskit, um de Schaltkreise für des Ziel-Backend z transpiliern.

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)

Schritt 3: Ausführung mit Qiskit-Primitiven

Nachdem mia de Schaltkreise für de Hardware-Ausführung optimiert ham, san mia bereit, sie auf da Ziel-Hardware auszuführn und Stichproben für de Grundzustandsenergieschätzung z sammeln. Nachdem mia de Sampler-Primitive verwendet ham, um Bitfolgen aus jedem Schaltkreis z ziagn, kombinieren mia alle Ergebnisse in an einzelnen Zähler-Wörterbuch und zeichnen de 20 am häufigsten gzognen Bitfolgen auf.

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)

Schritt 4: Nachbearbeitung und Rückgabe vum Ergebnis im gwünschtn klassischen Format

Jetzt führen mia den SQD-Algorithmus mit da Funktion diagonalize_fermionic_hamiltonian aus. Erläuterungen zu den Argumenten von derer Funktion findst in da API-Dokumentation.

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

De folgende Code-Zell zeichnet de Ergebnisse auf. De erste Grafik zeigt de berechnete Energie als Funktion von da Anzahl da Konfigurationswiederherstellungsiterationen, und de zweite Grafik zeigt de durchschnittliche Besetzung von jedem räumlichn Orbital nach da letztn Iteration. Für de Referenzenergie verwenden mia de Ergebnisse von ana DMRG-Berechnung, de separat durchgführt worn is.

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

Verifizierung von da Energie

De von SQD zurückgegebene Energie is garantiert a obere Schranke für de wahre Grundzustandsenergie. Da Wert von da Energie ko verifiziert wern, da SQD au de Koeffizienten vum Zustandsvektor zurückgibt, der den Grundzustand approximiert. Du konnst de Energie aus dem Zustandsvektor unter Verwendung von seine 1- und 2-Teilchen-reduzierten Dichtematrizen berechnen, wia in da folgndn Code-Zell demonstriert wird.

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

Referenzen

• 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)