Квантова діагоналізація Крилова на основі вибірок для ферміонної ґраткової моделі

Оцінка використання: дев'ять секунд на процесорі Heron r2 (ПРИМІТКА: Це лише оцінка. Ваш час виконання може відрізнятися.)

Передумови

Цей підручник демонструє, як використовувати квантову діагоналізацію на основі вибірок (SQD) для оцінки енергії основного стану ферміонної ґраткової моделі. Зокрема, ми досліджуємо одновимірну одноімпурітну модель Андерсона (SIAM), яка використовується для опису магнітних домішок, вбудованих у метали.

Цей підручник має подібний робочий процес до пов'язаного підручника Квантова діагоналізація на основі вибірок для хімічного гамільтоніана. Однак ключова відмінність полягає у способі побудови квантових схем. Інший підручник використовує евристичний варіаційний анзац, який є привабливим для хімічних гамільтоніанів із потенційно мільйонами членів взаємодії. З іншого боку, цей підручник використовує схеми, які апроксимують часову еволюцію за гамільтоніаном. Такі схеми можуть бути глибокими, що робить цей підхід кращим для застосування до ґраткових моделей. Вектори стану, підготовлені цими схемами, формують базис для підпростору Крилова, і завдяки цьому алгоритм доказово та ефективно збігається до основного стану за відповідних припущень.

Підхід, використаний у цьому підручнику, можна розглядати як комбінацію методів, що застосовуються в SQD та квантовій діагоналізації Крилова (KQD). Комбінований підхід іноді називають квантовою діагоналізацією Крилова на основі вибірок (SQKD). Дивіться Квантова діагоналізація Крилова для ґраткових гамільтоніанів для підручника з методу KQD.

Цей підручник базується на роботі "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", до якої можна звернутися за додатковими деталями.

Одноімпурітна модель Андерсона (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}$ є квадратичним ферміонним гамільтоніаном, його можна ефективно діагоналізувати за допомогою орбітального обертання.

Наближена часова еволюція за гамільтоніаном

Для наближення часової еволюції за гамільтоніаном ми використовуємо декомпозицію Троттера-Сузукі другого порядку,

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

За перетворенням Йордана-Віґнера часова еволюція за $H_2$ зводиться до одного вентиля CPhase між спін-вгору та спін-вниз орбіталями на вузлі домішки. Оскільки $H_1$ є квадратичним ферміонним гамільтоніаном, часова еволюція за $H_1$ зводиться до орбітального обертання.

Базисні стани Крилова $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, де $D$ — розмірність підпростору Крилова, формуються шляхом повторного застосування одного кроку Троттера, тому

$$|\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 або новіший, з підтримкою візуалізації
• Qiskit Runtime v0.22 або новіший (pip install qiskit-ibm-runtime)
• Додаток SQD для Qiskit 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

Далі ми генеруємо схеми для отримання базисних станів Крилова. Для кожного спінового різновиду початковий стан $\ket{\psi_0}$ задається суперпозицією всіх можливих збуджень трьох електронів, найближчих до рівня Фермі, у 4 найближчі порожні моди, починаючи зі стану $|00\cdots 0011 \cdots 11\rangle$, і реалізується застосуванням семи вентилів XXPlusYYGate. Стани з часовою еволюцією отримуються послідовним застосуванням кроку Троттера другого порядку.

Для більш детального опису цієї моделі та способу проєктування схем зверніться до "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 кубітами. Ознайомтеся з документацією 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 для транспіляції схем на цільовий бекенд.

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: Виконання за допомогою примітивів Qiskit

Після оптимізації схем для виконання на апаратному забезпеченні ми готові запустити їх на цільовому обладнанні та зібрати зразки для оцінки енергії основного стану. Після використання примітиву 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 також повертає коефіцієнти вектора стану, що апроксимує основний стан. Ти можеш обчислити енергію з вектора стану, використовуючи його одно- та двочастинкові зведені матриці густини, як продемонстровано в наступній комірці коду.

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)
• Chemistry Beyond Exact Solutions on a Quantum-Centric Supercomputer (стаття SQD)