Diagonalizarea cuantică Krylov bazată pe eșantionare a unui model de rețea fermionică

Estimare de utilizare: Nouă secunde pe un procesor Heron r2 (NOTĂ: Aceasta este doar o estimare. Timpul tău de execuție poate varia.)

Fundal

Acest tutorial arată cum să folosești diagonalizarea cuantică bazată pe eșantionare (SQD) pentru a estima energia stării fundamentale a unui model de rețea fermionică. În mod concret, studiem modelul Anderson unidimensional cu o singură impuritate (SIAM), care este utilizat pentru a descrie impuritățile magnetice încorporate în metale.

Acest tutorial urmează un flux de lucru similar cu tutorialul înrudit Diagonalizarea cuantică bazată pe eșantionare a unui Hamiltonian chimic. Totuși, o diferență cheie constă în modul în care sunt construite circuitele cuantice. Celălalt tutorial folosește un ansatz variațional euristic, atrăgător pentru Hamiltonieni chimici cu potențial milioane de termeni de interacțiune. Pe de altă parte, acest tutorial folosește circuite care aproximează evoluția temporală guvernată de Hamiltonian. Astfel de circuite pot fi adânci, ceea ce face ca această abordare să fie mai potrivită pentru aplicații la modele de rețea. Vectorii de stare pregătiți de aceste circuite formează baza unui subspațiu Krylov, iar ca rezultat, algoritmul converge în mod demonstrabil și eficient spre starea fundamentală, în condiții potrivite.

Abordarea utilizată în acest tutorial poate fi privită ca o combinație a tehnicilor folosite în SQD și diagonalizarea cuantică Krylov (KQD). Abordarea combinată este uneori denumită diagonalizare cuantică Krylov bazată pe eșantionare (SQKD). Consultă Diagonalizarea cuantică Krylov a Hamiltonienilor de rețea pentru un tutorial despre metoda KQD.

Acest tutorial se bazează pe lucrarea "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", la care poți apela pentru mai multe detalii.

Modelul Anderson cu o singură impuritate (SIAM)

Hamiltonianul SIAM unidimensional este suma a trei termeni:

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

unde

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

Aici, $c^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma}$ sunt operatorii fermionici de creare/anihilare pentru al $\mathbf{j}^{\textrm{th}}$ sit de baie cu spin $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ sunt operatorii de creare/anihilare pentru modul de impuritate, iar $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$ și $V$ sunt numere reale care descriu interacțiunile de hopping, on-site și hibridizare, iar $\varepsilon$ este un număr real care specifică potențialul chimic.

Observă că Hamiltonianul este o instanță specifică a Hamiltonianului generic de electroni în interacțiune,

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

unde $H_1$ constă din termeni cu un corp, care sunt cuadratici în operatorii de creare și anihilare fermionici, iar $H_2$ constă din termeni cu două corpuri, care sunt quartici. Pentru SIAM,

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

iar $H_1$ conține restul termenilor din Hamiltonian. Pentru a reprezenta Hamiltonianul în mod programatic, stocăm matricea $h_{pq}$ și tensorul $h_{pqrs}$.

Bazele de poziție și impuls

Din cauza simetriei de translație aproximative din $H_\textrm{bath}$, nu ne așteptăm ca starea fundamentală să fie rară în baza de poziție (baza orbitalilor în care este specificat Hamiltonianul de mai sus). Performanța SQD este garantată doar dacă starea fundamentală este rară, adică are pondere semnificativă pe doar un număr mic de stări ale bazei de calcul. Pentru a îmbunătăți raritatea stării fundamentale, efectuăm simularea în baza orbitalilor în care $H_\textrm{bath}$ este diagonal. Numim această bază baza de impuls. Deoarece $H_\textrm{bath}$ este un Hamiltonian fermion cuadratic, el poate fi diagonalizat eficient printr-o rotație orbitală.

Evoluția temporală aproximativă guvernată de Hamiltonian

Pentru a aproxima evoluția temporală guvernată de Hamiltonian, folosim o descompunere Trotter-Suzuki de ordinul doi,

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

Sub transformarea Jordan-Wigner, evoluția temporală prin $H_2$ corespunde unei singure porți CPhase între orbitalii de spin-up și spin-down de la situl de impuritate. Deoarece $H_1$ este un Hamiltonian fermion cuadratic, evoluția temporală prin $H_1$ corespunde unei rotații orbitale.

Stările bazei Krylov $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, unde $D$ este dimensiunea subspațiului Krylov, sunt formate prin aplicarea repetată a unui singur pas Trotter, astfel

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

În fluxul de lucru bazat pe SQD de mai jos, vom eșantiona din acest set de circuite și vom post-procesa setul combinat de șiruri de biți cu SQD. Această abordare contrastează cu cea utilizată în tutorialul înrudit Diagonalizarea cuantică bazată pe eșantionare a unui Hamiltonian chimic, unde eșantioanele erau extrase dintr-un singur circuit variațional euristic.

Cerințe

Înainte de a începe acest tutorial, asigură-te că ai instalate următoarele:

• Qiskit SDK v1.0 sau mai recent, cu suport de vizualizare
• Qiskit Runtime v0.22 sau mai recent (pip install qiskit-ibm-runtime)
• SQD Qiskit addon v0.11 sau mai recent (pip install qiskit-addon-sqd)
• ffsim (pip install ffsim)

Pasul 1: Maparea problemei pe un Circuit cuantic

Mai întâi, generăm Hamiltonianul SIAM în baza de poziție. Hamiltonianul este reprezentat de matricea $h_{pq}$ și tensorul $h_{pqrs}$. Apoi îl rotim în baza de impuls. În baza de poziție, plasăm impuritatea la primul sit. Totuși, când rotim în baza de impuls, mutăm impuritatea la un sit central pentru a facilita interacțiunile cu alți orbitali.

# 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

În continuare, generăm circuitele pentru a produce stările bazei Krylov. Pentru fiecare specie de spin, starea inițială $\ket{\psi_0}$ este dată de superpoziția tuturor excitațiilor posibile ale celor trei electroni mai apropiați de nivelul Fermi în cele 4 moduri goale cele mai apropiate, pornind de la starea $|00\cdots 0011 \cdots 11\rangle$, realizată prin aplicarea a șapte porți XXPlusYYGate. Stările evoluate în timp sunt produse prin aplicări succesive ale unui pas Trotter de ordinul doi.

Pentru o descriere mai detaliată a acestui model și a modului în care sunt proiectate circuitele, consultă "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)

Pasul 2: Optimizează problema pentru execuție cuantică

Acum că am creat circuitele, le putem optimiza pentru un hardware țintă. Selectăm cel mai puțin ocupat QPU cu cel puțin 127 de qubiți. Consultă documentația Qiskit IBM® Runtime pentru mai multe informații.

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

Acum folosim Qiskit pentru a transpila circuitele către backend-ul țintă.

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)

Pasul 3: Execută folosind primitivele Qiskit

După optimizarea circuitelor pentru execuție hardware, suntem gata să le rulăm pe hardware-ul țintă și să colectăm eșantioane pentru estimarea energiei stării fundamentale. După ce folosim primitiva Sampler pentru a eșantiona șiruri de biți din fiecare circuit, combinăm toate rezultatele într-un singur dicționar de numărători și reprezentăm grafic cele mai frecvente 20 de șiruri de biți eșantionate.

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)

Pasul 4: Post-procesează și returnează rezultatul în formatul clasic dorit

Acum rulăm algoritmul SQD folosind funcția diagonalize_fermionic_hamiltonian. Consultă documentația API pentru explicații privind argumentele acestei funcții.

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

Următoarea celulă de cod reprezintă grafic rezultatele. Primul grafic arată energia calculată în funcție de numărul de iterații de recuperare a configurației, iar al doilea grafic arată ocuparea medie a fiecărui orbital spațial după iterația finală. Ca referință energetică, folosim rezultatele unui calcul DMRG efectuat separat.

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

Verificarea energiei

Energia returnată de SQD este garantată să fie o limită superioară a energiei adevărate a stării fundamentale. Valoarea energiei poate fi verificată deoarece SQD returnează și coeficienții vectorului de stare care aproximează starea fundamentală. Poți calcula energia din vectorul de stare folosind matricele de densitate redusă de ordinul 1 și 2, după cum este demonstrat în celula de cod de mai jos.

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

Referințe

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