Sample-basierte Krylov-Quantendiagonalisiеrung von'm fermionischn Gittamodell

Nutzungsschätzung: Neun Sekunden uff'm Heron-r2-Prozessor (ACHTUNG: Dit is nur 'ne Schätzung. Deine Laufzeit kann abweichen.)

Hintajrund

Dit Tutorial zeigt dir, wie ma sample-basierte Quantendiagonalisierung (SQD) nutzt, um de Grundzustandsenergie von'm fermionischn Gittamodell zu schätzn. Konkret kiekn wa uns det eindimensionale Einzelstörstellen-Anderson-Modell (SIAM) an, det ma nutzt, um magnetische Störstellen in Metallen zu beschreibn.

Dit Tutorial folgt'm ähnlichn Workflow wie det verwandte Tutorial Sample-basierte Quantendiagonalisiеrung von'm Chemie-Hamiltonian. 'n wichtija Unjerschied liegt aba darin, wie de Quantenschaltkreise gebaut werdn. Det anda Tutorial nutzt'n heuristischn Variations-Ansatz, der for Chemie-Hamiltonians mit potenziell Millionen von Wechselwirkungstermen attraktiv is. Dit Tutorial hier nutzt dagegen Schaltkreise, die de Zeitentwicklung durch'n Hamiltonian approximieren. Solche Schaltkreise können tief sein, weshalb dissa Ansatz bessa für Anwendungen bei Gittamodellen jeeijnet is. De Zustandsvektoren, die diese Schaltkreise vorbereeten, bilden die Basis für'n Krylov-Unterraum, un als Folge davon konverjiert der Algorithmus nachweislich un effizient zum Grundzustand — unter jeeigneten Voraussetzungen.

Da Ansatz in dit Tutorial kann als'ne Kombination von Techniken aus SQD un Krylov-Quantendiagonalisierung (KQD) jesehn werdn. Die kombinierte Methode wird manchmal als sample-basierte Krylov-Quantendiagonalisiеrung (SQKD) bezeichnet. Kiek dir Krylov-Quantendiagonalisiеrung von Gitta-Hamiltonians for'n Tutorial zur KQD-Methode an.

Dit Tutorial basiert uff da Arbeit "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", wo du meha Details findn kannst.

Einzelstörstellen-Anderson-Modell (SIAM)

Da eindimensionale SIAM-Hamiltonian is'ne Summe von 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*}$$

Hiea sind $c^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma}$ die fermionischn Erzeugungs- und Vernichtungsoperatoren for'n $\mathbf{j}^{\textrm{ten}}$ Badplatz mit Spin $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ sind Erzeugungs- und Vernichtungsoperatoren for de Störstellenmode, un $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$ un $V$ sind reelle Zahlen, die de Hüpf-, Vor-Ort- un Hybridisierungswechselwirkungen beschreibn, un $\varepsilon$ is'ne reelle Zahl, die det chemische Potenzial anjibt.

Beachte, datt da Hamiltonian'n speziellen Fall vom allgemeinen Wechselwirkungs-Elektronen-Hamiltonian darstellt,

$$\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$ Einköpastermen enthält, die quadratisch in den fermionischn Erzeugungs- un Vernichtungsoperatoren sind, un $H_2$ Zweiköpasterme enthält, die quartisch sind. For's SIAM jilt:

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

un $H_1$ enthält de restlichn Terme im Hamiltonian. Um den Hamiltonian programmatisch darzustelln, speichern wa die Matrix $h_{pq}$ un den Tensor $h_{pqrs}$.

Orts- un Impulsraumbasis

Wegen da approximativen Translationssymmetrie in $H_\textrm{bath}$ erwarten wa nicht, datt da Grundzustand in da Ortsbasis (die Orbitalbasis, in der der Hamiltonian oben anjejebn is) dünn besetzt is. Die Leistung von SQD is nur jarantiert, wenn da Grundzustand dünn besetzt is, det heißt, er uff nur wenijen Berechnungsbasisnzuständen erhebliches Jewicht hat. Um de Dünnbesetzung vom Grundzustand zu verbessern, führn wa die Simulation in da Orbitalbasis durch, in der $H_\textrm{bath}$ diagonal is. Wa nennse diese Basis die Impulsraumbasis. Weil $H_\textrm{bath}$'n quadratischn fermionischn Hamiltonian is, kann er durch'ne Orbitalrotation effizient diagonalisiegt werdn.

Approximative Zeitentwicklung durch'n Hamiltonian

Um die Zeitentwicklung durch'n Hamiltonian zu approximieren, nutzn wa 'ne Trotter-Suzuki-Zerlegung zweita 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 der Jordan-Wigner-Transformation entspricht die Zeitentwicklung durch $H_2$ 'nem einzelnen CPhase-Gate zwischn den Spin-Aufwärts- un Spin-Abwärts-Orbitalen am Störstellenplatz. Da $H_1$ 'n quadratischn fermionischn Hamiltonian is, entspricht de Zeitentwicklung durch $H_1$ 'na Orbitalrotation.

Die Krylov-Basiszustände $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, wobei $D$ die Dimension vom Krylov-Unterraum is, werdn durch wiederholte Anwendung von 'nem einzelnen Trotter-Schritt jebildet, also:

$$|\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 folgenden SQD-basierten Workflow sampeln wa aus diesem Satz von Schaltkreisen un vaarbeiten die kombinierte Menge von Bitstrings mit SQD. Dissa Ansatz unjerscheidet sich von dem im verwandten Tutorial Sample-basierte Quantendiagonalisiеrung von'm Chemie-Hamiltonian, wo Samples aus 'nem einzelnen heuristischn Variations-Schaltkreis jezogn wordn sind.

Voraussetzungen

Bevor du mit dit Tutorial anfängst, stell sicher, datte foljendes installiert hast:

• Qiskit SDK v1.0 oda höha, mit Visualisierungs-Unterstützung
• Qiskit Runtime v0.22 oda höha (pip install qiskit-ibm-runtime)
• SQD Qiskit Addon v0.11 oda höha (pip install qiskit-addon-sqd)
• ffsim (pip install ffsim)

Schritt 1: Problem uff'n Quantenschaltkreis abbilden

Erst jeneriern wa den SIAM-Hamiltonian in da Ortsbasis. Der Hamiltonian wird durch die Matrix $h_{pq}$ un den Tensor $h_{pqrs}$ darstellt. Dann drehen wa ihn in die Impulsraumbasis. In da Ortsbasis plazieren wa die Störstelle am ersten Platz. Wenn wa aba in die Impulsraumbasis wechseln, verschieben wa die Störstelle an 'n mittlern Platz, um Wechselwirkungen mit andren Orbitalen zu 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 jeneriern wa die Schaltkreise, um die Krylov-Basiszustände zu erzeujen. For jede Spinspezies is da Anfangszustand $\ket{\psi_0}$ durch die Überlagerung aller möglichen Anregungen der drei Elektronen, die dem Fermi-Niveau am nächsten sind, in die 4 nächsten leeren Moden jejebn, ausgehend vom Zustand $|00\cdots 0011 \cdots 11\rangle$, un wird durch die Anwendung von sieben XXPlusYYGates realisiert. Die zeitentwickelten Zustände werdn durch sukzessive Anwendungen von 'm Trotter-Schritt zweita Ordnung erzeugt.

Fur'ne ausführlichere Beschreibung von diesem Modell un wie die Schaltkreise entworfen werdn, kiek dir "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization" an.

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 for die Quantenausführung optimiern

Jetzt, wo wa die Schaltkreise erstellt ham, können wa se for Zielhardware optimiern. Wa wähln den am wenigstn ausjelasteten QPU mit mindestns 127 Qubits. Kiek dir de Qiskit IBM® Runtime Docs an for meha Infos.

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 nutzn wa Qiskit, um die Schaltkreise for'n Ziel-Backend zu transpilieren.

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ührn mit Qiskit Primitives

Nachdem wa die Schaltkreise for die Hardwareausführung optimiert ham, sind wa parat, se uff da Zielhardware laufen zu lassn un Samples for die Grundzustandsenergieschätzung zu sammeln. Nachdem wa det Sampler-Primitive jenutzt ham, um Bitstrings aus jedem Schaltkreis zu sampeln, kombinieren wa alle Erjebnisse in 'n einzelnen Counts-Dictionary un plobn die 20 häufigsten sampelten Bitstrings.

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 un Rückgabe des Erjebnisses im jewünschten klassischn Format

Jetzt führn wa den SQD-Algorithmus mit da Funktion diagonalize_fermionic_hamiltonian aus. Kiek dir de API-Dokumentation an for Erklärungen zu den Argumenten dieser Funktion.

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

Die folgende Codezelle plobt die Erjebnisse. Da erste Plot zeigt die berechnete Energie als Funktion von da Anzahl der Konfijurationswiederherstellungs-Iterationen, un da zweite Plot zeigt die durchschnittliche Besetzung von jedem räumlichen Orbital nach da letztn Iteration. Als Referenzenergie nutzn wa die Erjebnisse von 'na DMRG-Berechnung, die separat durchjeführt wordn 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

Überprüfung der Energie

Die von SQD zurückjegebene Energie is jarantiert 'ne obere Schranke zur echtn Grundzustandsenergie. Da Wert von da Energie kann überprüft werdn, weil SQD auch die Koeffizienten vom Zustandsvektor zurückjibt, der den Grundzustand approximiert. Du kannst die Energie aus'm Zustandsvektor mit seinen 1- un 2-Teilchen-Dichtematrizen berechnen, wie in da folgenden Codezelle jezeigt.

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

Quelln

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