Diagonalisasi kuantum Krylov berbasis sampel untuk model kisi fermionik

Estimasi penggunaan: Sembilan detik pada prosesor Heron r2 (CATATAN: Ini hanya estimasi. Waktu eksekusi kamu bisa berbeda.)

Hasil belajar

Setelah mengikuti tutorial ini, pengguna diharapkan memahami:

• Cara menggunakan SQD Qiskit addon untuk mendekati energi ground state dari model kisi menggunakan bitstring yang diambil dari unit pemrosesan kuantum (QPU).
• Cara menggunakan ffsim untuk membangun Circuit evolusi waktu untuk simulasi fermionik.
• Cara menggabungkan sampel dari beberapa sirkuit untuk pasca-pemrosesan dengan algoritma diagonalisasi Krylov berbasis sampel (SKQD).

Prasyarat

Kami menyarankan pengguna untuk familiar dengan topik-topik berikut sebelum mengikuti tutorial ini:

• Diagonalisasi kuantum berbasis sampel untuk Hamiltonian kimia
• Diagonalisasi kuantum Krylov untuk Hamiltonian kisi
• Qiskit primitives

Latar Belakang

Tutorial ini menunjukkan cara menggunakan diagonalisasi kuantum berbasis sampel (SQD) untuk memperkirakan energi keadaan dasar dari model kisi fermionik. Secara khusus, kita mempelajari model Anderson single-impurity (SIAM) satu dimensi, yang digunakan untuk mendeskripsikan impuritas magnetik yang tertanam dalam logam.

Tutorial ini mengikuti alur kerja yang mirip dengan tutorial terkait Diagonalisasi kuantum berbasis sampel untuk Hamiltonian kimia. Namun, perbedaan utamanya terletak pada bagaimana Circuit kuantum dibangun. Tutorial lain menggunakan ansatz variasional heuristik, yang menarik untuk Hamiltonian kimia dengan potensi jutaan suku interaksi. Di sisi lain, tutorial ini menggunakan Circuit yang mendekati evolusi waktu oleh Hamiltonian. Circuit semacam itu bisa sangat dalam, yang membuat pendekatan ini lebih baik untuk aplikasi pada model kisi. Vektor state yang dipersiapkan oleh Circuit-Circuit ini membentuk basis untuk subruang Krylov, dan hasilnya, algoritma terbukti dan secara efisien konvergen ke keadaan dasar, di bawah asumsi yang sesuai.

Pendekatan yang digunakan dalam tutorial ini dapat dipandang sebagai kombinasi teknik yang digunakan dalam SQD dan diagonalisasi kuantum Krylov (KQD). Pendekatan gabungan ini kadang-kadang disebut sebagai diagonalisasi kuantum Krylov berbasis sampel (SQKD). Lihat Diagonalisasi kuantum Krylov untuk Hamiltonian kisi untuk tutorial tentang metode KQD.

Tutorial ini didasarkan pada karya "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", yang dapat dirujuk untuk detail lebih lanjut.

Model Anderson single-impurity (SIAM)

Hamiltonian SIAM satu dimensi adalah jumlah dari tiga suku:

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

di mana

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

Di sini, $c^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma}$ adalah operator kreasi/anihilasi fermionik untuk situs bath ke-$\mathbf{j}^{\textrm{th}}$ dengan spin $\sigma$, $\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma}$ adalah operator kreasi/anihilasi untuk mode impuritas, dan $\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}$. $t$, $U$, dan $V$ adalah bilangan real yang mendeskripsikan interaksi hopping, on-site, dan hibridisasi, serta $\varepsilon$ adalah bilangan real yang menentukan potensial kimia.

Perhatikan bahwa Hamiltonian ini adalah instans spesifik dari Hamiltonian elektron interaksi generik,

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

di mana $H_1$ terdiri dari suku satu-badan, yang kuadratik dalam operator kreasi dan anihilasi fermionik, dan $H_2$ terdiri dari suku dua-badan, yang kuartik. Untuk SIAM,

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

dan $H_1$ memuat sisa suku-suku dalam Hamiltonian. Untuk merepresentasikan Hamiltonian secara terprogram, kita menyimpan matriks $h_{pq}$ dan tensor $h_{pqrs}$.

Basis posisi dan momentum

Karena simetri translasi aproksimasi dalam $H_\textrm{bath}$, kita tidak mengharapkan keadaan dasar menjadi jarang dalam basis posisi (basis orbital di mana Hamiltonian ditentukan di atas). Performa SQD hanya terjamin jika keadaan dasar bersifat jarang, artinya memiliki bobot signifikan hanya pada sejumlah kecil state basis komputasi. Untuk meningkatkan kejarangan keadaan dasar, kita melakukan simulasi dalam basis orbital di mana $H_\textrm{bath}$ diagonal. Kita sebut basis ini sebagai basis momentum. Karena $H_\textrm{bath}$ adalah Hamiltonian fermionik kuadratik, ia dapat didiagonalisasi secara efisien melalui rotasi orbital.

Aproksimasi evolusi waktu oleh Hamiltonian

Untuk mendekati evolusi waktu oleh Hamiltonian, kita menggunakan dekomposisi Trotter-Suzuki orde kedua,

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

Di bawah transformasi Jordan-Wigner, evolusi waktu oleh $H_2$ setara dengan satu Gate CPhase antara orbital spin-up dan spin-down pada situs impuritas. Karena $H_1$ adalah Hamiltonian fermionik kuadratik, evolusi waktu oleh $H_1$ setara dengan rotasi orbital.

State basis Krylov $\{ |\psi_k\rangle \}_{k=0}^{D-1}$, di mana $D$ adalah dimensi subruang Krylov, dibentuk oleh aplikasi berulang dari satu langkah Trotter, sehingga

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

Dalam alur kerja berbasis SQD berikut, kita akan mengambil sampel dari kumpulan Circuit ini dan memproses pasca-pemrosesan kumpulan bitstring yang digabungkan dengan SQD. Pendekatan ini berbeda dengan yang digunakan dalam tutorial terkait Diagonalisasi kuantum berbasis sampel untuk Hamiltonian kimia, di mana sampel diambil dari satu sirkuit variasional heuristik tunggal.

Persyaratan

Sebelum memulai tutorial ini, pastikan kamu telah menginstal hal-hal berikut:

• Qiskit SDK v1.0 atau lebih baru, dengan dukungan visualisasi
• Qiskit Runtime v0.22 atau lebih baru (pip install qiskit-ibm-runtime)
• SQD Qiskit addon v0.11 atau lebih baru (pip install qiskit-addon-sqd)
• ffsim v0.0.72 atau lebih baru (pip install ffsim)

Contoh simulator skala kecil

Langkah 1: Petakan masalah ke Circuit kuantum

Pertama, kita membuat Hamiltonian SIAM dalam basis posisi. Hamiltonian direpresentasikan oleh matriks $h_{pq}$ dan tensor $h_{pqrs}$. Kemudian, kita merotasinya ke basis momentum. Dalam basis posisi, kita menempatkan impuritas di situs pertama. Namun, ketika kita merotasi ke basis momentum, kita memindahkan impuritas ke situs tengah untuk memudahkan interaksi dengan orbital lain.

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

import numpy as np
import pyscf.fci

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 = 8

# 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

# Use PySCF to compute the exact ground state energy
reference_energy, _ = pyscf.fci.direct_spin1.kernel(h1e, h2e, norb, nelec)

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."""
    assert norb >= 8
    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())

Selanjutnya, kita membuat Circuit untuk menghasilkan state basis Krylov. Untuk setiap spesies spin, state awal $\ket{\psi_0}$ diberikan oleh superposisi dari semua kemungkinan eksitasi tiga elektron terdekat ke level Fermi ke dalam 4 mode kosong terdekat mulai dari state $|00\cdots 0011 \cdots 11\rangle$, dan diwujudkan melalui aplikasi tujuh XXPlusYYGate. State yang berevolusi terhadap waktu dihasilkan oleh aplikasi berturut-turut dari langkah Trotter orde kedua.

Untuk deskripsi lebih rinci tentang model ini dan bagaimana Circuit dirancang, lihat "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization".

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

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

from qiskit.providers.fake_provider import GenericBackendV2

backend = GenericBackendV2(
    2 * norb, basis_gates=["cp", "xx_plus_yy", "p", "x"]
)

Langkah 2: Optimalkan masalah untuk eksekusi kuantum

Selanjutnya, kita mengoptimalkan Circuit untuk perangkat keras target. Untuk saat ini, kita akan membuat Backend generik dengan jumlah qubit tertentu dan set Gate yang cocok untuk dekomposisi alami Circuit evolusi waktu.

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)

from qiskit.visualization import plot_histogram
from qiskit.primitives import StatevectorSampler

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

Sekarang, kita menggunakan Qiskit untuk men-transpile Circuit ke Backend target.

from qiskit.primitives import BitArray

# Combine the shots 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)

Langkah 3: Eksekusi menggunakan Qiskit primitives

Setelah mengoptimalkan Circuit untuk eksekusi hardware, kita siap menjalankannya di hardware target dan mengumpulkan sampel untuk estimasi energi ground state. Setelah menggunakan Sampler primitive untuk mengambil sampel bitstring dari setiap sirkuit, kita menggabungkan semua hasil ke dalam satu kamus counts dan memplot 20 bitstring yang paling sering diambil.

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: -13.4222953188441
		Subspace dimension: 529
	Subsample 1
		Energy: -13.42237556285828
		Subspace dimension: 784
	Subsample 2
		Energy: -13.422045397387413
		Subspace dimension: 529
Iteration 2
	Subsample 0
		Energy: -13.422379583305478
		Subspace dimension: 900
	Subsample 1
		Energy: -13.422376197704326
		Subspace dimension: 841
	Subsample 2
		Energy: -13.422421162849295
		Subspace dimension: 1089
Iteration 3
	Subsample 0
		Energy: -13.422421164670345
		Subspace dimension: 1156
	Subsample 1
		Energy: -13.422421492737689
		Subspace dimension: 1156
	Subsample 2
		Energy: -13.422421205869572
		Subspace dimension: 1156
Iteration 4
	Subsample 0
		Energy: -13.422421494558726
		Subspace dimension: 1225
	Subsample 1
		Energy: -13.422421492737689
		Subspace dimension: 1156
	Subsample 2
		Energy: -13.422421492737689
		Subspace dimension: 1156

Langkah 4: Pasca-pemrosesan dan kembalikan hasil ke format klasik yang diinginkan

Sekarang, kita menjalankan algoritma SQD menggunakan fungsi diagonalize_fermionic_hamiltonian. Lihat dokumentasi API untuk penjelasan argumen-argumen fungsi ini.

import matplotlib.pyplot as plt

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=reference_energy,
    color="#BF5700",
    linestyle="--",
    label="reference 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 energy: {reference_energy:.5f}")
print(f"SQD energy: {min_e:.5f}")
print(f"Absolute error: {abs(min_e - reference_energy):.5f}")
plt.tight_layout()
plt.show()

Reference energy: -13.42249
SQD energy: -13.42242
Absolute error: 0.00007

Sel kode berikut memplot hasilnya. Plot pertama menampilkan energi yang dihitung sebagai fungsi dari jumlah iterasi pemulihan konfigurasi, dan plot kedua menampilkan rata-rata hunian setiap orbital spasial setelah iterasi terakhir. Karena ini adalah masalah yang sangat kecil, iterasi pertama sudah membawa kita sangat dekat ke energi eksak (perhatikan skala sumbu y).

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: -13.42242

Verifikasi energi

Energi yang dikembalikan oleh SQD dijamin merupakan batas atas terhadap energi ground state yang sebenarnya. Nilai energi dapat diverifikasi karena SQD juga mengembalikan koefisien dari vektor state yang mendekati ground state. Kamu dapat menghitung energi dari vektor state menggunakan matriks densitas tereduksi satu- dan dua-partikel, seperti yang ditunjukkan pada sel kode berikut.

Contoh hardware skala besar

Sekarang, kita menjalankan contoh yang lebih besar pada QPU nyata. Sebagai energi referensi, kita menggunakan hasil perhitungan DMRG yang dilakukan secara terpisah.

from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import QiskitRuntimeService

# Model parameters
norb = 20
nelec = (norb // 2, norb // 2)
n_bath = norb - 1
hybridization = 1.0
hopping = 1.0
onsite = 10.0
chemical_potential = -0.5 * onsite

# Generate Hamiltonian and orbital rotation
h1e, h2e = siam_hamiltonian(
    norb=norb,
    hopping=hopping,
    onsite=onsite,
    hybridization=hybridization,
    chemical_potential=chemical_potential,
)
orbital_rotation = momentum_basis(norb)
h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())
impurity_index = n_bath // 2

# Set reference energy to DMRG value computed separately
reference_energy = -28.70659686

# Algorithm parameters
time_step = 0.2
krylov_dim = 8

# Construct circuits
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)
for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):
    circuit.append(instruction)
circuit.measure_all()
circuits = [circuit.copy()]
one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)
for i in range(krylov_dim - 1):
    circuit.remove_final_measurements()
    for instruction in trotter_step(
        qubits,
        time_step,
        one_body_evolution,
        h2e_momentum,
        impurity_index,
        norb,
    ):
        circuit.append(instruction)
    circuit.measure_all()
    circuits.append(circuit.copy())

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

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

# Sample from the circuits
sampler = Sampler(backend)
sampler.options.environment.job_tags = ["TUT_SKQD"]
job = sampler.run(isa_circuits, shots=500)

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

# Run configuration recovery and diagonalization
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,
)

# Plot results
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])
x1 = range(len(result_history))
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
axs[0].plot(x1, min_es, label="energy", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].axhline(
    y=reference_energy,
    color="#BF5700",
    linestyle="--",
    label="reference 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()
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 energy: {reference_energy:.5f}")
print(f"SQD energy: {min_e:.5f}")
print(f"Absolute error: {abs(min_e - reference_energy):.5f}")
plt.tight_layout()
plt.show()

Using backend ibm_boston
Iteration 1
	Subsample 0
		Energy: -28.63965951544449
		Subspace dimension: 9801
	Subsample 1
		Energy: -28.625588929202006
		Subspace dimension: 9409
	Subsample 2
		Energy: -28.647371834135498
		Subspace dimension: 8281
Iteration 2
	Subsample 0
		Energy: -28.67213260849567
		Subspace dimension: 29584
	Subsample 1
		Energy: -28.670340686158816
		Subspace dimension: 27225
	Subsample 2
		Energy: -28.669976379525988
		Subspace dimension: 31329
Iteration 3
	Subsample 0
		Energy: -28.68622875601382
		Subspace dimension: 36100
	Subsample 1
		Energy: -28.698569623143126
		Subspace dimension: 34225
	Subsample 2
		Energy: -28.694848533971882
		Subspace dimension: 33856
Iteration 4
	Subsample 0
		Energy: -28.69883392844593
		Subspace dimension: 42025
	Subsample 1
		Energy: -28.701289495200996
		Subspace dimension: 38025
	Subsample 2
		Energy: -28.699319594978245
		Subspace dimension: 45369
Iteration 5
	Subsample 0
		Energy: -28.701936886834154
		Subspace dimension: 51076
	Subsample 1
		Energy: -28.702468711812013
		Subspace dimension: 53824
	Subsample 2
		Energy: -28.702298147575938
		Subspace dimension: 52900
Reference energy: -28.70660
SQD energy: -28.70247
Absolute error: 0.00413

Langkah selanjutnya

Rekomendasi

Jika kamu menemukan karya ini menarik, kamu mungkin tertarik dengan materi berikut:

• Diagonalisasi kuantum berbasis sampel untuk Hamiltonian kimia - tutorial terkait yang menggunakan ansatz variasional heuristik alih-alih Circuit Trotter
• Diagonalisasi kuantum Krylov untuk Hamiltonian kisi - tutorial tentang metode KQD
• Dokumentasi API SQD addon - referensi untuk fungsi diagonalize_fermionic_hamiltonian
• Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization - makalah yang menjadi dasar tutorial ini